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

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

[摘要]牡丹皮為我國傳統常用中藥，安徽省銅陵地區栽培的“鳳丹”根皮加工而成的藥材牡丹皮被譽為道地藥材，藥用活性成分豐富多樣，但目前尚不清楚 “鳳丹”藥用部位次生代謝過程中活性物質合成的遺傳學基礎。研究采用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 improvementendprint

[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/）進行比對，分析相關的代謝通路。endprint

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條。endprint

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條，說明五年生“鳳丹”的植物根皮中參與次生代謝產物合成與代謝的基因豐富，為進一步研究“鳳丹”道地藥材品質的形成提供了依據。來自于非道地產區樣品與傳統道地產區“鳳丹”根皮中差異表達基因數目的多少可以較好地反映在

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

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

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

[參考文獻]

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

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

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

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

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

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

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

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

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

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

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

[12]王曉玥，宋經元，謝彩香，等RNASeq與道地藥材研究[J]. 藥學學報，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]郝大程，馬培，穆軍，等中藥植物虎杖根的高通量轉錄組測序及轉錄組特性分析[J]. 中國科學：生命科學，2012，42（5）：398.

[18]王興春，譚河林，陳釗，等基于RNA Seq技術的連翹轉錄組組裝與分析及SSR分子標記的開發[J]. 中國科學：生命科學，2015，45（3）：301.

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

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

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

[22]常越，閆嵩，劉振鵬，等膜莢黃芪的轉錄組測序質量評估及其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]黎曉英，劉勝貴，王丹，等魚腥草轉錄組SSR位點信息分析及其多態性研究[J]. 中草藥，2016，47（10）：1762.

[25]黃海燕，杜紅巖，烏云塔娜，等基于杜仲轉錄組序列的 SSR 分子標記的開發[J]. 林業科學， 2013，49（5）：176.

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

[責任編輯呂冬梅]endprint