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

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

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

[關(guān)鍵詞]牡丹皮； 轉(zhuǎ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 improvementendprint

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

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

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

隨著中藥現(xiàn)代研究的逐步推進(jìn)，高通量測序（next generation sequencing，NGS）技術(shù)的不斷進(jìn)步與實驗流程的不斷完善，藥用植物轉(zhuǎn)錄組研究成為一個新的研究熱點。自采用新一代高通量測序技術(shù)對二年生丹參根的轉(zhuǎn)錄組進(jìn)行測序研究以來[10]，陸續(xù)有30余種藥用植物開展了轉(zhuǎn)錄組研究[11]。利用轉(zhuǎn)錄組測序技術(shù)（RNA sequencing，RNASeq）可以較好地應(yīng)用于捕捉道地藥材與非道地藥材不同樣品間差異表達(dá)的基因，從整體上對參與生長發(fā)育和次生代謝基因進(jìn)行系統(tǒng)的研究，進(jìn)而更好地闡釋藥材道地性形成的分子機(jī)制[12]。一些研究者雖然對牡丹低溫刺激下轉(zhuǎn)錄組的變化進(jìn)行了研究[13]，而對于道地產(chǎn)區(qū)與非道地產(chǎn)區(qū)“鳳丹”根皮中差異表達(dá)的基因缺少相關(guān)研究。本研究利用第二代測序技術(shù)對不同產(chǎn)地栽培的 “鳳丹”根皮的轉(zhuǎn)錄組進(jìn)行了從頭組裝（de novo assembly），豐富藥用牡丹的基因數(shù)據(jù)庫；同時，利用比較轉(zhuǎn)錄組學(xué)，分析不同產(chǎn)地、五年生、供采收的“鳳丹”根皮部基因差異表達(dá)的富集情況，為藥用牡丹根皮中活性成分的生物合成研究以及藥材道地性品質(zhì)的形成提供生物學(xué)依據(jù)；通過SSR位點的信息分析，為開展分子標(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)∑涓ぶ靡旱欣鋬?，后轉(zhuǎn)入-80 ℃保存?zhèn)溆谩?/p>

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

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

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

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

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

15基因的表達(dá)水平及表達(dá)差異分析Illumina HiSeq 4000測序平臺得到的reads一般較短、插入刪除錯誤較少，因此選擇短序列比對軟件Bowtie （http：//bowtiebio. sourceforge. net/index. shtml）完成兩兩樣本間的轉(zhuǎn)錄組數(shù)據(jù)比對；利用bowtie的比對結(jié)果，使用軟件RSEM（http：//deweylab. biostat. wisc. edu/rsem/）計算3個樣品比對到每個基因上的read count數(shù)目，然后對其極性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）對樣品間差異表達(dá)基因進(jìn)行GO庫注釋的分類統(tǒng)計和功能富集分析，并使用Fisher精確檢驗法對差異表達(dá)基因富集的顯著性進(jìn)行評價；使用KOBAS（http：//kobas. cbi. pku. edu. cn/home. do）進(jìn)行KEGG pathway富集分析，使用Fisher精確檢驗法對差異表達(dá)基因KEGG通路富集的顯著性進(jìn)行評價。

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

2結(jié)果與分析

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

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

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

GO注釋：在總共組裝的72 997條unigene中，有24 952條unigene能被GO 數(shù)據(jù)庫成功注釋，分為生物學(xué)過程（biological process）、細(xì)胞成分（cellular component）及分子功能（molecular function）三大類

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

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

KEGG注釋：共20 016條unigene能夠在KEGG數(shù)據(jù)庫成功注釋，涉及代謝（metabolism）、遺傳信息處理（genetic information processing）、環(huán)境信息處理（environmental information processing）、細(xì)胞進(jìn)程（cellular processes）及生物系統(tǒng)（organismal systems）5個大類及34個中類、352條代謝通路，其中代謝大類中有6 469條unigene被成功注釋，其次是次生代謝產(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序列分別是苯丙素類合成（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差異表達(dá)基因的分析3個栽培產(chǎn)地樣品間“鳳丹”根皮總AllUnigenes比對的結(jié)果分別為BZDP（8318%），TLDP（8178%），WHDP（8621%），不同樣品的轉(zhuǎn)錄本/基因的表達(dá)量概率（FPKM scores）密度分布見圖3，其中以TLDP基因的密度最高；顯著差異表達(dá)基因的篩選標(biāo)準(zhǔn)為FDR<005且log2|FC|≥1，3個樣品在72 997條根皮AllUnigenes轉(zhuǎn)錄本數(shù)據(jù)庫能被識別到的差異基因總數(shù)為3 430條，不同樣品間分組BZDP vs. TLDP，WHDP vs. TLDP，BZDP vs. WHDP的表達(dá)差異基因數(shù)分別為 1 502，878，1 050，各組間差異基因的數(shù)量分布狀況見圖4?？梢钥闯觯堑赖禺a(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）的差異基因功能類型和GO庫的ID號見表2。

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

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

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

3討論與結(jié)論

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

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

經(jīng)過 BlastX 比對分析，436%的unigene在NR庫中不能匹配到已知基因，可能由于NR 數(shù)據(jù)庫中由于缺少牡丹基因組信息；COG分類與注釋中，發(fā)現(xiàn)有318條unigene（367%）參與了次生代謝物質(zhì)的生物合成、轉(zhuǎn)運和降解，深入研究這些基因的功能將有助于揭示“鳳丹”根皮中有效成分合成的調(diào)控機(jī)制；KEGG注釋中，與次生代謝產(chǎn)物的合成途徑相關(guān)的unigene有2 481條，與不同環(huán)境的微生物代謝相關(guān)的有1 735條，說明五年生“鳳丹”的植物根皮中參與次生代謝產(chǎn)物合成與代謝的基因豐富，為進(jìn)一步研究“鳳丹”道地藥材品質(zhì)的形成提供了依據(jù)。來自于非道地產(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)境、氣候相近，兩地來源樣本間在次生代謝物合成中差異表達(dá)基因的數(shù)量遠(yuǎn)遠(yuǎn)少于環(huán)境、氣候距離較遠(yuǎn)的亳州栽培“鳳丹”根皮樣本。endprint

轉(zhuǎn)錄組測序技術(shù)的應(yīng)用同時也加速了 SSR 位點的挖掘。本研究中“鳳丹”根皮的AllUnigenes轉(zhuǎn)錄本中SSR出現(xiàn)的頻率為136%，低于連翹（2515%）[18]、半夏（1624%）[18]和獨行菜（1564%）[20]，但高于黨參（1222%）[21]、膜莢黃芪（924%）[22]、魚腥草（751%）[2324]、人參（73%）[16]和杜仲（29%）[25]等藥用植物中SSR出現(xiàn)的頻率，可能與物種以及供試材料的來源部位多少有關(guān)；“鳳丹”根皮轉(zhuǎn)錄組中SSR種類較為豐富，主要以單核苷酸重復(fù)為主，其次是二、三核苷酸重復(fù)，二核苷酸重復(fù)單元中AG/TC（681條，占SSR總數(shù)的685%）占主導(dǎo)地位；此外，三核苷酸的重復(fù)單元類型較為豐富，這與多數(shù)植物SSR標(biāo)記的特性相似；與基于NCBI 數(shù)據(jù)庫的觀賞牡丹ESTSSR研究結(jié)果[26]不同的是，本研究中發(fā)現(xiàn)了“鳳丹”根皮轉(zhuǎn)錄本中有11條unigene序列中存在五核苷酸重復(fù)，可能與“鳳丹”根皮供藥用的活性成分代謝與合成有關(guān)，為深入挖掘并利用“鳳丹”轉(zhuǎn)錄組數(shù)據(jù)開發(fā)用于遺傳多樣性分析、分子標(biāo)記輔助育種的SSR標(biāo)記奠定了基礎(chǔ)。

藥用植物的主要有效成分是次生代謝產(chǎn)物，其合成與代謝途徑一般受到多個基因的調(diào)控，通過與NR，COG，GO，KEGG等公共數(shù)據(jù)庫進(jìn)行序列比對，注釋的結(jié)果將有望大大推動道地藥材品質(zhì)的形成以及根皮中有效成分合成途徑的相關(guān)基因提供了依據(jù)；完成“鳳丹”根皮的轉(zhuǎn)錄組測序僅是優(yōu)質(zhì)藥用植物種質(zhì)資源功能基因組研究的第一步，只有對藥用植物的基因功能分化進(jìn)行深入分析，并與根皮部活性成分種類、數(shù)量和積累變化相結(jié)合，進(jìn)而從分子生藥學(xué)角度揭示由于不同生長時間和生長環(huán)境基因表達(dá)的差異與藥用部位中次生代謝成分差異的機(jī)制，才能為藥材質(zhì)量控制與提高提供理論依據(jù)。

[參考文獻(xiàn)]

[1]中國藥典一部[S]2015：172.

[2]王祝舉，唐力英，赫炎牡丹皮的化學(xué)成分和藥理作用[J]. 國外醫(yī)藥·植物藥分冊，2006，21（4）：155.

[3]胡云飛，許國兵牡丹皮及其主要成分丹皮酚的藥理作用研究進(jìn)展[J]. 安徽醫(yī)藥，2014，18（4）：589.

[4]方成武安徽道地藥材牡丹皮的采收及產(chǎn)地加工方法考察[J]. 中藥材，2000，23（2）：82.

[5]侯宇榮，劉煒，鄭艷安徽南陵丫山產(chǎn)牡丹皮的道地性研究[J]. 中藥材，2014，37（8）：1488.

[6]孔昭琰，巢建國，谷巍，等不同采收期鳳丹皮中3種藥效成分動態(tài)變化研究[J]. 中國中醫(yī)藥信息雜志，2011，18（1）：63.

[7]紀(jì)開明，方成武，管玉云，等銅陵和亳州牡丹皮中丹皮酚和芍藥苷的含量比較研究[J]. 中醫(yī)藥導(dǎo)報，2014，20（10）：78.

[8]胡云飛基于藥物分析組合技術(shù)研究鳳丹藥材的道地性[D]. 合肥：安徽中醫(yī)藥大學(xué)，2015.

[9]鄭艷，劉煒，黃軍祥，等基于牡丹根際土壤微生物的中藥材道地性研究[J]. 藥學(xué)學(xué)報，2016，51（8）：1325.

[10]李瀅，孫超，羅紅梅基于高通量測序454 GS FLX的丹參轉(zhuǎn)錄組學(xué)研究[J]. 藥學(xué)學(xué)報， 2010，45（4）：524.

[11]王堯龍，黃璐琦，袁媛，等藥用植物轉(zhuǎn)錄組研究進(jìn)展[J]. 中國中藥雜志，2015，40（11）：2055.

[12]王曉玥，宋經(jīng)元，謝彩香，等RNASeq與道地藥材研究[J]. 藥學(xué)學(xué)報，2014，49（12）：1650.

[13]Gai S，Zhang Y，Mu P，et al. Transcriptome analysis of tree peony during chilling requirement fulfillment： assembling， annotation and markers discovering[J]. Gene，2012， 497（2）：256.

[14]Grabherr M G， Haas B J，Yassour M，et al. Fulllength transcriptome assembly from RNAseq data without a reference genome[J]. Nat Biotechnol，2011，29（7）：644.

[15]Li B，Dewey C N RSEM： accurate transcript quantification from RNASeq data with or without a reference genome[J]. BMC Bioinformatics，2011，12（1）：93.

[16]Li C，Zhu Y，Guo X，et al. Transcriptome analysis reveals ginsenosides biosynthetic genes， microRNAs and simple sequence repeats in Panax ginseng CA Meyer[J]. BMC Genomics，2013，14（1）： 245.

[17]郝大程，馬培，穆軍，等中藥植物虎杖根的高通量轉(zhuǎn)錄組測序及轉(zhuǎn)錄組特性分析[J]. 中國科學(xué)：生命科學(xué)，2012，42（5）：398.

[18]王興春，譚河林，陳釗，等基于RNA Seq技術(shù)的連翹轉(zhuǎn)錄組組裝與分析及SSR分子標(biāo)記的開發(fā)[J]. 中國科學(xué)：生命科學(xué)，2015，45（3）：301.

[19]王森，張震，姜倪皓，等半夏轉(zhuǎn)錄組中的SSR位點信息分析[J]. 中藥材，37（9）：1566.

[20]周茜，趙惠新，李萍萍，等獨行菜種子轉(zhuǎn)錄組的高通量測序及分析[J]. 中國生物工程雜志，2016，36（1）：38.

[21]王東，曹玲亞，高建平黨參轉(zhuǎn)錄組中SSR位點信息分析[J].中草藥， 2016，45（16）：2390.

[22]常越，閆嵩，劉振鵬，等膜莢黃芪的轉(zhuǎn)錄組測序質(zhì)量評估及其SSR位點信息分析的研究[J]. 中國中藥雜志，2016，41（8）：1430.

[23]Wei L，Li S， Liu S，et al. Transcriptome analysis of Houttuynia cordata Thunb by Illumina pairedend RNA sequencing and SSR marker discovery[J]. PLoS ONE，2014，9（1）： e84105.

[24]黎曉英，劉勝貴，王丹，等魚腥草轉(zhuǎn)錄組SSR位點信息分析及其多態(tài)性研究[J]. 中草藥，2016，47（10）：1762.

[25]黃海燕，杜紅巖，烏云塔娜，等基于杜仲轉(zhuǎn)錄組序列的 SSR 分子標(biāo)記的開發(fā)[J]. 林業(yè)科學(xué)， 2013，49（5）：176.

[26]侯小改，王娟，郭大龍，等牡丹EST資源的SSR信息分析[J]. 湖南農(nóng)業(yè)大學(xué)學(xué)報：自然科學(xué)版，2011，37（2）：172

[責(zé)任編輯呂冬梅]endprint