基于高通量測(cè)序的藥用植物“鳳丹”根皮的轉(zhuǎn)錄組分析

謝冬梅 俞年軍 黃璐琦 彭代銀 +劉叢彬 朱月健 黃浩

[摘要]牡丹皮為我國(guó)傳統(tǒng)常用中藥，安徽省銅陵地區(qū)栽培的“鳳丹”根皮加工而成的藥材牡丹皮被譽(yù)為道地藥材，藥用活性成分豐富多樣，但目前尚不清楚 “鳳丹”藥用部位次生代謝過(guò)程中活性物質(zhì)合成的遺傳學(xué)基礎(chǔ)。研究采用Illumina HiSeq 4000高通量測(cè)序平臺(tái)對(duì)五年生“鳳丹”根皮轉(zhuǎn)錄組進(jìn)行測(cè)序，對(duì)測(cè)序結(jié)果進(jìn)行de novo 拼接和功能注釋?zhuān)瑴y(cè)序后獲得72 997條unigene。進(jìn)一步利用公共數(shù)據(jù)庫(kù)進(jìn)行同源比對(duì)，其中41 139條unigene被Nr數(shù)據(jù)庫(kù)成功注釋?zhuān)?4 952條unigene能被GO數(shù)據(jù)庫(kù)成功注釋?zhuān)?0 016條unigene被KEGG數(shù)據(jù)庫(kù)成功舒注釋?zhuān)采婕暗?個(gè)大類(lèi)、34個(gè)種類(lèi)、352條代謝通路；在次生物質(zhì)合成與代謝途徑中，其中苯丙素類(lèi)化合物、萜類(lèi)化合物骨架合成、各種類(lèi)型萜類(lèi)化合物、生物堿類(lèi)化合物以及黃酮類(lèi)成分生物合成途徑中的unigene分別有214，104，152，55，36個(gè)；不同產(chǎn)地樣本間差異表達(dá)基因的富集性比較顯示不同產(chǎn)地樣本間存在明顯差異；此外，在72 997條unigene中共檢測(cè)到9 939個(gè)SSR序列，其中二核苷酸重復(fù)的SSR標(biāo)記占2075%。研究的結(jié)果不僅為挖掘“鳳丹”次生代謝物生物合成關(guān)鍵基因提供了基礎(chǔ)數(shù)據(jù)信息，也為藥用牡丹的遺傳多樣性研究和分子標(biāo)記開(kāi)發(fā)奠定了分子基礎(chǔ)。

[關(guān)鍵詞]牡丹皮； 轉(zhuǎn)錄組； 次生代謝； 差異基因； 簡(jiǎn)單重復(fù)序列

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來(lái)源于毛茛科植物牡丹Paeonia suffruticosa Andr的干燥根皮，具有清熱涼血，活血化瘀的功效[1]。現(xiàn)代中藥化學(xué)和藥理研究表明，牡丹皮主要化學(xué)成分為酚類(lèi)、酚苷類(lèi)、單萜及單萜苷類(lèi)以及三萜、甾醇及其苷類(lèi)、黃酮、有機(jī)酸、香豆素等，具有抗炎抗菌、降血糖以及心血管系統(tǒng)保護(hù)、中樞神經(jīng)保護(hù)、增強(qiáng)免疫、止血、凝血等藥理作用[23]。安徽為藥材牡丹皮的主產(chǎn)區(qū)，其中銅陵地區(qū)“鳳丹”P(pán) suffruticosa Andr cv Feng Dan的干燥根皮加工而成的牡丹皮歷來(lái)被為譽(yù)為道地藥材[45]。

隨著中藥現(xiàn)代化的發(fā)展，關(guān)于牡丹皮的道地性研究也逐漸深入。采用藥物分析技術(shù)，研究人員不僅對(duì)其最佳采收期、不同產(chǎn)地來(lái)源藥材中主要活性成分丹皮酚和芍藥苷的含量進(jìn)行了研究[67]，也對(duì)不同產(chǎn)地“鳳丹”栽培后所獲得藥材的特征圖譜、QTOF圖譜、揮發(fā)油成分、無(wú)機(jī)元素等進(jìn)行了比較研究，探討藥材道地性的形成機(jī)制[8]。此外，利用生物技術(shù)研究土壤微生物的整體性和多樣性也逐步用于揭示中藥材牡丹皮的道地性[9]。

隨著中藥現(xiàn)代研究的逐步推進(jìn)，高通量測(cè)序（next generation sequencing，NGS）技術(shù)的不斷進(jìn)步與實(shí)驗(yàn)流程的不斷完善，藥用植物轉(zhuǎn)錄組研究成為一個(gè)新的研究熱點(diǎn)。自采用新一代高通量測(cè)序技術(shù)對(duì)二年生丹參根的轉(zhuǎn)錄組進(jìn)行測(cè)序研究以來(lái)[10]，陸續(xù)有30余種藥用植物開(kāi)展了轉(zhuǎn)錄組研究[11]。利用轉(zhuǎn)錄組測(cè)序技術(shù)（RNA sequencing，RNASeq）可以較好地應(yīng)用于捕捉道地藥材與非道地藥材不同樣品間差異表達(dá)的基因，從整體上對(duì)參與生長(zhǎng)發(fā)育和次生代謝基因進(jìn)行系統(tǒng)的研究，進(jìn)而更好地闡釋藥材道地性形成的分子機(jī)制[12]。一些研究者雖然對(duì)牡丹低溫刺激下轉(zhuǎn)錄組的變化進(jìn)行了研究[13]，而對(duì)于道地產(chǎn)區(qū)與非道地產(chǎn)區(qū)“鳳丹”根皮中差異表達(dá)的基因缺少相關(guān)研究。本研究利用第二代測(cè)序技術(shù)對(duì)不同產(chǎn)地栽培的 “鳳丹”根皮的轉(zhuǎn)錄組進(jìn)行了從頭組裝（de novo assembly），豐富藥用牡丹的基因數(shù)據(jù)庫(kù)；同時(shí)，利用比較轉(zhuǎn)錄組學(xué)，分析不同產(chǎn)地、五年生、供采收的“鳳丹”根皮部基因差異表達(dá)的富集情況，為藥用牡丹根皮中活性成分的生物合成研究以及藥材道地性品質(zhì)的形成提供生物學(xué)依據(jù)；通過(guò)SSR位點(diǎn)的信息分析，為開(kāi)展分子標(biāo)記輔助育種、基因工程技術(shù)選育創(chuàng)制新的藥用牡丹優(yōu)良品種奠定基礎(chǔ)。

1材料和方法

11樣品于藥材最佳采收期，分別從位于藥材牡丹皮的傳統(tǒng)道地產(chǎn)區(qū)安徽省銅陵市鳳凰山（TLDP）及其周邊地區(qū)蕪湖市南陵縣（WHDP） 的北京同仁堂安徽中藥材有限公司GAP栽培基地、非傳統(tǒng)道地產(chǎn)區(qū)亳州市（BZDP）的安徽濟(jì)人藥業(yè)有限公司牡丹皮GAP栽培基地采挖五年生藥用植物“鳳丹”的主根，迅速剝?nèi)∑涓ぶ靡旱欣鋬觯筠D(zhuǎn)入-80 ℃保存?zhèn)溆谩?/p>

12總RNA的提取、文庫(kù)構(gòu)建和轉(zhuǎn)錄組測(cè)序使用TRIzol@Reagent試劑盒、Plant RNA Purification Reagent試劑盒（Invitrogen，美國(guó)） 提取、純化根皮總RNA，瓊脂糖凝膠電泳檢測(cè)RNA完整性，Nanodrop 2000超微量分光光度計(jì)（Thmermo，美國(guó)）對(duì)所提總RNA的濃度和純度進(jìn)行檢測(cè)，Agilent 2100生物分析儀（Agilent Technology，美國(guó)）測(cè)定RIN值，以RNA條帶清晰、28/23S亮度大于18/16S，60≤RIN<80為質(zhì)量合格。為了獲得盡可能豐富的“鳳丹”根皮采收期轉(zhuǎn)錄組信息，選取不同栽培基地各6株“鳳丹”根皮RNA等量混合為1份，得到3份“鳳丹”根皮混合總RNA樣品后送美吉生物完成反轉(zhuǎn)合成cDNA、連接adaptor以及HiSeq 4000（Illumina，美國(guó)）上機(jī)測(cè)序等后續(xù)工作。

13測(cè)序數(shù)據(jù)過(guò)濾及轉(zhuǎn)錄組的組裝使用SeqPrep （https：//githubcom/jstjohn/SeqPrep）、Sickle （https：//githubcom/najoshi/sickle） 軟件將測(cè)序所得到的原始數(shù)據(jù)（raw reads）進(jìn)行過(guò)濾以去除rRNA、含接頭的低質(zhì)量reads以及含N比率超過(guò)10%的reads，使得Q30（堿基的測(cè)序錯(cuò)誤率小于01%）達(dá)到80%以上，獲得去掉接頭后的測(cè)序序列（clean reads）用于后續(xù)分析。

由于牡丹基因組測(cè)序尚未完成，無(wú)參考序列，因此使用軟件Trinity[14]進(jìn)行“鳳丹”根皮部的轉(zhuǎn)錄組從頭組裝（de novo assembly）：通過(guò)測(cè)序序列之間的重疊（overlap）信息組裝得到重疊群（contigs），依據(jù)序列的雙端信息（pairedend）信息和重疊群之間的相似性對(duì)其進(jìn)行聚類(lèi)整合和延長(zhǎng)，進(jìn)而組裝得到轉(zhuǎn)錄本（transcripts），并從中選取最長(zhǎng)的轉(zhuǎn)錄本作為非重復(fù)序列基因（unigenes）。Trinity參數(shù)設(shè)置為Kmer=25，序列延伸成contig時(shí)overlap為Kmer1。

差異表達(dá)基因的分析與篩選是建立在同一套參考基因的基礎(chǔ)上，因此對(duì)3個(gè)樣品得到的unigene數(shù)據(jù)通過(guò)cdhit聚類(lèi)去除冗余，利用TGIGL的聚類(lèi)組裝策略最終得到非冗余的“鳳丹”采收期根皮的Allunigenes數(shù)據(jù)庫(kù)。

14功能基因注釋對(duì)“鳳丹”根皮Allunigenes轉(zhuǎn)錄數(shù)據(jù)使用BlastX軟件（Evalue<1×10-5）對(duì)unigene序列與NR庫(kù)（nonredundant database，http：//wwwncbinlmnihgov/），COG數(shù)據(jù)庫(kù)（clusters of orthologous group，http：//www. ncbi. nlm. nih. gov/COG/）和SwissProt數(shù)據(jù)庫(kù)（http：//wwwuniprotorg/）比對(duì)，進(jìn)行功能注釋和分類(lèi)處理；再對(duì)unigene序列進(jìn)行GO（http：//www. geneontology. org/）功能注釋和分類(lèi)，并用軟件對(duì)GO注釋結(jié)果分類(lèi)作圖，并將unigene與KEGG數(shù)據(jù)（kyoto encyclopedia of genes and genomes，http：//www. genome. jp/kegg/）進(jìn)行比對(duì)，分析相關(guān)的代謝通路。

15基因的表達(dá)水平及表達(dá)差異分析Illumina HiSeq 4000測(cè)序平臺(tái)得到的reads一般較短、插入刪除錯(cuò)誤較少，因此選擇短序列比對(duì)軟件Bowtie （http：//bowtiebio. sourceforge. net/index. shtml）完成兩兩樣本間的轉(zhuǎn)錄組數(shù)據(jù)比對(duì)；利用bowtie的比對(duì)結(jié)果，使用軟件RSEM（http：//deweylab. biostat. wisc. edu/rsem/）計(jì)算3個(gè)樣品比對(duì)到每個(gè)基因上的read count數(shù)目，然后對(duì)其極性FRKM轉(zhuǎn)換，進(jìn)而獲得不同樣本基因的表達(dá)水平[15]；根據(jù)RSEM得到的gene read count，使用軟件edgeR（http：//www. bioconductor. org/packages/2. 12/bioc/html/edgeR. html）進(jìn)行樣品間的差異基因分析。

為了進(jìn)一步研究道地與非道地產(chǎn)區(qū)“鳳丹”根皮中差異表達(dá)基因的富集水平，使用軟件Goatools（http：//github. com/tanghaibao/GOatools）對(duì)樣品間差異表達(dá)基因進(jìn)行GO庫(kù)注釋的分類(lèi)統(tǒng)計(jì)和功能富集分析，并使用Fisher精確檢驗(yàn)法對(duì)差異表達(dá)基因富集的顯著性進(jìn)行評(píng)價(jià)；使用KOBAS（http：//kobas. cbi. pku. edu. cn/home. do）進(jìn)行KEGG pathway富集分析，使用Fisher精確檢驗(yàn)法對(duì)差異表達(dá)基因KEGG通路富集的顯著性進(jìn)行評(píng)價(jià)。

16SSR分析簡(jiǎn)單重復(fù)序列（simple sequence repeats，SSR）在真核生物基因組中的重復(fù)次數(shù)表現(xiàn)為個(gè)體間高度變異，數(shù)量豐富，具有廣泛的應(yīng)用性。因此，采用MISA（http：//pgrc. ipkgatersleben. de/misa/misa. html）對(duì)unigene進(jìn)行SSR檢測(cè)，檢測(cè)中對(duì)應(yīng)的單核苷酸（mononucleotide）、二核苷酸（dinucleotide）、三核苷酸（trinucleotide）、四核苷酸（tetranucleotide）、五核苷酸（pentanucleotide）以及六核苷酸（hexanucleotide）等SSR重復(fù)類(lèi)型的最少重復(fù)次數(shù)分別設(shè)置為10，6，5，5，5，5。

2結(jié)果與分析

21“鳳丹”根皮RNASeq測(cè)序及質(zhì)量評(píng)估利用Illumina HiSeq 4000高通量測(cè)序技術(shù)對(duì)3個(gè)產(chǎn)地來(lái)源的“鳳丹”根皮轉(zhuǎn)錄組進(jìn)行測(cè)序，得到raw reads，clean reads的數(shù)量及測(cè)序長(zhǎng)度等見(jiàn)表1。3個(gè)樣本的轉(zhuǎn)錄組clean reads中Q20，Q30均大于90%，GC量均在45%左右，表明測(cè)序質(zhì)量較好，數(shù)據(jù)可用于后續(xù)AllUnigenes數(shù)據(jù)庫(kù)的建立和基因功能注釋、分類(lèi)及樣本間異同性分析。

22“鳳丹”根皮轉(zhuǎn)錄本組裝與分析利用Trinity軟件對(duì)上述3份獲得的clean reads去除重復(fù)部分，并根據(jù)序列之間的重疊部分進(jìn)行混合組裝，并通過(guò)相似性比較后組裝獲得AllUnigenes數(shù)據(jù)庫(kù)的轉(zhuǎn)錄本72 997條，轉(zhuǎn)錄本總長(zhǎng)度達(dá)到71 528 823 bp，最短的僅為201 bp，最長(zhǎng)的達(dá)到18 123 bp，平均長(zhǎng)度78207 bp，長(zhǎng)度在200～400 bp的有40 123條，占4387%；獲得unigene72 997條，總長(zhǎng)度達(dá)51 290 894 bp，平均長(zhǎng)度70264 bp。其中，長(zhǎng)度在200～400 bp的有37 320條，占5113%；轉(zhuǎn)錄本和unigene的N50分別達(dá)到1 135，1 227 bp，均大于800 bp，表明組裝片段的完整性較高，見(jiàn)圖1。

23“鳳丹”根皮轉(zhuǎn)錄組AllUnigenes基因功能注釋及分類(lèi)經(jīng)過(guò) BlastX 序列比對(duì)分析，AllUnigenes在NR，COG，GO，KEGG數(shù)據(jù)庫(kù)中獲得注釋unigene的數(shù)量分別為41 139，8 627，24 952，20 016。其中，NR 數(shù)據(jù)庫(kù)中搜索unigene 相似序列最多，為41 139條（564%）；通過(guò)COG注釋?zhuān)? 672條unigene獲得了分類(lèi)信息，這些基因分屬于25個(gè)功能分類(lèi)，見(jiàn)圖2。1 063條unigene屬一般功能預(yù)測(cè)項(xiàng)（1226%），是獲得注釋數(shù)量最多的部分；1 057條unigene（1219%）與翻譯、核糖體結(jié)構(gòu)與生源相關(guān)。此外，還有394條（247%）功能未知的unigene基因可能是“鳳丹”特有的基因，與其供藥用的所具有的發(fā)揮藥效的次生代謝物質(zhì)有關(guān)。

GO注釋?zhuān)涸诳偣步M裝的72 997條unigene中，有24 952條unigene能被GO 數(shù)據(jù)庫(kù)成功注釋?zhuān)譃樯飳W(xué)過(guò)程（biological process）、細(xì)胞成分（cellular component）及分子功能（molecular function）三大類(lèi)

ARNA 加工與修飾；B染色質(zhì)結(jié)構(gòu)與動(dòng)力學(xué)；C能量生產(chǎn)和轉(zhuǎn)換；D細(xì)胞周期控制、細(xì)胞分裂、染色體分區(qū)；E氨基酸運(yùn)輸和代謝；F核苷酸運(yùn)輸和代謝；G碳水化合物的運(yùn)輸和代謝；H輔酶運(yùn)輸和代謝；I脂質(zhì)運(yùn)輸和代謝；J翻譯；核糖體結(jié)構(gòu)和生物轉(zhuǎn)化；K轉(zhuǎn)錄；L復(fù)制、重組和修復(fù)；M細(xì)胞被膜生物起源；N細(xì)胞運(yùn)動(dòng)性；O翻譯后修飾，蛋白質(zhì)轉(zhuǎn)換與分子伴侶；P無(wú)機(jī)離子轉(zhuǎn)運(yùn)和代謝；Q次生代謝產(chǎn)物合成、運(yùn)輸與分解；R一般功能預(yù)測(cè)；S未知功能；T信號(hào)傳導(dǎo)機(jī)制；U胞內(nèi)運(yùn)輸、分泌及膜泡轉(zhuǎn)運(yùn)；V防御機(jī)制；W細(xì)胞外結(jié)構(gòu)；Y核結(jié)構(gòu)；Z細(xì)胞骨架。

共60個(gè)功能組。進(jìn)一步研究發(fā)現(xiàn)： 24 952條unigene在GO庫(kù)中歸屬于生物學(xué)過(guò)程中的24個(gè)功能組。其中，涉及代謝過(guò)程的unigene17 113條、細(xì)胞過(guò)程14 817條、單生物過(guò)程11 975條、定位3 413條、定位建立3 350條、生物調(diào)節(jié)3 297條、應(yīng)激反應(yīng)3 202條、生物學(xué)過(guò)程的調(diào)控3 098條，以參與代謝過(guò)程和細(xì)胞過(guò)程獲得序列注釋最多；同時(shí)，另有31 331條 unigene在GO庫(kù)中歸屬于分子功能的16個(gè)功能組，其中獲得注釋較多的unigene分別為催化活性13 589條、結(jié)合功能13 207條，這些基因可能與牡丹皮中化學(xué)成分的催化與合成有關(guān)。此外，44 873條unigenen在GO庫(kù)中歸屬于細(xì)胞成分的20個(gè)功能組。其中，被注釋到較多的unigene分別為細(xì)胞9 669條、細(xì)胞器6 737條、細(xì)胞膜5 500條和復(fù)雜大分子4 363條。

KEGG注釋?zhuān)汗?0 016條unigene能夠在KEGG數(shù)據(jù)庫(kù)成功注釋?zhuān)婕按x（metabolism）、遺傳信息處理（genetic information processing）、環(huán)境信息處理（environmental information processing）、細(xì)胞進(jìn)程（cellular processes）及生物系統(tǒng)（organismal systems）5個(gè)大類(lèi)及34個(gè)中類(lèi)、352條代謝通路，其中代謝大類(lèi)中有6 469條unigene被成功注釋?zhuān)浯问谴紊x產(chǎn)物的合成（biosynthesis of secondary metabolites）2 481條和不同環(huán)境的微生物代謝（microbial metabolism in diverse environment）1 735條。進(jìn)一步研究還發(fā)現(xiàn)，在代謝途徑中，與中藥效活性成分生物合成相關(guān)的unigene序列分別是苯丙素類(lèi)合成（phenylpropanoid biosynthesis）214條、萜類(lèi)化合物骨架合成（terpenoid backbone biosynthesis）104條、各種萜類(lèi)化合物合成（sesquiterpenoid，triterpenoid，monoterpenoid and diterpenoid biosynthesis）152條、（莨菪類(lèi)、哌啶類(lèi)、吡啶類(lèi)）生物堿合成（tropane，piperidine and pyridine alkaloid biosynthesis）55條、黃酮類(lèi)成分合成（flavonoid，flavone，flavonol and isoflavonoid biosynthesis）36條。

24差異表達(dá)基因的分析3個(gè)栽培產(chǎn)地樣品間“鳳丹”根皮總AllUnigenes比對(duì)的結(jié)果分別為BZDP（8318%），TLDP（8178%），WHDP（8621%），不同樣品的轉(zhuǎn)錄本/基因的表達(dá)量概率（FPKM scores）密度分布見(jiàn)圖3，其中以TLDP基因的密度最高；顯著差異表達(dá)基因的篩選標(biāo)準(zhǔn)為FDR<005且log2|FC|≥1，3個(gè)樣品在72 997條根皮AllUnigenes轉(zhuǎn)錄本數(shù)據(jù)庫(kù)能被識(shí)別到的差異基因總數(shù)為3 430條，不同樣品間分組BZDP vs. TLDP，WHDP vs. TLDP，BZDP vs. WHDP的表達(dá)差異基因數(shù)分別為 1 502，878，1 050，各組間差異基因的數(shù)量分布狀況見(jiàn)圖4。可以看出，非道地產(chǎn)區(qū)樣品與傳統(tǒng)道地產(chǎn)區(qū)樣品（BZDP vs. TLDP）間顯著差異表達(dá)基因數(shù)目最多，傳統(tǒng)道地產(chǎn)區(qū)銅陵及其周邊南陵地區(qū)樣品（WHDP vs. TLDP）差異表達(dá)基因數(shù)目最少；具有顯著性富集（P<0001）的差異基因功能類(lèi)型和GO庫(kù)的ID號(hào)見(jiàn)表2。

由表2可以看出，BZDPTLDP間差異基因顯著富集于生物學(xué)過(guò)程和分子功能的9個(gè)功能組，BZDPWHDP間差異基因顯著富集于生物學(xué)過(guò)程、分子功能和細(xì)胞組成的22個(gè)功能組，TLDPWHDP間差異基因顯著富集于生物學(xué)過(guò)程、分子功能和細(xì)胞組成的20個(gè)功能組。不同樣品間差異基因在KEGG通路富集具有特別顯著性差異（P<0001）的通路描述和ID號(hào)見(jiàn)表3。可以看出BZDPTLDP間差異基因顯著富集主要是遺傳信息和代謝的8個(gè)代謝通路，BZDPWHDP間差異基因顯著富集于遺傳信息、生物系統(tǒng)、代謝和疾病的8個(gè)代謝通路，TLDPWHDP間差異基因顯著富集于代謝和疾病的6個(gè)代謝通路。

25SSR分子標(biāo)記在72 997條AllUnigene庫(kù)中，共檢測(cè)到9 939條unigene中含有SSR，發(fā)生頻率（含有SSR的unigene數(shù)量與總unigene數(shù)量之比）為136%。其中6 060條unigene序列只含有1個(gè)SSR位點(diǎn)，518條為復(fù)雜重復(fù)類(lèi)型SSR，含有SSR標(biāo)記的unigene序列數(shù)量及類(lèi)型統(tǒng)計(jì)見(jiàn)表4。

統(tǒng)計(jì)結(jié)果顯示，除1個(gè)SSR位點(diǎn)外，“鳳丹”根皮AllUnigenes轉(zhuǎn)錄組中SSR的主要類(lèi)型是二核苷酸，占SSR總數(shù)的2075%，出現(xiàn)頻率最高的3個(gè)重復(fù)類(lèi)型是CT/GA，AG/TC，AT/TA；其次是三核苷酸，占SSR總數(shù)的1223%；四、五、六核苷酸重復(fù)類(lèi)型的數(shù)量很少，占SSR總數(shù)的085%；SSR位點(diǎn)的序列分布從10～219 bp不等，平均長(zhǎng)度17 bp。

3討論與結(jié)論

中藥丹參[10]、人參[16]和虎杖[17]等的轉(zhuǎn)錄組測(cè)序多數(shù)是利用藥用部位根進(jìn)行，分別獲得了18 235，31 741，86 418條unigenes，為最大限度的獲得“鳳

丹”根皮的轉(zhuǎn)錄組信息，研究選擇安徽省傳統(tǒng)道地產(chǎn)區(qū)和非道地產(chǎn)區(qū)、五年生“鳳丹”根皮為試驗(yàn)材料，進(jìn)行混合樣品測(cè)序，共組裝得到72 997條“鳳丹”根皮的基因序列，加快了藥用植物“鳳丹”的分子生物學(xué)研究。

經(jīng)過(guò) BlastX 比對(duì)分析，436%的unigene在NR庫(kù)中不能匹配到已知基因，可能由于NR 數(shù)據(jù)庫(kù)中由于缺少牡丹基因組信息；COG分類(lèi)與注釋中，發(fā)現(xiàn)有318條unigene（367%）參與了次生代謝物質(zhì)的生物合成、轉(zhuǎn)運(yùn)和降解，深入研究這些基因的功能將有助于揭示“鳳丹”根皮中有效成分合成的調(diào)控機(jī)制；KEGG注釋中，與次生代謝產(chǎn)物的合成途徑相關(guān)的unigene有2 481條，與不同環(huán)境的微生物代謝相關(guān)的有1 735條，說(shuō)明五年生“鳳丹”的植物根皮中參與次生代謝產(chǎn)物合成與代謝的基因豐富，為進(jìn)一步研究“鳳丹”道地藥材品質(zhì)的形成提供了依據(jù)。來(lái)自于非道地產(chǎn)區(qū)樣品與傳統(tǒng)道地產(chǎn)區(qū)“鳳丹”根皮中差異表達(dá)基因數(shù)目的多少可以較好地反映在

藥用種質(zhì)遺傳背景一致的情況下，道地產(chǎn)區(qū)的環(huán)境、氣候?qū)τ谒幉牡赖匦云焚|(zhì)的形成具有較大的影響，傳統(tǒng)道地產(chǎn)區(qū)銅陵及其相鄰地區(qū)南陵因生態(tài)環(huán)境、氣候相近，兩地來(lái)源樣本間在次生代謝物合成中差異表達(dá)基因的數(shù)量遠(yuǎn)遠(yuǎn)少于環(huán)境、氣候距離較遠(yuǎn)的亳州栽培“鳳丹”根皮樣本。