基于高通量測序的藥用植物“鳳丹”根皮的轉錄組分析

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

[摘要]牡丹皮為我國傳統常用中藥，安徽省銅陵地區栽培的“鳳丹”根皮加工而成的藥材牡丹皮被譽為道地藥材，藥用活性成分豐富多樣，但目前尚不清楚 “鳳丹”藥用部位次生代謝過程中活性物質合成的遺傳學基礎。研究采用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?？梢钥闯觯堑赖禺a區樣品與傳統道地產區樣品（BZDP vs. TLDP）間顯著差異表達基因數目最多，傳統道地產區銅陵及其周邊南陵地區樣品（WHDP vs. TLDP）差異表達基因數目最少；具有顯著性富集（P<0001）的差異基因功能類型和GO庫的ID號見表2。

由表2可以看出，BZDPTLDP間差異基因顯著富集于生物學過程和分子功能的9個功能組，BZDPWHDP間差異基因顯著富集于生物學過程、分子功能和細胞組成的22個功能組，TLDPWHDP間差異基因顯著富集于生物學過程、分子功能和細胞組成的20個功能組。不同樣品間差異基因在KEGG通路富集具有特別顯著性差異（P<0001）的通路描述和ID號見表3?？梢钥闯鯞ZDPTLDP間差異基因顯著富集主要是遺傳信息和代謝的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條，說明五年生“鳳丹”的植物根皮中參與次生代謝產物合成與代謝的基因豐富，為進一步研究“鳳丹”道地藥材品質的形成提供了依據。來自于非道地產區樣品與傳統道地產區“鳳丹”根皮中差異表達基因數目的多少可以較好地反映在

藥用種質遺傳背景一致的情況下，道地產區的環境、氣候對于藥材道地性品質的形成具有較大的影響，傳統道地產區銅陵及其相鄰地區南陵因生態環境、氣候相近，兩地來源樣本間在次生代謝物合成中差異表達基因的數量遠遠少于環境、氣候距離較遠的亳州栽培“鳳丹”根皮樣本。