基于高通量测序的药用植物“凤丹”根皮的转录组分析

谢冬梅 俞年军 黄璐琦 彭代银 +刘丛彬 朱月健 黄浩

[摘要]牡丹皮為我国传统常用中药，安徽省铜陵地区栽培的“凤丹”根皮加工而成的药材牡丹皮被誉为道地药材，药用活性成分丰富多样，但目前尚不清楚 “凤丹”药用部位次生代谢过程中活性物质合成的遗传学基础。研究采用Illumina HiSeq 4000高通量测序平台对五年生“凤丹”根皮转录组进行测序，对测序结果进行de novo 拼接和功能注释，测序后获得72 997条unigene。进一步利用公共数据库进行同源比对，其中41 139条unigene被Nr数据库成功注释，34 952条unigene能被GO数据库成功注释，20 016条unigene被KEGG数据库成功舒注释，共涉及到5个大类、34个种类、352条代谢通路；在次生物质合成与代谢途径中，其中苯丙素类化合物、萜类化合物骨架合成、各种类型萜类化合物、生物碱类化合物以及黄酮类成分生物合成途径中的unigene分别有214，104，152，55，36个；不同产地样本间差异表达基因的富集性比较显示不同产地样本间存在明显差异；此外，在72 997条unigene中共检测到9 939个SSR序列，其中二核苷酸重复的SSR标记占2075%。研究的结果不仅为挖掘“凤丹”次生代谢物生物合成关键基因提供了基础数据信息，也为药用牡丹的遗传多样性研究和分子标记开发奠定了分子基础。

[关键词]牡丹皮； 转录组； 次生代谢； 差异基因； 简单重复序列

Next generation sequencing and transcriptome analysis of

root bark from Paeonia suffruticosa cv Feng Dan

XIE Dongmei1，2， YU Nianjun1*， HUANG Luqi2*， PENG Daiyin1， LIU Congbin1， ZHU Yuejian3， HUANG Hao4

（1 Institute of Traditional Chinese Medicine Resources Protection and Development， Anhui Academy of Chinese

Medicine， Anhui University of Chinese Medicine， Hefei 230012， China；

2 State Key Laboratory of Daodi Herbs， National Resource Center for Chinese Materia Medica， China

Academy of Chinese Medical Sciences， Beijing 100700， China；

3 Anhui Jiren Pharmaceutical Co， Ltd， Bozhou 236800， China；

4 Beijing Tongrentang Anhui Traditional Chinese Medicinal Materials Co， Ltd， Tongling 244000， China）

[Abstract]Moutan Cortex is an important traditional Chinese medicine， “Fengdan Pi” was known as Daodi herbs from the root bark of Paeonia suffruticosa cv Feng Dan for its extracted various active components However， the genetic basis for their activity is virtually unknown The transcriptome of the root bark from “Fengdan” was sequenced using the Illumina HiSeq 4000 sequencing platform The clean reads were then de novo assembled into 72 997 unigenes Among them， the number of unigenes which could been annotated by dataset Nr and GO was 41 139 and 34 592. The 20 016 unigenes could been annotated by KEGG dataset， which were involved in 5 major categories， 34 middle categories， and 352 metabolism pathways. The number of unigenes which were mapped to the phenylpropanoid biosynthesis pathway， terpenoid backbone biosynthesis pathway， terpenoid biosynthesis pathway， alkaloid biosynthesis pathway， and flavonoid biosynthesis pathway was 214， 104， 152， 55 and 36 respectively， suggesting that they are involves in these pathways of pharmaceutically important Furthermore， there also showed remarkable differences in groups which enrichment ratio of the different expressed gene compared. In addition， a total of 9 939 SSRs were identified from the sequence of 72 997 unigenes This study not only provides many valuable basal data which was important gene in the synthesis pathway of secondary metabolites with gene searching， but also has important significance to find molecular marker in germplasm for breeding and improvement

[Key words]Moutan Cortex； transcriptome； secondary metabolism； different expressed gene； simple sequence repeat

中药牡丹皮Moutan Cortex来源于毛茛科植物牡丹Paeonia suffruticosa Andr的干燥根皮，具有清热凉血，活血化瘀的功效[1]。现代中药化学和药理研究表明，牡丹皮主要化学成分为酚类、酚苷类、单萜及单萜苷类以及三萜、甾醇及其苷类、黄酮、有机酸、香豆素等，具有抗炎抗菌、降血糖以及心血管系统保护、中枢神经保护、增强免疫、止血、凝血等药理作用[23]。安徽为药材牡丹皮的主产区，其中铜陵地区“凤丹”P suffruticosa Andr cv Feng Dan的干燥根皮加工而成的牡丹皮历来被为誉为道地药材[45]。

随着中药现代化的发展，关于牡丹皮的道地性研究也逐渐深入。采用药物分析技术，研究人员不仅对其最佳采收期、不同产地来源药材中主要活性成分丹皮酚和芍药苷的含量进行了研究[67]，也对不同产地“凤丹”栽培后所获得药材的特征图谱、QTOF图谱、挥发油成分、无机元素等进行了比较研究，探讨药材道地性的形成机制[8]。此外，利用生物技术研究土壤微生物的整体性和多样性也逐步用于揭示中药材牡丹皮的道地性[9]。

随着中药现代研究的逐步推进，高通量测序（next generation sequencing，NGS）技术的不断进步与实验流程的不断完善，药用植物转录组研究成为一个新的研究热点。自采用新一代高通量测序技术对二年生丹参根的转录组进行测序研究以来[10]，陆续有30余种药用植物开展了转录组研究[11]。利用转录组测序技术（RNA sequencing，RNASeq）可以较好地应用于捕捉道地药材与非道地药材不同样品间差异表达的基因，从整体上对参与生长发育和次生代谢基因进行系统的研究，进而更好地阐释药材道地性形成的分子机制[12]。一些研究者虽然对牡丹低温刺激下转录组的变化进行了研究[13]，而对于道地产区与非道地产区“凤丹”根皮中差异表达的基因缺少相关研究。本研究利用第二代测序技术对不同产地栽培的 “鳳丹”根皮的转录组进行了从头组装（de novo assembly），丰富药用牡丹的基因数据库；同时，利用比较转录组学，分析不同产地、五年生、供采收的“凤丹”根皮部基因差异表达的富集情况，为药用牡丹根皮中活性成分的生物合成研究以及药材道地性品质的形成提供生物学依据；通过SSR位点的信息分析，为开展分子标记辅助育种、基因工程技术选育创制新的药用牡丹优良品种奠定基础。

1材料和方法

11样品于药材最佳采收期，分别从位于药材牡丹皮的传统道地产区安徽省铜陵市凤凰山（TLDP）及其周边地区芜湖市南陵县（WHDP） 的北京同仁堂安徽中药材有限公司GAP栽培基地、非传统道地产区亳州市（BZDP）的安徽济人药业有限公司牡丹皮GAP栽培基地采挖五年生药用植物“凤丹”的主根，迅速剥取其根皮置液氮中冷冻，后转入-80 ℃保存备用。

12总RNA的提取、文库构建和转录组测序使用TRIzol@Reagent试剂盒、Plant RNA Purification Reagent试剂盒（Invitrogen，美国） 提取、纯化根皮总RNA，琼脂糖凝胶电泳检测RNA完整性，Nanodrop 2000超微量分光光度计（Thmermo，美国）对所提总RNA的浓度和纯度进行检测，Agilent 2100生物分析仪（Agilent Technology，美国）测定RIN值，以RNA条带清晰、28/23S亮度大于18/16S，60≤RIN<80为质量合格。为了获得尽可能丰富的“凤丹”根皮采收期转录组信息，选取不同栽培基地各6株“凤丹”根皮RNA等量混合为1份，得到3份“凤丹”根皮混合总RNA样品后送美吉生物完成反转合成cDNA、连接adaptor以及HiSeq 4000（Illumina，美国）上机测序等后续工作。

13测序数据过滤及转录组的组装使用SeqPrep （https：//githubcom/jstjohn/SeqPrep）、Sickle （https：//githubcom/najoshi/sickle） 软件将测序所得到的原始数据（raw reads）进行过滤以去除rRNA、含接头的低质量reads以及含N比率超过10%的reads，使得Q30（碱基的测序错误率小于01%）达到80%以上，获得去掉接头后的测序序列（clean reads）用于后续分析。

由于牡丹基因组测序尚未完成，无参考序列，因此使用软件Trinity[14]进行“凤丹”根皮部的转录组从头组装（de novo assembly）：通过测序序列之间的重叠（overlap）信息组装得到重叠群（contigs），依据序列的双端信息（pairedend）信息和重叠群之间的相似性对其进行聚类整合和延长，进而组装得到转录本（transcripts），并从中选取最长的转录本作为非重复序列基因（unigenes）。Trinity参数设置为Kmer=25，序列延伸成contig时overlap为Kmer1。

差异表达基因的分析与筛选是建立在同一套参考基因的基础上，因此对3个样品得到的unigene数据通过cdhit聚类去除冗余，利用TGIGL的聚类组装策略最终得到非冗余的“凤丹”采收期根皮的Allunigenes数据库。

14功能基因注释对“凤丹”根皮Allunigenes转录数据使用BlastX软件（Evalue<1×10-5）对unigene序列与NR库（nonredundant database，http：//wwwncbinlmnihgov/），COG数据库（clusters of orthologous group，http：//www. ncbi. nlm. nih. gov/COG/）和SwissProt数据库（http：//wwwuniprotorg/）比对，进行功能注释和分类处理；再对unigene序列进行GO（http：//www. geneontology. org/）功能注释和分类，并用软件对GO注释结果分类作图，并将unigene与KEGG数据（kyoto encyclopedia of genes and genomes，http：//www. genome. jp/kegg/）进行比对，分析相关的代谢通路。

15基因的表达水平及表达差异分析Illumina HiSeq 4000测序平台得到的reads一般较短、插入删除错误较少，因此选择短序列比对软件Bowtie （http：//bowtiebio. sourceforge. net/index. shtml）完成两两样本间的转录组数据比对；利用bowtie的比对结果，使用软件RSEM（http：//deweylab. biostat. wisc. edu/rsem/）计算3个样品比对到每个基因上的read count数目，然后对其极性FRKM转换，进而获得不同样本基因的表达水平[15]；根据RSEM得到的gene read count，使用软件edgeR（http：//www. bioconductor. org/packages/2. 12/bioc/html/edgeR. html）进行样品间的差异基因分析。

为了进一步研究道地与非道地产区“凤丹”根皮中差异表达基因的富集水平，使用软件Goatools（http：//github. com/tanghaibao/GOatools）对样品间差异表达基因进行GO库注释的分类统计和功能富集分析，并使用Fisher精确检验法对差异表达基因富集的显著性进行评价；使用KOBAS（http：//kobas. cbi. pku. edu. cn/home. do）进行KEGG pathway富集分析，使用Fisher精确检验法对差异表达基因KEGG通路富集的显著性进行评价。

16SSR分析简单重复序列（simple sequence repeats，SSR）在真核生物基因组中的重复次数表现为个体间高度变异，数量丰富，具有广泛的应用性。因此，采用MISA（http：//pgrc. ipkgatersleben. de/misa/misa. html）对unigene进行SSR检测，检测中对应的单核苷酸（mononucleotide）、二核苷酸（dinucleotide）、三核苷酸（trinucleotide）、四核苷酸（tetranucleotide）、五核苷酸（pentanucleotide）以及六核苷酸（hexanucleotide）等SSR重復类型的最少重复次数分别设置为10，6，5，5，5，5。

2结果与分析

21“凤丹”根皮RNASeq测序及质量评估利用Illumina HiSeq 4000高通量测序技术对3个产地来源的“凤丹”根皮转录组进行测序，得到raw reads，clean reads的数量及测序长度等见表1。3个样本的转录组clean reads中Q20，Q30均大于90%，GC量均在45%左右，表明测序质量较好，数据可用于后续AllUnigenes数据库的建立和基因功能注释、分类及样本间异同性分析。

22“凤丹”根皮转录本组装与分析利用Trinity软件对上述3份获得的clean reads去除重复部分，并根据序列之间的重叠部分进行混合组装，并通过相似性比较后组装获得AllUnigenes数据库的转录本72 997条，转录本总长度达到71 528 823 bp，最短的仅为201 bp，最长的达到18 123 bp，平均长度78207 bp，长度在200～400 bp的有40 123条，占4387%；获得unigene72 997条，总长度达51 290 894 bp，平均长度70264 bp。其中，长度在200～400 bp的有37 320条，占5113%；转录本和unigene的N50分别达到1 135，1 227 bp，均大于800 bp，表明组装片段的完整性较高，见图1。

23“凤丹”根皮转录组AllUnigenes基因功能注释及分类经过 BlastX 序列比对分析，AllUnigenes在NR，COG，GO，KEGG数据库中获得注释unigene的数量分别为41 139，8 627，24 952，20 016。其中，NR 数据库中搜索unigene 相似序列最多，为41 139条（564%）；通过COG注释，8 672条unigene获得了分类信息，这些基因分属于25个功能分类，见图2。1 063条unigene属一般功能预测项（1226%），是获得注释数量最多的部分；1 057条unigene（1219%）与翻译、核糖体结构与生源相关。此外，还有394条（247%）功能未知的unigene基因可能是“凤丹”特有的基因，与其供药用的所具有的发挥药效的次生代谢物质有关。

GO注释：在总共组装的72 997条unigene中，有24 952条unigene能被GO 数据库成功注释，分为生物学过程（biological process）、细胞成分（cellular component）及分子功能（molecular function）三大类

ARNA 加工与修饰；B染色质结构与动力学；C能量生产和转换；D细胞周期控制、细胞分裂、染色体分区；E氨基酸运输和代谢；F核苷酸运输和代谢；G碳水化合物的运输和代谢；H辅酶运输和代谢；I脂质运输和代谢；J翻译；核糖体结构和生物转化；K转录；L复制、重组和修复；M细胞被膜生物起源；N细胞运动性；O翻译后修饰，蛋白质转换与分子伴侣；P无机离子转运和代谢；Q次生代谢产物合成、运输与分解；R一般功能预测；S未知功能；T信号传导机制；U胞内运输、分泌及膜泡转运；V防御机制；W细胞外结构；Y核结构；Z细胞骨架。

共60个功能组。进一步研究发现： 24 952条unigene在GO库中归属于生物学过程中的24个功能组。其中，涉及代谢过程的unigene17 113条、细胞过程14 817条、单生物过程11 975条、定位3 413条、定位建立3 350条、生物调节3 297条、应激反应3 202条、生物学过程的调控3 098条，以参与代谢过程和细胞过程获得序列注释最多；同时，另有31 331条 unigene在GO库中归属于分子功能的16个功能组，其中获得注释较多的unigene分别为催化活性13 589条、结合功能13 207条，这些基因可能与牡丹皮中化学成分的催化与合成有关。此外，44 873条unigenen在GO库中归属于细胞成分的20个功能组。其中，被注释到较多的unigene分别为细胞9 669条、细胞器6 737条、细胞膜5 500条和复杂大分子4 363条。

KEGG注释：共20 016条unigene能够在KEGG数据库成功注释，涉及代谢（metabolism）、遗传信息处理（genetic information processing）、环境信息处理（environmental information processing）、细胞进程（cellular processes）及生物系统（organismal systems）5个大类及34个中类、352条代谢通路，其中代谢大类中有6 469条unigene被成功注释，其次是次生代谢产物的合成（biosynthesis of secondary metabolites）2 481条和不同环境的微生物代谢（microbial metabolism in diverse environment）1 735条。进一步研究还发现，在代谢途径中，与中药效活性成分生物合成相关的unigene序列分别是苯丙素类合成（phenylpropanoid biosynthesis）214条、萜类化合物骨架合成（terpenoid backbone biosynthesis）104条、各种萜类化合物合成（sesquiterpenoid，triterpenoid，monoterpenoid and diterpenoid biosynthesis）152条、（莨菪类、哌啶类、吡啶类）生物碱合成（tropane，piperidine and pyridine alkaloid biosynthesis）55条、黄酮类成分合成（flavonoid，flavone，flavonol and isoflavonoid biosynthesis）36条。

24差异表达基因的分析3个栽培产地样品间“凤丹”根皮总AllUnigenes比对的结果分别为BZDP（8318%），TLDP（8178%），WHDP（8621%），不同样品的转录本/基因的表达量概率（FPKM scores）密度分布见图3，其中以TLDP基因的密度最高；顯著差异表达基因的筛选标准为FDR<005且log2|FC|≥1，3个样品在72 997条根皮AllUnigenes转录本数据库能被识别到的差异基因总数为3 430条，不同样品间分组BZDP vs. TLDP，WHDP vs. TLDP，BZDP vs. WHDP的表达差异基因数分别为 1 502，878，1 050，各组间差异基因的数量分布状况见图4。可以看出，非道地产区样品与传统道地产区样品（BZDP vs. TLDP）间显著差异表达基因数目最多，传统道地产区铜陵及其周边南陵地区样品（WHDP vs. TLDP）差异表达基因数目最少；具有显著性富集（P<0001）的差异基因功能类型和GO库的ID号见表2。

由表2可以看出，BZDPTLDP间差异基因显著富集于生物学过程和分子功能的9个功能组，BZDPWHDP间差异基因显著富集于生物学过程、分子功能和细胞组成的22个功能组，TLDPWHDP间差异基因显著富集于生物学过程、分子功能和细胞组成的20个功能组。不同样品间差异基因在KEGG通路富集具有特别显著性差异（P<0001）的通路描述和ID号见表3。可以看出BZDPTLDP间差异基因显著富集主要是遗传信息和代谢的8个代谢通路，BZDPWHDP间差异基因显著富集于遗传信息、生物系统、代谢和疾病的8个代谢通路，TLDPWHDP间差异基因显著富集于代谢和疾病的6个代谢通路。

25SSR分子标记在72 997条AllUnigene库中，共检测到9 939条unigene中含有SSR，发生频率（含有SSR的unigene数量与总unigene数量之比）为136%。其中6 060条unigene序列只含有1个SSR位点，518条为复杂重复类型SSR，含有SSR标记的unigene序列数量及类型统计见表4。

统计结果显示，除1个SSR位点外，“凤丹”根皮AllUnigenes转录组中SSR的主要类型是二核苷酸，占SSR总数的2075%，出现频率最高的3个重复类型是CT/GA，AG/TC，AT/TA；其次是三核苷酸，占SSR总数的1223%；四、五、六核苷酸重复类型的数量很少，占SSR总数的085%；SSR位点的序列分布从10～219 bp不等，平均长度17 bp。

3讨论与结论

中药丹参[10]、人参[16]和虎杖[17]等的转录组测序多数是利用药用部位根进行，分别获得了18 235，31 741，86 418条unigenes，为最大限度的获得“凤

丹”根皮的转录组信息，研究选择安徽省传统道地产区和非道地产区、五年生“凤丹”根皮为试验材料，进行混合样品测序，共组装得到72 997条“凤丹”根皮的基因序列，加快了药用植物“凤丹”的分子生物学研究。

经过 BlastX 比对分析，436%的unigene在NR库中不能匹配到已知基因，可能由于NR 数据库中由于缺少牡丹基因组信息；COG分类与注释中，发现有318条unigene（367%）参与了次生代谢物质的生物合成、转运和降解，深入研究这些基因的功能将有助于揭示“凤丹”根皮中有效成分合成的调控机制；KEGG注释中，与次生代谢产物的合成途径相关的unigene有2 481条，与不同环境的微生物代谢相关的有1 735条，说明五年生“凤丹”的植物根皮中参与次生代谢产物合成与代谢的基因丰富，为进一步研究“凤丹”道地药材品质的形成提供了依据。来自于非道地产区样品与传统道地产区“凤丹”根皮中差异表达基因数目的多少可以较好地反映在

药用种质遗传背景一致的情况下，道地产区的环境、气候对于药材道地性品质的形成具有较大的影响，传统道地产区铜陵及其相邻地区南陵因生态环境、气候相近，两地来源样本间在次生代谢物合成中差异表达基因的数量远远少于环境、气候距离较远的亳州栽培“凤丹”根皮样本。