大青转录组测序及生物信息学分析

刘敏，黄炜忠，何孟璐，梅瑜，王继华

（1.广东省农作物遗传改良重点实验室/广东省农业科学院作物研究所，广东广州510640；2.广东罗浮山国药股份有限公司，广东惠州516100）

大青（Clerodendron cyrtophyllumTurez）是马鞭草科大青属（Verbenaceae clerodendronL.）灌木或小乔木植物，又名大青、大青木、山大青和羊咪青等，多生于海拔1 700米以下的平原、山地林下或溪谷旁，分布于两湖、两广及云贵高原等地，朝鲜、越南、马来西亚也有分布[1-2]。大青全株皆可入药，具有清热解毒、祛风除湿的功效，主治乙脑、流脑、外感热病热盛烦渴、咽喉肿痛、口疮、黄疸、热毒痢、急性肠炎、痈疽肿毒、衄血、血淋、外伤出血等病症。现代药理学研究表明，大青具有抑菌抗炎、抗病毒、利尿、镇痛降压的作用，是多种重要的中成药及保健品的主要成分[3-4]。大青叶在我国南方也常作为食疗药材，其干燥叶也是凉茶的主要原料之一，在广东多被制作成大青叶茶等保健产品[5]。

据《本草纲目》考证，民间药用大青为马鞭草科大青属植物大青，但其种类繁多。现根据地域区分就有广东大青、江西大青、广西大青、浙江大青等不同品种，这些不同品种之间的形态特征极其相似[6]。而现代中药材的资源鉴定多以形态鉴定及气味辨别为主，但不同大青的种类大多形态及气味极其相近，极难辨别鉴定区分；且不同品种之间虽然化学成分相近，但部分药效成分的含量却差别极大[7]。随着生物学技术的发展，基于高通量测序技术的转录组和生物信息学分析不受基因组信息的限制，在解析药用活性成分合成途径、挖掘相关功能基因、开发中药分子标记等方面得到了广泛的应用[8-12]。目前，许多药用植物的转录组学研究已经开展，如金银花[12]、地黄[13]、甘草[14]、茯苓[15]等，并鉴定到大量药效活性物质的合成途径相关基因，并开发了可靠的分子标记。目前，大青的研究多集中于化学成分鉴定及药理活性研究等方面，而关于大青有效化学成分的合成代谢途径关键基因的挖掘与鉴定尚未见报道。因此，本研究开展大青的转录组测序及生物信息学分析，以期为下一步解析其有效活性物质的合成代谢通路，挖掘其调控关键基因，开发分子标记提供基础数据，现将研究结果报道如下。

1 材料与方法

1.1 实验材料供试材料种植于广东省农业科学院作物研究所特色作物与南药资源圃。选取健康的大青植株的根、茎、叶，迅速用锡箔纸包裹并在液氮中冷冻，存储于-80℃冰箱中备用。

1.2 大青总RNA的提取用TRIzol（上海生工）方法分别提取大青根、茎、叶总RNA，采用Qubit 2.0 RNA检测试剂盒（Fluorometer Life Tech Invitrogen）分别对大青的根、茎、叶总RNA质量进行检测和定量，选取光密度OD（260 nm）/OD（230 nm）值大于1.8的RNA样品等量混合，建库。

1.3 大青转录组测序与拼接组装将检测合格的RNA委托北京百迈客生物科技有限公司采用Illumina HiSeqTM2500的高通量测序平台对大青的转录组进行测序。测序获得原始数据后，按照标准分析流程（通过FastQC软件进行质量评估，再使用Trimmomatic软件去掉接头和低质量序列）对数据进行处理，得到Clean Data用于de novo组装，然后采用Trinily软件将转录本拼接组装获得Unigene。

1.4 基因功能注释采用Blastx程序将组装的Unigene与保守域数据库（Conserved Domain Database，CDD），非冗余蛋白序列（NCBI non-redundant protein sequences，NR），核酸序列数据库（NCBI nonredundant nucleotide sequences，NT），蛋白结构域预测（Protein Families Database of Alignments and Hidden Markov Models，PFAM），真核生物蛋白质同源簇数据库/蛋白质聚类（eu Karyotic ortholog groups/clusters of orthologous groups，KOG/COG），SwissProt，TrEMBL，KEGG直系同源数据库（KEGGortholog database，KO），基因本体论（Gene Ontology，GO）等多个数据库比对，获得Unigene功能得到注释信息、GO功能注释信息。

1.5简单重复序列（SSR）、单核苷酸多态性（SNP）检测和引物设计采用微卫星识别工具（Microsatellite Identification Tool，MISA）软 件 对大青转录组所有Unigene存在的SSR微点进行鉴定及分析，并使用Primer 3软件（http://primer3.sourceforge.net/releases.php）设计SSR引物。

2 结果与分析

2.1 大青转录组测序与de novo组装采用Illumina HiSeqTM2500高通量测序平台进行大青转录组测序，获得Raw Reads后，通过软件去除接头盒低质量序列后，共得到26 394 223条Clean Reads，长度为7 859 604 312 bp，GC含量为45.5%，Q20达98.17%，Q30达94.58%。该结果表明，大青转录组数据质量较高，能够开展后续生物信息学分析。利用Trinity软件将大青Clean Reads进行de novo组装成转录本，共得到172 984条转录本，共计183 501 164 bp，平均长度1 060.8 bp，N50值为1 726 bp，其中序列长度大于500 bp的有105 429个，占总序列数的60.94%，序列长度大于1 000 bp的有64 345个，占总序列数的37.19%；去除冗余处理后得到100 191个Unigene，平均长度724.4 bp，N50值为1 055 bp，其中序列长度大于500 bp的Unigene数有42 108个，占总序列数的42.03%，序列长度大于1 000 bp的Unigene数有18 421个，占总序列数的18.39%。见表1、图1。从组装结果来看，转录本和Unigene的N50值均大于平均值，表明大青转录组的测序深度大，组装效果较好。

图1 大青Unigene序列长度分布Figure 1 Sequence length distribution of of unigenes of Clerodendron cyrtophyllum Turez

表1 大青转录组测序的结果Table 1 Results of transcriptome sequencing of Clerodendron cyrtophyllum Turez

2.2 大青转录组功能注释及其分类使用Blast软件将组装的Unigene与九大生物信息学数据库进行比对分析，得到大青转录组中的基因功能、代谢通路和转录因子等数据。通过比对分析，共有59 690（59.58%）个Unigene至少在一个数据库中获得信息注释，长度介于300 bp到1 000 bp之间的Unigene有31 333个，长度大于1 000 bp的Unigene有16 325个。除COG数据库外，在其他八大数据库中能注释的大青Unigene数目均在30%以上。其中：在NR和TrEMBL数据库中注释到的Unigene数目最多，分别达到58 952（58.84%）个和58 311（58.2%）个；在COG数据库中注释到的Unigene数目最少，仅为13 668（13.64%）个；49 278（49.18%）个Unigene在GO数据库中得到注释；38 260（38.19%）个Unigene在KEEG数据库中得到注释；32 453（32.39%）个Unigene在KOG数据库中得到注释；尚有40 501（40.42%）个Unigene尚未得到注释。见表2。结果表明，大青转录组获得的数据可以用于后续开展相应分析。

表2 大青Unigene注释的统计结果Table 2 Statistical results of annotated unigenes of Clerodendron cyrtophyllum Turez

2.2.1 NR功能注释大青的Unigene在NR数据库中的比对结果表明，其与芝麻（Sesamum indicum）、玫瑰木属（Rhodamnia argentea）、紫花风铃木（Handroanthus impetiginosus）的基因组有一定的相似度。其中与芝麻的匹配度最高，共有11 140条Unigene与之相匹配，占整个NR数据库注释Unigene的18.9%，其次分别为玫瑰木属和紫花风铃木，分别占比为18.41%和10.25%，见图2。从物种注释数据可以看出，大青的Unigene序列全部可以与植物相匹配，且与唇形目（管状花目）物种的序列相似度最高。但由于马鞭草科大青属植物的参考基因组数据比较少，因此并未从比对数据库中发现同科属的比对结果。

图2 大青Unigene的NR注释物种相似度分布Figure 2 Similarity distribution of NR annotated species of unigenes of Clerodendron cyrtophyllum Turez

2.2.2 KOG功能分类将组装的大青Unigene与KOG数据库比对分析，结果表明共有32 453条（32.39%）Unigene被注释到25个KOG功能分类中。其中涉及一般功能预测的基因最多，共有6 359个（占19.59%），其后依次是翻译后修饰，蛋白质转运，伴侣基因（3 945个，12.16%），信号转导机制基因（3 382个，10.42%），细胞内运输，分泌和囊泡运输（2 072个，6.38%），碳水化合物运输与代谢（2 110个，6.5%）及翻译，核糖体结构与生物起源（1 935个，5.96%），仅有11个基因注释到细胞运动性，占比0.03%。此外，还有1 862个Unigene被注释到未知功能。共有1 252个Unigene被注释到次生代谢物生物合成、运输和分解代谢，占总基因的比例为3.86%，这表明这些基因可能与大青中次生代谢物的生物合成、运输和积累相关。见图3。

图3 大青Unigene的KOG功能分类Figure 3 KOG functional classification of unigenes of Clerodendron cyrtophyllum Turez

2.2.3 GO功能注释将组装的大青Unigene与GO数据库比对分析，结果显示共有49 278个（49.18%）Unigene被注释到43个GO功能分类中。其中，生物过程注释到的Unigene最多，共有76 354个，占76.2%，其中注释较多的功能分别为细胞过程（26 146个，26.1%），代谢过程（22 970个，22.93%），生物调节（8 166个，8.15%），定位（5 359个，5.35%）及刺激响应（4 681个，4.67%）。注释到分子功能中的Unigene数目为55 108个，占比为55%，其中注释到较多功能的依次为结合（25 094个，25.05%），催化活性（22 737个，22.69%）及转运活性（2 509个，2.5%）。注释到细胞组分中的Unigene的数目为47 555个，占比为47.46%，其中注释到较多的功能依次为细胞结构物质（27 133个，27.08%），细胞内物质（15 645个，15.62%）及蛋白质复合物（4 771个，4.76%）。见图4。

图4 大青Unigene的GO功能分类Figure 4 GO functional classification of unigenes of Clerodendron cyrtophyllum Turez

2.2.4 KEGG代谢通路分析大青中主要含有萜类、黄酮类、酚酸类及氨基酸等药效活性成分，但大青缺少基因组信息作为参考，因此，对其基因功能注释具有较高的应用价值。将组装的大青Unigene与KEGG数据库比对分析，结果显示，共有38 260个Unigene得到注释，并参与到细胞过程、环境信息过程、生物体系统、人类疾病、遗传信息过程及代谢等六大类共136个代谢通路中。其中386个Unigene涉及苯丙烷生物合成（ko00940），165个Unigene涉及类黄酮生物合成（ko00941），37个Unigene参与异黄酮生物合成（ko00943），29个Unigene参与黄酮和黄酮醇生物合成（ko00944）；206个Unigene涉及萜类骨架生物合成（ko00900），130个Unigene涉及类胡萝卜素生物合成（ko00906），80个Unigene涉及玉米素生物合成（ko00908），50个Unigene参 与 单 萜 生 物 合 成（ko00902），71个Unigene参 与 二 萜 生 物 合 成（ko00904），80个Unigene参与倍半萜和三萜生物合成（ko00909）。这些与次生代谢产物相关的Unigene的鉴定结果为进一步解析大青药用成分物质的生物合成打下了基础。见图5。

图5 大青Unigene的KEGG功能分类Figure 5 KEGG functional classification of unigenes of Clerodendron cyrtophyllum Turez

2.2.5 转录因子分析植物转录因子主要调节各种生物过程的基因表达模式，且普遍存在于植物中，如植物中常见的有bHLH、MYB和NAC等。通过研究这些TFs的作用机制是功能基因组学的重要环节。根据大青的对比结果显示，共有3 543个Unigene被分为210个转录因子家族，其中C2H2转录因子类的Unigene数量最多，达到了115个，占比为3.24%，其次是bHLH、RLK-Pelle_DLSV、MYB-related、GRAS、C3H等。这些转录因子涉及到大青的生长发育、抗逆、次生代谢和合成等多种生物学途径。大青转录因子的分析为进一步研究其基因调控与基因互作提供了一定的数据支持，也为提高其类黄酮、萜类酚酸类等有效成分的生物合成和抗逆性提供了理论依据。见图6。

图6 大青转录因子分类Figure 6 Classification of transcription factors of Clerodendron cyrtophyllum Turez

2.3 大青转录组中遗传标记点位的分布SSR是检测遗传多样性和构建遗传图谱的有效分子标记之一。利用MISA软件从大青Unigene序列中鉴定潜在的SSR标记，并且进行统计。结果表明，在6 680条Unigene上共检测到8 640个SSR位点。其中，400条Unigene中检测到513个（5.8%）复杂重复类型的SSR位点。最丰富的重复类型是单碱基重复，共检测到5 135个位点，占59.4%；其次为双碱基重复（2 195个，25.4%），三碱基重复（745个，8.6%），混合碱基重复（500个，5.7%）和四碱基重复（42个，0.4%）；最少的为五碱基重复和六碱基重复，分别仅检测到5个和5个位点，分别占0.06%和0.06%。根据8 640个SSR的位点和Unigene序列，使用Premier 3.0设计出相应的SSR扩增引物，为采用分子标记辅助育种技术的大青遗传改良、群体结构的研究提供了基础数据。见图7。

图7 大青Unigene的SSR位点分析Figure 7 SSR locus analysis of unigenes of Clerodendron cyrtophyllum Turez

3 讨论

近年来，本草基因组学、转录组学技术的快速发展使得更多药用植物的基因资源得以被保护及利用。通过转录组数据的分析研究，可挖掘重要植物功能基因和代谢通路，进一步构建中药指纹图谱，为天然药物来源新途径、种质资源鉴定、保存及良种选育提供分子基础[16]；其次，通过对代谢通路关键基因的研究，可为中药资源活性成分的生物合成与调控提供新的思路和方法；第三，通过转录水平的调节，提高药用成分的产量与活性，寻找最佳中药状态，为重要的良种选育、规范种植、质量控制提供技术支撑[17]。

目前，对大青化学成分的分析和植株药理作用的研究已有一定的经验，但对大青的研究多限于整体水平实验，还未见关于其基因组、转录组学的研究报道，这制约了大青次生代谢产物的生物合成途径解析，不能充分挖掘其新的药用价值以使其资源得以充分利用[18-19]。本研究采用Illumina HiSeq 2500的高通量测序平台对大青的不同组织进行转录组测序并构建大青的Unigene库。通过de novo组装共获得100 191条Unigene，N50为1 055 bp，长度大于1 000 bp的Unigene有18 421个，占总序列数的18.39%。N50的长度是评价转录组组装质量的重要指标，大青转录组组装的N50长度与茶树（1 081 bp）的长度相近，但比中药黄芩（797.64 bp）要长，说明大青的转录组组装序列质量较高，能够满足后续数据分析的要求[20-21]。本研究通过与多个公共数据库进行比对分析获得Unigene的功能注释，结果表明，大青组装的Unigene共有59 690个（59.58%），Unigene在至少1个数据库中得到注释，其中以NR数据库中比对到的Unigene最多，共有58 952个Unigene得到注释，其结果显示与管状花目物种芝麻的相似度最高。这可能是由于马鞭草科大青属植物的种类少、研究基础薄弱导致大青转录组Unigene并未比对到同科属物种。另外，还有40 501个（40.42%）Unigene没有在一个数据库中获得注释，可能是由于组装的Unigene序列太短而缺乏保守区域，或者数据库中缺乏大青基因组、转录组信息导致未能匹配得到注释。在KEGG数据库中共注释到38 260个Unigene参与136个代谢途径，包括萜类、类胡萝卜素、生物碱类、黄酮类等次级代谢的生物合成途径。其中，407个Unigene参与萜类化合物的生物合成，231个Unigene参与黄酮类化合物的生物合成，386个Unigene参与苯丙烷素的生物合成。这些数据为进一步解析大青中活性成分的生物合成途径提供了研究基础。SSR标记操作简单、重复性好。转录组数据中包含着大量SSR位点，随着高通量测序技术的快速发展，利用这些位点开发的功能ESTSSR标记在群体遗传学及标记辅助育种研究中显得更加方便快捷。随着大量药用植物的转录组数据产生及功能基因得到注释，也必然为药用植物的功能EST-SSR标记开发奠定了基础[22]。本研究利用MISA软件共检测到8 640个SSR位点，其中单碱基重复最为丰富（5 135个，59.4%），有很大一部分SSR位点的Unigene在各大数据库中得到功能注释，这也为大青后续功能EST-SSR分子标记开发提供了数据基础。

本研究对大青转录组进行了初步的探究，弥补了大青基因组信息的不足，为解析大青次级代谢物质合成通路及分子生物学方面的研究奠定了基础。大青SSR位点的发掘，可为其分子标记的开发、遗传多样性分析、种质资源鉴定与优选、分子标记育种等提供理论基础，为利用分子手段鉴定、区分大青品种及评价其质量提供了依据。