石好宇，王慧芸，赵艺，朱甫臻，尹艳燕，缪剑华，王灵芝*（.北京中医药大学生命科学学院，北京 0009；.广西药用植物园，广西药用植物资源保护与遗传改良重点实验室，南宁 53003）

中药鸦胆子是苦木科植物鸦胆子[

Brucea javanica

（L.）Merr.]的干燥成熟果实，又名鸦蛋子、老鸦胆、苦榛子、苦参子，主产于我国广东、广西、云南、福建、海南和台湾等地，始记于清代《本草纲目拾遗》。鸦胆子味苦、性寒，有小毒，归大肠、肝经，有清热解毒、截疟止痢等功效。鸦胆子含有多种化学成分，主要包括苦木素类、类固醇、三萜类、生物碱、木质素、黄酮类、甾体类和脂肪酸等。苦木素类化合物是鸦胆子主要的药理活性物质，是一类由苦楝烷型三萜降解得到的高度氧合三萜及其苷类化合物，多为四环三萜及五环三萜内酯，是鸦胆子的代表性药理活性成分，其中鸦胆子苦醇（brusatol）和鸦胆丁（bruceantin）是该药材中代表性抗肿瘤成分；苦木素还具有抗炎和抗医学原虫等药理活性。鸦胆子油是鸦胆子中的脂肪油，抗肿瘤疗效确切且毒副作用小，被广泛应用于各种肿瘤的辅助治疗。鸦胆子中蛋白质含量为17.47%，鸦胆子多肽具有抗菌活性，本课题组前期研究发现鸦胆子球蛋白酶解物具有显著的抑制肿瘤细胞增殖的活性。

目前关于鸦胆子的研究主要集中于新化学成分的发现及作用，对于其次生代谢物合成途径的研究较少，缺乏鸦胆子基因和蛋白质序列信息，极大地限制了对鸦胆子生长发育研究以及功能分析和利用。转录组测序（RNA-seq）是利用高通量测序技术将细胞或组织中mRNA、small RNA或no-codingRNA 进行测序分析的技术，可用于发现新基因、提供基因表达和调控信息。该技术已广泛用于药用植物，如半夏、黄芪以及薏苡仁等研究中。因此本文对鸦胆子进行转录组测序，并分析苦木素类化合物萜类骨架生物合成以及其衍生化修饰途径，为今后药用成分开发和药用植物的遗传改良提供理论支撑。

1 材料

1.1 药材

鸦胆子采集于广西南宁市武鸣区（108.27719°N，23.15643°E），采集不同发育时期的根、茎、叶、花序及果实等组织部位，液氮速冻后，保存于－80℃备用。

1.2 仪器与试药

Aglient RNA 6000 Nano Reagents 试剂盒、Agilent 2100 生物分析仪（美国Aglient 公司）；BGISEQ-500 测序仪（深圳华大智造科技有限公司）；NanoDrop 超微量分光光度计（美国Nanodrop 公司）。

2 方法

2.1 RNA 提取与转录组文库构建

利用Aglient RNA 6000 试剂盒提取鸦胆子各组织总RNA。采用NanoDrop 超微量分光光度计进行样本RNA 纯度检测，通过Agilent 2100生物分析仪进行样本浓度及完整性RIN（RNA integrity number）进行检测。质检合格后，等量混合样品，用带有Oligo（dT）的磁珠富集mRNA，采用fragmentation buffer 将mRNA 片段化，随后使用随机N6 引物进行反转录合成双链cDNA，然后将双链DNA 黏性末端补平、5'端磷酸化、添加并连接测序接头，经PCR 扩增后构建cDNA 文库，质检合格的文库使用BGISEQ-500进行测序。

2.2 转录组测序分析及功能注释

测序原始序列（raw reads）使用SOAPunke软件进行过滤，去除接头污染、未知碱基N 含量大于5%、低质量的reads（质量低于10 的碱基占总碱基数20%以上），获得clean reads。由于鸦胆子基因组测序尚未完成，缺少参考序列，因此使用Trinity 软件对clean reads 进行

de novo

组装，将组装获得的转录本去除冗余，得到非重复序列基因（Unigene）用于后续分析。使用单拷贝直系同源数据库BUSCO（benchmarking universal single-copy orthologs）对组装的转录本进行质量评估。

组装得到转录组数据后，使用Blast 软件将Unigene 序列与NR（non-redundant protein sequence database）、NT（nucleotide sequence database）、Swissprot、KEGG（kyoto encyclopedia of genes and genomes）和KOG（clusters of orthologous groups for eukaryotic complete genomes）数据库进行比对，进行功能注释。基于NR 注释，使用Blast2GO 软件对Unigene 进行GO（gene ontology）注释，使用Hmmscan 软件进行Pfam 注释。

2.3 生物信息学分析

Getorf 软件检测Unigene 的开放阅读框（ORF）后使用Hmmsearch 将ORF 比对到转录因子蛋白结构域（数据来源于TF），然后根据PlantTFDB 描述的转录因子家族特征对Unigene进行能力鉴定。使用TransDecoder 软件识别Unigene 中的候选编码区域，通过Blast 比对SwissProt 数据库和Hmmscan 搜索Pfam 蛋白同源序列，从而进行编码序列（CDS）预测。经MISA软件对Unigene 进行分析后使用Primer3 对检测到的简单重复序列（SSR）进行引物设计。

2.4 基因表达量计算

基因表达量用FPKM（fragments per kilobase of transcript per million mapped reads）表示，即每百万reads 中来自比对到某一基因每千碱基长度的reads 数目。本实验使用Bowtie2 软件将clean reads 比对到基因序列上，然后使用RSEM 计算各个样品的基因表达水平。计算公式如下：

FPKM ＝（cDNA fragments）/[mapped fragments（millions）/transcript length（kd）]

式中，cDNA fragments 表示比对到某一转录本上的片段数目，即双端reads 数目；mapped fragments（Millions）表示比对到转录本上的片段总数，以10为单位；transcript length（kd）表示转录本的长度，以10个碱基为单位。

3 结果与分析

3.1 转录组测序与组装

采用NanoDrop 超微量分光光度计和Agilent 2100 生物分析仪对各样本RNA 进行质量检测，结果RIN ＞5.8，基线平整，表明完整性较好，样品质量合格。使用BGISEQ-500 平台共获得鸦胆子114.55 Mb Raw reads，过滤去除低质量、接头污染以及未知碱基N 含量过高的reads后，获得110.67 Mb（96.61%）clean reads，其中质量值大于20（Q20）的碱基数目占总碱基数目的98.73%，且无明显AT 或GC 分离现象（见图1A），表明测序质量较高，能满足后续分析。

对clean reads 进行组装获得转录本，使用BUSCO 软件进行质量评估，其完全匹配的BUSCOs占比98.02%，表明组装效果良好。使用Tgicl 软件对转录本进行聚类并去除冗余得到85 828 条Unigenes，平均长度为1775 bp，GC 含量为40.07%，N50 值为2608 bp，N70 值为1887 bp，N90 值1005 bp。对Unigenes 的长度进行统计（见图1B），其中长度≥3000 bp 的序列数目最多，为14 663 条。

图1 Clean reads 的碱基含量（A）与Unigenes 长度（B）分布图Fig 1 Distribution of clean reads base content（A）and Unigenes length（B）

3.2 鸦胆子Unigenes 的功能注释

为获得更全面的基因功能信息，将Unigenes与NR、NT、Swissprot、KEGG、KOG、Pfam和GO 七大功能数据库进行比对，共有68 806 条Unigenes 得到注释，占总数的80.17%。如表1 所示，有67 087（78.16%）条Unigenes 在NR 数据库匹配成功，数量最多，而在GO 数据库中，获得注释的Unigenes 数量最少，仅有43 042 个，占50.15%；有25 734 条Unigens 在七大功能注释库同时得到注释，占Unigenes 总数的29.98%。

表1 Unigenes 功能注释表
Tab 1 Function annotation of Unigenes

数据库数量占比/%NR67 08778.16 NT59 85069.73 Swissprot49 80158.02 KEGG52 41661.07 KOG53 75562.63 Pfam 52 21460.84 GO43 04250.15 Intersection25 73429.98

3.2.1 NR 数据库对比分析 将Unigenes 与NR 数据库进行相似序列匹配，结果如图2A 所示，与鸦胆子Unigenes 匹配数最多的物种是甜橙（

Citrus sinensis

），约占47.07%。此外，克莱门柚（

Citrus clementina

）和可可树（

Theobroma cacao

）同源性分别占27.7%和2.08%。

3.2.2 GO 数据库功能分类注释 将鸦胆子Unigenes 与GO 数据库进行比对，共有43 042 条Unigenes 得到注释，获得了99 940 条对应关系，结果如图2B 所示，分属于生物学过程、细胞组分和分子功能三大类。在生物学过程分类中，与细胞过程（cellular process）相关的Unigenes 为11 377 条，所占比例最高；在细胞组分分类中，与膜部分（membrane part）相关的Unigenes 最多，为12 940 条；在分子功能分类中，注释上的Unigenes 涉及抗氧化、结合、催化、分子功能调节、转运、转录调节等多方面功能活性，其中，涉及到结合（binding）功能的Unigenes 有21 235条，所占比例最高；对应蛋白标签（protein tag）功能的Unigenes 最少，为7 条。

3.2.3 KEGG 代谢通路分析 鸦胆子52 416 条Unigenes 被注释到KEGG 数据库的134 条代谢通路中（见图2C），其中涉及包含细胞过程（cellular processes，2282 条）、环境信息加工（environmental information processing，3191 条）、遗传信息加工（genetic information processing，10 076条）、新陈代谢（metabolism，29 789 条）及生物体系统（organismal systems，1775 条）五大类。其中，与新陈代谢相关的Unigenes 数量最多，占总数的63.23%，主要参与影响氨基酸代谢、糖代谢、脂代谢等多个方面。

对鸦胆子药用活性成分合成相关信号通路进行了分析，涉及不饱和脂肪酸生物合成（Ko01040）、萜类骨架生物合成（Ko00900）、类固醇生物合成（Ko00100）相关Unigenes 数量分别为206、276、8 条。根据鸦胆子KEGG 注释结果，筛选276 条基因参与萜类骨架生物合成，其中涉及甲羟戊酸（MVA）以及非甲羟戊酸（MEP）途径的分别有68、83 条Unigenes（见表2）。萜类化合物的生物合成涉及甲羟戊酸途径和脱氧木酮糖-5-磷酸途径两条生物合成途径，羟甲基戊二酰 CoA 还原酶（HMGR）是萜类代谢过程中的关键调控位点，有8 条Unigenes 调控此位点，其FPKM 值为85.37。1-脱氧木糖-5-磷酸合酶（DXS）与1-脱氧木糖-5-磷酸还原异构酶（DXR）是MEP 途径关键限速酶，分别有53 和4条基因编码，其FPKM 值分别为121.74、64.03。

表2 鸦胆子萜类骨架生物合成中编码关键酶的Unigenes
Tab 2 Unigenes encoding the key enzymes related with terpenoid backbone biosynthesis in (L.) Merr

合成途径酶的名称酶编号Unigenes 数目FPKM 值MVAacetyl-CoA acetyltransferase（ACAT）2.3.1.96124.71 MVAhydroxy methylglutaryl-CoA synthase（HMGS）2.3.3.10540.20 MVAhydroxymethylglutaryl-CoA reductase（HMGR）1.1.1.34885.37 MVAmevalonate kinase（MVK）2.7.1.36217.02 MVAphosphomevalonate kinase（PMK）2.7.4.2311.56 MVAmevalonate diphosphate decarboxylase（MVD）4.1.1.33333.39 MVAisopentenyl-diphosphate delta-isomerase（IPPs）5.3.3.211339.72 MVAfarnesyl diphosphate synthase（FDPS）2.5.1.130239.09 MEP1-deoxy-D-xylulose-5-phosphate synthase（DXS）2.2.1.753121.74 MEP1-deoxy-D-xylulose-5-phosphate reductoisomerase（DXR）1.1.1.267464.03 MEP2-C-methyl-D-erythritol 4-phosphate cytidylyltransferase（ispD）2.7.7.60820.36 MEP4-diphosphocytidyl-2-C-methyl-D-erythritol kinase（ispE）2.7.1.148548.69 MEP2-C-methyl-D-erythritol 2，4-cyclodiphosphate synthase（ispF）4.6.121120.13 MEP（E）-4-hydroxy-3-methylbut-2-enyl-diphosphate synthase（ispG）1.17.7.1188.64 MEP4-hydroxy-3-methylbut-2-en-1-yl diphosphate reductase（ispH）1.17.7.411203.29

萜类骨架形成后需要衍生修饰，细胞色素P450（cytochrome P450，CYP450）和糖基转移酶（UDP-glycosyltransferase，UGT）主要起氧化、置换和糖基化修饰的重要作用。根据鸦胆子转录组SwissProt 数据库注释结果，有284 条 Unigenes被注释为 CYP450，隶属于20 个CYP450 家族，其中 CYP71 家族的Unigenes 为71 条，占比最多，其次是CYP82 和CYP704，分别为34 条和22 条。共找到属于18 个UGT 亚家族的129 个UGTs，包括22 个 UGT74、17 个 UGT83 等。对各个样品中基因表达量进行聚类表达分析，其中CYP450家族 FPKM 值较高的有 CYP82A1（FPKM ＝284.88）、CYP71A1（FPKM ＝112.44）和CYP71A2（FPKM ＝105.64），UGT 家族FPKM 值较高的有UGT73C3（FPKM ＝65.93）、UGT71A15（FPKM ＝53.73）和UGT71A16（FPKM ＝36.44）。

3.2.4 KOG 分类 利用Blast 软件，共有53 755条Unigenes 被注释到KOG 数据库中，得到基因同源物的分类信息（见图2D）。KOG 根据其功能可分为25 类，其中，一般功能预测（general function prediction only）类别中Unigenes 数量最多，占总量的21.94%；细胞运动（cell motility）类Unigenes 数量最少， 为1 个； 未知功能（function unknown），翻译后修饰、蛋白质转换、分子伴侣（posttranslational modification，protein turnover，chaperones）、信号转导机制（signal transduction mechanisms）和转录（transcription）等类别的Unigenes 表达丰度也相对较高，分别占总量的9.32%、7.92%、10.78%和6.50%。

图2 Unigene 功能注释图Fig 2 Functional annotation of Unigenes

3.3 CDS 预测

CDS 预测可为鸦胆子基因组图谱绘制及基因功能研究提供关键依据，结果共预测得到56 624个CDS 序列，总长度为64 922 145 bp，N50 值为1446 bp，GC 含量为43.48%。课题组前期工作表明，鸦胆子球蛋白酶解物具有显著的体外抗肿瘤活性，通过注释，共获得10 个与球蛋白相关的CDS序列，结果见表3，可为今后高球蛋白鸦胆子药材的筛选及抗肿瘤生物肽的开发提供文本信息。

表3 鸦胆子球蛋白相关CDS 序列
Tab 3 Globulin related CDS of (L.) Merr

编号基因IDFPKM 值功能注释数据库来源1 Unigene 1078627.63Basic 7S globulinSwissprot 2 Unigene 123901.69Basic 7S globulinSwissprot 3 Unigene 1523759.38Basic 7S globulinSwissprot 4 Unigene 1564050.42Basic 7S globulinSwissprot 5 Unigene 20499.46Basic 7S globulinSwissprot 6 Unigene 2744610.94Basic 7S globulin 2Swissprot 7 Unigene 718641.50Basic 7S globulinSwissprot 8 Unigene1305220.96Citrus sinensis 12S seed storage globulin 1-likeNT 9 Unigene 108290.632S globulin [Corchorus olitorius]NR 10Unigene 308863.03basic 7S globulin-like [Manihot esculenta]NR

3.4 鸦胆子转录组SSR 分析

鸦胆子转录组共鉴定出6 种SSR 重复类型，共检测出36 045 个SSR 分布于25 747 个Unigenes中，单核苷酸至六核苷酸重复类型均存在，其中有14 631 个单核苷酸，11 594 个二核苷酸，7554个三核苷酸，649 个四核苷酸，574 个五核苷酸，1043 个六核苷酸，单核苷酸重复类型数量最多。

4 讨论

RNA-seq 技术能够对真核生物复杂的转录本结构和表达水平进行分析，还可用于发现非编码RNA 和全新的基因转录本，有助于揭示生物生长发育、逆境应答机制和次生代谢产物的富集调控机制。此外，该技术还在低丰度基因检测、基因家族鉴定等方面发挥重要作用。有学者利用Illumina 高通量测序发现三江源地区灌木亚菊中与药用活性相关的基因主要富集在甾体、黄酮类、苯丙素类及萜类化合物；通过与乌头转录组测序分析比较，证实铁棒锤是乌头碱生物合成途径的理想材料；梅瑜等利用高通量测序PacBio Sequel 平台深入了解甘葛藤转录组的整体水平，挖掘其中黄酮类生物合成相关候选基因；林江波等利用Illumina HiSeq 4000 高通量测序技术对铁皮石斛茎和叶进行转录组测序，并分析了植物甾醇生物合成关键酶基因的表达。

本研究首次采用BGISEQ-500 测序平台，进行了鸦胆子转录组序列分析，共获得85 828 条Unigenes，数据量大且组装效果良好，填补了鸦胆子转录组信息的空白，确定了参与萜类骨架生物合成的基因，文本信息将有助于进一步研究苦木素类化合物及其衍生化在鸦胆子中生物合成的分子机制，并促进鸦胆子功能基因组学的研究。MEP 和MVA 通路在植物萜类化合物的生物合成中发挥着重要作用。在次级代谢物的衍生化修饰过程中，需要CYP450 和糖基转移酶的参与。Zheng 等研究表明CYP450 参与了真菌甲基的羟基化反应；UGT73AD1 可使羟基积雪草酸的C-28 羧基部分糖基化。本研究发现鸦胆子转录组中有较多CYP450、UGT 序列，推测其可能在苦木素类衍生修饰中起重要作用。

本研究共获得56 624 个CDS 序列，构建鸦胆子基因编码数据库，丰富了该物种蛋白质编码基因的注释，为鸦胆子及其近缘物种蛋白序列鉴定提供了更可靠的质谱数据库信息。后期可通过蛋白质组学等手段来进行数据校正，并确定新的编码基因；与蛋白质组学联用还可以有效地量化生物体中重要的蛋白质，确定基因在植物生长发育过程中的功能。鸦胆子球蛋白及其酶解物具有明确的抗肿瘤活性，在获得的CDS 序列中有10 个与球蛋白相关，今后可借助分子标记辅助育种技术，进行高球蛋白亲本的筛选，增强鸦胆子球蛋白抗肿瘤活性。通过植物转录组测序，还可反向预测和开发多肽类药物。Rodríguez-Decuadro 等根据

Peltophorum dubium

幼苗转录组测序结果，采用同源序列搜索和半胱氨酸残基模式匹配的方法，对其中的抗菌肽进行了预测。

鸦胆子SSR 位点重复类型以单核苷酸和双核苷酸为主，占72.7%，与苹果、川穹等植物研究结果一致。短重复序列的大量存在表明该物种的进化水平相对较高，表明鸦胆子处于相对较高的生物进化与分类地位，其基因组经历了较长的进化时间或具有较高的突变频率，还可为近缘物种SSR 标记的开发及其遗传分析提供便利。

5 结论

本研究基于高通量测序平台进行了鸦胆子全长转录组分析，获得了高质量转录本，共有68 806 条Unigenes 注释到七大功能数据。萜类骨架是苦木素生物合成的重要前体，KEGG 分析共有276 条Unigenes 参与，其中涉及MEP 和MVA途径的15 种酶；还发现涉及鸦萜类骨架修饰的CYP450 和UGT 的Unigenes 分别有284 条和129条，为提高苦木素类化合物产量的分子生药学研究提供了重要文本信息。通过CDS 预测共获得10 个CDS 序列与球蛋白相关，进一步丰富了鸦胆子球蛋白源抗肿瘤肽的数据库。本研究获得了鸦胆子转录组的数据，可为今后鸦胆子药用价值开发与利用提供理论基础。