基于RNA-seq技术的楮头红转录组分析

金红焦根林陈刚

（1. 深圳市中国科学院仙湖植物园，深圳　518004；2. 肇庆学院生命科学学院，肇庆　526061）

基于RNA-seq技术的楮头红转录组分析

金红1焦根林1陈刚2

（1. 深圳市中国科学院仙湖植物园，深圳518004；2. 肇庆学院生命科学学院，肇庆526061）

利用RNA-seq技术对所构建的楮头红叶片的转录组进行测定，对原始reads进行过滤和组装，得到了51 305条质量较高的Unigenes，平均长度为921 nt，N50为1 490 nt。利用BLAST和BLAST2GO 软件对这些从头组装的Unigenes进行注释。用NCBI蛋白质数据库（Nr）、非冗余核苷酸数据库（Nt）、基因本体论（GO）、直系同源基因簇（COG）和京都基因与基因组百科全书（KEGG）数据库做参考，共注释了40 532条Unigenes。注释到Nr、Nt、Swiss-Prot、KEGG、COG和GO库中的比例相对较高，分别为77.53%、56.18%、53.14%、46.58%、29.69%和60.72%。在蛋白质数据库中对所有的Unigenes进行blast以后，发现有39 302个CDS，用ESTscan预测了2 065个CDS。KEGG通路分析显示，参与次生代谢物生物合成的Unigenes有2 323条，占全部Unigenes的9.72%。其中有78条Unigenes编码了细胞色素P450家族蛋白，这些信息为药用植物次生代谢物生物合成关键基因的挖掘提供了理论参考。

楮头红转录组；RNA-seq

楮头红（Sarcopyramis nepalensis Wall.），民间又称风柜斗草，为野牡丹科肉穗草属植物，无亚种名。主要分布于我国台湾、福建、江西、云南和四川等地。该植物全草可入药，其提取物中主要成分为黄酮类、多酚、脂肪酸和酚酸类［1，2］，常被用于治疗急、慢性肝炎、肺热咳嗽、风湿骨痛、无名肿毒等疾病［3］，是较有价值的珍稀名贵中草药，在我国台湾和福建地区广为使用。然而，目前对楮头红的研究主要停留在药理活性及化学成分等方面［4-6］，缺乏对其分子生物学方面的研究，其基因信息也极少，给楮头红次生代谢途径的研究带来了极大困难。

转录组是在特定的发育阶段和不同的生理条件下，所有转录出来的RNA的集合。转录组学研究是解读基因组功能元件、揭示细胞和组织中的分子组成所必需的［7］。目前由于蛋白质组学技术的限制，转录组学成为研究基因表达调控的主要方法之一，是连接基因组学和蛋白质组学的桥梁。通过高通量的转录组分析，可以获得机体在生命过程中基因的表达模式。RNA-seq是2008年建立起来的基于深度测序的转录组分析新技术，它能够在单核苷酸水平上对任何物种进行整体转录活动的检测［8］。与依赖基因组背景信息的基因芯片技术相比，该方法更适合基因组图谱尚未完成或遗传背景信息较为匮乏的物种［9，10］，是目前在全基因组水平上研究基因表达模式的主导技术。本研究拟采用RNA-seq技术对所构建的楮头红叶片的转录组进行测定，并对此全基因组水平转录组进行全局分析，为重要性状相关基因的克隆及功能分析、鉴定次级代谢物生物合成相关基因奠定基础，对于未来通过基因工程调控有用成分的含量以及相关药物的生产具有重要意义。

1 材料与方法

1.1 材料

试验材料为楮头红的新鲜叶片，叶片全展开，宽约1.5 cm，于2013年11月采自湖北省利川市毛坝新华村，采集后液氮速冻，并转移至-80℃保存备用。

1.2 方法

1.2.1 叶片总RNA的提取和检测 6片新鲜的叶片混合用于RNA的提取，采用Trizol（Invitrogen，USA）法提取叶片总RNA，提取过程中所用器具和耗材都经过RNase free处理，提取步骤参照说明书进行。琼脂糖凝胶电泳检测总RNA完整性，微量紫外分光光度计（NaNoDrop1000）检测RNA的纯度和浓度。

1.2.2 楮头红叶片转录组测序 转录组测序由深圳华大基因研究院提供技术服务，测序平台为Illumina HiSeqTM2000。

1.2.3 测序数据的组装和注释 测序后得到的原始reads，去除含有接头、重复的和测序质量很低的reads后，得到clean reads。使用短reads组装软件Trinity对clean reads进行从头组装［11］。过滤和组装以后得到的高质量的Unigenes，对这些从头组装的Unigenes进行注释。用NCBI蛋白质数据库（Nr）、非冗余核苷酸数据库（Nt）、基因本体论（GO）、直系同源基因簇（COG）和京都基因与基因组百科全书（KEGG）数据库作参考，得到Unigene的功能注释信息。根据Nr注释信息，使用Blast2GO［12］和WEGO［13］软件对Unigene进行GO注释和功能分类。

2 结果

2.1 转录组测序和组装

通过Illumina Hiseq2000平台测序，总计产出58 187 602条reads，去除低质量的和含有接头的reads以后，得到54 795 874条clean reads，共计4 931 628 660个核苷酸（Nucleotides，nt）。G+C含量、Q20和未知碱基序列分别为52.59%、98.33%和0.00%（表1）。

Raw reads总数量　Clean reads总数量　Clean reads碱基数/nt　Q20/%　N值/%　G+C含量/% 58 187 602　54 795 874　4 931 628 660　98.33　　0.00　52.59

利用Trinity软件对这些reads进行组装得到平均长度为347 nt（N50=648 nt）的Contigs 120 934条，共计41 909 136 nt。对Contigs进行去冗余以后最终得到51 305条Unigenes，平均长度为921 nt， N50为1 490 nt。Unigenes的长度分布显示（图1），长度大于1 000 nt的Unigenes有16 999条，占全部Unigenes的33.13%。说明本研究中转录组文库的测序和组装结果都较好，能够进行后续生物信息学分析。

图1 Unigenes的长度分布图

2.2 Unigenes的功能注释

通过blastx将Unigene序列比对到NCBI上的蛋白数据库Nr、Swiss-Prot以及KEGG（http：//www. genome.jp/kegg/pathway.html） 和COG（http：//www. ncbi.nlm.nih.gov/COG）（E-value＜1e-5），并通过blastn将Unigene比对到核酸数据Nt（E-value＜1e-5），得到与给定Unigene具有最高序列相似性的蛋白，从而得到该Unigene的蛋白功能注释信息。

注释结果显示共有40 532（79.0%）的Unigenes是有注释的（表2），其中匹配到Nr数据库中的有39 781条，占全部Unigenes的77.53%。此外，有21%的Unigenes无法注释，这可能与楮头红的遗传背景信息较少有关。注释到Nr数据库中Unigenes的E-value分布显示（图2-A），比对到的物种序列E-value均小于1e-5，其中E-value小于1e-100的有32.8%，说明比对结果的可信度较高。对Unigene的注释结果进行物种相似性分析显示（图2-B），序列比对相似性在15%-100%之间，其中大部分序列相似性在40%-80%之间，序列相似度大于80%的有21.4%，说明楮头红转录组的功能注释结果较好。

表2 楮头红叶片转录组中Unigenes的注释统计

图2 Unigenes在Nr库中的E-value分布（A）、相似性分布（B）及物种分布（C）

注释基因的同源序列的物种分布情况见图2-C，注释到葡萄（Vitis vinifera）的序列有21.6%；其次是蓖麻（Ricinus communis），有15.0%。这是因为葡萄和蓖麻具有丰富的基因组信息，为本研究中转录组的注释提供了参考序列［14，15］。

Gene Ontology（简称GO）是一个国际标准化的基因功能分类体系，提供了一套动态更新的标准词汇表（controlled vocabulary）来全面描述生物体中基因和基因产物的属性。据NCBI 数据库注释信息，使用Blast2GO软件得到Unigene的GO注释信息，然后用WEGO［13］对所有Unigene做GO功能分类统计，从宏观上认识楮头红的基因功能分布特征。对楮头红转录组Unigenes进行GO分析发现，有31 153条Unigenes注释到GO数据库，注释比例为60.72%。注释到生物学过程的基因最多，为129 573个，其次是细胞组成（96 967个），分子功能的最少，只有34 741个。GO分析的这3个ontology又分为56个亚类，如在生物过程中，细胞过程和代谢过程所占比例较高，细胞和细胞器部分在细胞组成所占比例较高，连接和催化活性在分子功能中占有较高比例（图3）。

COG注释表明，本研究中得到的注释到COG中的15 234条Unigenes分布于25个基因家族（图4），如RNA加工与修饰、染色体结构和动力学、能量产生与运输、细胞周期控制、细胞分裂及染色体分裂等。在25类基因家族中，注释最多的是一般功能预测（R），其次是转录（K）。值得注意的是，有4 564条Unigenes被注释到次生代谢物的生物合成、运输及分解代谢，为后续研究楮头红中有用化学成分相关的基因奠定了良好的基础。

2.3 Unigenes的代谢通路分析

对楮头红叶片的转录组进行KEGG分析发现，有23 898条Unigenes注释到KEGG数据库中，分布于128条已知的通路中，包括类黄酮生物合成（169条，ko00941）、类胡萝卜素生物合成（160条，ko00906）、环烯醚萜生物合成（144条，ko00900）。注释序列数目较多的3个通路分别是代谢途径、次生代谢物生物合成和植物激素信号转导（表3）。此外，从KEGG分析中，鉴定出了78个编码细胞色素P450的Unigenes。这些注释为后续次生代谢物的合成和代谢提供了很多有价值的信息。

2.4 预测编码蛋白框

按照Nr、Swiss-Prot、KEGG和COG的优先级顺序对楮头红叶片转录组中51 305条Unigenes进行blastx比对（E-value＜1e-5），以确定该Unigene的编码区序列。然后将编码区序列翻译成氨基酸序列，得到该Unigene编码区的氨基酸序列和核酸序列（序列方向5'→3'）。最后，用ESTScan软件预测无法比对到以上蛋白质库的Unigenes编码区，得到其编码区的氨基酸序列和核酸序列（序列方向5'→3'）。将所有的Unigenes在蛋白质数据库中进行blast比对以后，发现有39 302个Coding Sequence（CDS），其中有74.32%（29 211个）的序列长度大于1 000 nt（图5）。在ESTscan预测了2 065个CDS（图6），序列长度范围分布在200-2 000 nt之间，其中片段长度为500-1 000 nt的有233条，比例约为11.28%。说明Unigenes的序列质量较好，共预测了41 367个CDS。

3 讨论

近年来，新一代高通量测序技术的建立和发展，为非模式植物功能基因的研究提供了有效手段［16］，也为植物基因组的研究提供了便利。然而，与模式植物和重要的农作物相比，目前对药用植物基因组的研究相对较少［17］，RNA-seq技术为植物转录组研究提供了新的发展契机［18］。通过药用植物的研究，可以清楚地了解其生物学和生物医学功能的基本信息［19，20］。本研究首次利用RNA-seq高通量测序技术对楮头红叶片的转录组进行测序和转录组结果分析，为其次生代谢物合成相关基因的挖掘、研究基因功能提供了基础信息，同时为其有效成分的生物合成途径和调控机制的研究奠定了基础。

本研究中共得到51 305条质量较高的Unigenes，平均长度为921 nt，注释比例高达79.0%，其Unigenes数目、平均长度和注释比例都远高于基于Roche 454的丹参（Salvia miltiorrhiza）转录组（46 722条Unigenes，平均长度为414 nt，注释比例约为28.48%）［21］，美国西洋参（American ginseng）的转录组（31 088条序列，平均长度为427 nt，注释比例为69.8%）［22］。在丹参和西洋参的研究中，测序平台为Roche 454，而在本研究中，测序平台为Illumina。Roche 454平均读长较长，但是准确率低，Illumina性价比较高，覆盖度较深，经过技术升级后，平均读长也得到了改善，是目前使用较为广泛的测序平台之一［23］。Unigenes的片段长度对注释效率也有显著的影响，本研究中Unigenes注释到Nr、Nt、Swiss-Prot、KEGG、COG和GO库中的比例也相对较高，分别为77.53%、56.18%、53.14%、46.58%、29.69%和60.72%。在蛋白质数据库中对所有的Unigenes进行blast发现，有39 302个CDS，用ESTscan预测了2 065个CDS。通过与以上两种药用植物的测序结果相比，进一步说明了本研究转录组测序的产量更高、质量更好。

图3 Unigene在GO库中的分布

图4 Unigene在COG库中的分布

编号　　途径　　注释基因数/23 898　　占总基因数的比例/%　　代谢通路ID 1新陈代谢通路　5810 　24.31　ko01100 2次生代谢物生物合成　2323 　9.72　ko01110 3植物激素信号转导　1532　6.41　ko04075 4内吞作用　1461 　6.11　ko04144 5植物-病原菌感染　1453 　6.08　ko04626 6甘油磷脂代谢　1350　5.65　ko00564 7醚脂类代谢　1193　4.99　ko00565 8 RNA运输　1084　4.54　ko03013 9剪接体844 　3.53　ko03040 10　　淀粉和蔗糖代谢　800 　3.35　ko00500 11　　内质网蛋白加工　740　3.1　ko04141 12　mRNA监测途径　681 　2.85　ko03015 13　　核糖体　650 　2.72　ko03010 14　　泛素介导的蛋白水解作用　603　2.52　ko04120 15　　戊糖和葡萄糖醛酸酯互变　511　2.14　ko00040 16　　嘌呤代谢　501 　2.1　ko00230 17　　嘧啶代谢　434　1.82　ko00240 18　　氧化磷酸化　383　1.6　ko00190 19　RNA降解　354　1.48　ko03018 20　　真核生物核糖体合成　334　1.4　ko03008

图5 Blast预测的Unigenes的CDS长度分布图

图6 ESTscan预测的Unigenes的CDS长度分布图

细胞色素P450（Cytochrome P450，CYP450）是广泛存在于几乎所有生物体内的代谢酶，参与了内源性生理物质如类固醇、细胞激素和脂肪酸等物质的代谢［24］，是植物中所占比例最大的一类蛋白家族，负责大多数次生代谢物的氧化过程［25，26］。本研究中KEGG通路分析显示，参与次生代谢物生物合成的Unigenes有2 323条，占全部Unigenes的9.72%，其中有78条Unigenes编码了细胞色素P450家族蛋白。除此之外，还有参与二萜类生物合成、异黄酮生物合成、生物碱合成等次生代谢物合成的Unigenes。

研究表明，楮头红中主要的活性物质为黄酮类化合物，黄酮是植物体内具有重要作用的次生代谢物，具有抗氧化、清除自由基、抗肿瘤等广泛的药理活性［27］。查尔酮异构酶（Chalconeisomerase，CHI）是黄酮类化合物代谢途径中的一个关键酶，能够催化分子内环化反应，从而使双环的查尔酮转化为有生物活性的三环（2S）-黄烷酮［28］。目前已在灯盏菊［29］、水母雪兔子［30］、金荞麦［31］等药用植物中分离出了CHI基因。本研究中有4条Unigenes编码了CHI，这对于调控楮头红中黄酮类化合物的合成提供了可靠的依据。

植物生物碱是一种广泛存在于高等和低等植物中的含氮有机化合物，是许多中草药和药用植物的有效成分［32］。据报道，大约有20%的有花植物都能产生生物碱［33］。由于目前对楮头红的有效成分研究较少，尚未发现楮头红中有生物碱成分的报道。然而，基于转录组测序的优势，我们发现在KEGG通 路 中 有7条（KO00901、KO00950、KO00960、KO00900、KO00904、KO00909、KO00902） 是 与生物碱合成相关的代谢途径，包括萜类合成途径、吲哚合成途径、异喹啉合成途径等，共有362条Unigenes，这从侧面反映了楮头红中含有不同类型的生物碱，为其药用成分的开发和应用领域的拓展奠定了基础。

4 结论

利用Illumina测序技术对楮头红叶片的转录组进行测序和分析，揭示了其转录组的整体表达模式，初步获得了一些参与次生代谢物合成的基因序列信息，共获得51 305条质量较高的Unigenes，其中有注释的40 532条，注释比例高达79%。初步获得了一些与次生代谢物合成相关的基因序列，为深入开展药用植物生物活性成分的合成和鉴定以及功能基因的筛选提供了丰富了数据资料。在利用生物技术提高药用植物有效成分含量以及相关药物研发和生产方面，具有重要的应用价值。

［1］ Wang XM， Wan CP， Zhou SR， Qiu Y. Two new flavonol glycosides from Sarcopyramis bodinieri var. delicate［J］. Molecules， 2008， 13（6）：1399-1405.

［2］ Wan C， Zheng X， Chen H， et al. Flavonoid constituents from herbs of Sarcopyramis bodinieri var. delicata［J］. China Journal of Chinese Materia Medica， 2009， 34（2）：172-174.

［3］ 陈桂菲， 叶淑玲. 风柜斗草药用探讨［J］. 时珍国医国药，2000， 11（2）：129.

［4］ Wei WB， Huang YJ， Cong HJ， et al. Chemical constituents from Sarcopyramis nepalensis Wall［J］. Phytochemistry Letters， 2014， 8：101-104.

［5］ Zhang JW， Liao M， Chen HD， Zhang YH. Structural elucidation of a new flavone from Sarcopyramis nepalensis［J］. Journal of Asian Natural Products Research， 2011， 13（3）：256-259.

［6］ 陈豪， 何丽君， 林丽芳. 民间草药风柜斗草研究进展［J］. 海峡药学， 2012， 24（10）：55-56.

［7］ Wang Z， Gerstein M， Snyder M. RNA-Seq：a revolutionary tool for transcriptomics［J］. Nature Reviews Genetics， 2009， 10（1）：57-63.

［8］ Costa V， Angelini C， De Feis I， Ciccodicola A. Uncovering the complexity of transcriptomes with RNA-Seq［J］. Journal of Biomedicine and Biotechnology， 2010， 2010， doi：10.1155/2010/1853916.

［9］ Birzele F， Schaub J， Rust W， et al. Into the unknown：expression profiling without genome sequence information in CHO by next generation sequencing［J］. Nucleic Acids Research， 2010， 38（12）：3999-4010.

［10］ Wang ET， Sandberg R， Luo S， et al. Alternative isoform regulation in human tissue transcriptomes［J］. Nature， 2008， 456（7221）：470-476.

［11］ Grabherr MG， Haas BJ， Yassour M， et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome［J］. Nature Biotechnology， 2011， 29（7）：644-652.

［12］ Conesa A， Götz S， García-Gómez JM， et al. Blast2GO：a universal tool for annotation， visualization and analysis in functional genomics research［J］. Bioinformatics， 2005， 21（18）：3674-3676.

［13］ Ye J， Fang L， Zheng H， et al. WEGO：a web tool for plotting GO annotations［J］. Nucleic Acids Research， 2006， 34（suppl. 2）：W293-W297.

［14］ Chan AP， Crabtree J， Zhao Q， et al. Draft genome sequence of the oilseed species Ricinus communis［J］. Nature Biotechnology，2010， 28（9）：951-956.

［15］ Jaillon O， Aury JM， Noel B， et al. The grapevine genome sequence suggests ancestral hexaploidization in major angiosperm phyla［J］. Nature， 2007， 449（7161）：463-467.

［16］ 刘玉林， 李伟， 张志翔. 基于高通量测序的辽东栎转录组学研究［J］. 生物技术通报， 2014（7）：119-124.

［17］ Chen S， Xiang L， Guo X， Li Q. An introduction to the medicinal plant genome project［J］. Frontiers of Medicine， 2011， 5（2）：178-184.

［18］ 付畅， 黄宇. 转录组学平台技术及其在植物抗逆分子生物学中的应用［J］. 生物技术通报， 2011（6）：40-46.

［19］ Collins FS， Green ED， Guttmacher AE， Guyer MS. A vision for the future of genomics research. pdf［J］. Nature， 2003， 422（6934）：835-847.

［20］ Lockhart DJ， Winzeler EA. Genomics， gene expression and DNA arrays［J］. Nature， 2000， 405（6788）：827-836.

［21］ 李滢， 孙超， 罗红梅， 等. 基于高通量测序 454 GS FLX 的丹参转录组学研究［J］. 药学学报， 2010， 45（4）：524-529.

［22］ Sun C， Li Y， Wu Q， et al. De novo sequencing and analysis of the American ginseng root transcriptome using a GS FLX Titanium platform to discover putative genes involved in ginsenoside biosynthesis［J］. BMC Genomics， 2010， 11：262.

［23］ Zhang J， Chiodini R， Badr A， Zhang G. The impact of nextgeneration sequencing on genomics［J］. Journal of Genetics and Genomics， 2011， 38（3）：95-109.

［24］ Jørgensen A， Giessing A， Rasmussen LJ， Andersen O. Biotransformation of polycyclic aromatic hydrocarbons in marine polychaetes［J］. Marine Environmental Research， 2008， 65（2）：171-186.

［25］ Coon MJ. Cytochrome P450：nature’s most versatile biological catalyst［J］. Annual Review of Pharmacology and Toxicology，2005， 45：1-25.

［26］ Morant M， Bak S， Møller BL， Werck-Reichhart D. Plant cytochromes P450：tools for pharmacology， plant protection and phytoremediation［J］. Current Opinion in Biotechnology， 2003，14（2）：151-162.

［27］ 程秋月， 郭菁， 张成义. 黄酮类化合物药理作用的研究［J］.北华大学学报：自然科学版， 2012， 12（2）：180-183.

［28］ 李琳玲， 程华， 许锋， 程水源. 植物查尔酮异构酶研究进展［J］.生物技术通讯， 2008（6）：935-937.

［29］ Zhang YF， Zhang Y， Cui MK， Yan SQ. The Clone of cDNA enco-ding chalcone isomerases and its information analyze in E. breviscapus（vant）Hand. -Mazz［J］. Procedia Engineering， 2011， 18：194-199.

［30］ Li FX， Jin ZP， Zhao DX， et al. Overexpression of the Saussurea medusa chalcone isomerase gene in S. involucrata hairy root cultures enhances their biosynthesis of apigenin［J］. Phytochemistry， 2006， 67（6）：553-560.

［31］ 雒晓鹏， 白悦辰， 高飞， 等. 金荞麦查尔酮异构酶的基因克隆及其花期表达与黄酮量分析［J］. 中草药， 2013， 11：1481-1485.

［32］ Kuete V. 21 -Health effects of alkaloids from african medicinal plants//［M］. Toxicological Survey of African Medicinal Plants. 2014：611-633.

［33］ 邢世海， 王荃， 潘琪芳， 等. 长春花萜类吲哚生物碱的生物合成途径［J］. 西北植物学报， 2012， 32（9）：1917-1127.

（责任编辑 李楠）

Transcriptome Analysis of Sarcopyramis nepalensis via RNA-seq Technology

Jin Hong1Jiao Genlin1Chen Gang2
（1. Fairy Lake Botanical Garden，Shenzhen and CAS，Shenzhen518004；2. College of Life Science，Zhaoqing University，Zhaoqing526061）

Transcriptome analysis of Sarcopyramis nepalensis leaves was performed via a newly developed high-throughput sequencing technology（Illumina RNA-seq）. A total of 51 305 unigenes were generated with 921 nt of average length and 1 490 nt of unigene N50 after filtering and assembly of original reads. These unigenes from the de novo assembly were further annotated using BLAST and BLAST2GO softwares. A total of 40 532 unigenes annotated with databases of non-redundant protein sequence（Nr）， non-redundant nucleotide（Nt），Swiss-Prot， Gene Ontology database（GO）， Clusters of Orthologous Groups（COG）and Kyoto Encyclopedia of Genes and Genomes（KEGG）databases available at NCBI as

. The proportion of unigenes annotated in Nr， Nt， Swiss-Prot， KEGG， COG and GO databases were 77.53%， 56.18%， 53.14%， 46.58%， 29.69% and 60.72%， respectively. Total 39 302 CDSs were obtained using blast in protein databases，and 2 065 CDSs were predicted using ESTscan software. KEGG pathway parsing revealed that 2 323（9.72%）unigenes were involved in biosynthesis of secondary metabolites（KO01110）， and 78 unigenes encoding the cytochrome P450 family proteins were identified. These annotated information provided theoretical foundation for determining the vital genes involved in biosynthesis of secondary metabolites of medicinal plants.

Sarcopyramis nepalensis transcriptome；RNA-seq

10.13560/j.cnki.biotech.bull.1985.2015.05.014

2014-09-11

深圳市城市管理局项目（201318）

金红，博士，高级工程师，研究方向：植物多样性保护及利用；E-mail：jinhong@szum.gov.cn