藥用植物冬凌草高通量轉(zhuǎn)錄組測(cè)序與分析△

陳延清，胡志剛，劉大會(huì)，黃必勝，劉迪

(湖北中醫(yī)藥大學(xué) 藥學(xué)院，湖北 武漢 430065)

大多數(shù)唇形科藥用植物味甘苦，性微寒，中國(guó)民間習(xí)用于治療失眠、消化不良、風(fēng)濕病、腫瘤等[1-2]?，F(xiàn)代化學(xué)和藥理研究表明，唇形科植物主要的生物活性成分為二萜化合物，由于其復(fù)雜多樣的碳骨架和具有多種抗腫瘤藥理活性成分，二萜類(lèi)化合物已經(jīng)成為抗癌藥物候選藥物[3-4]。冬凌草Isodonrubescens(Hemsl.)H.Hara屬于唇形科中藥用價(jià)值和經(jīng)濟(jì)價(jià)值較為重要的藥用植物，其中全草中含量較高二萜類(lèi)化合物冬凌草甲素、冬凌草乙素具有良好殺菌消炎、抗腫瘤、抗氧化等作用。目前，冬凌草中的二萜類(lèi)化合物的研究與市場(chǎng)開(kāi)發(fā)成為了熱點(diǎn)[5-6]。隨著對(duì)冬凌草資源與品質(zhì)的研究，其道地性和種質(zhì)資源的遺傳規(guī)律探索也逐漸深入[7]。采用現(xiàn)代化學(xué)分析技術(shù)，對(duì)冬凌草最佳采收期[8]，以及對(duì)不同產(chǎn)地來(lái)源主要二萜活性成分冬凌草甲素和冬凌草乙素的含量差異進(jìn)行了研究[9]，并建立了不同產(chǎn)地冬凌草藥材的HPLC-PDA指紋圖譜特征圖譜[10]，闡述冬凌草藥材的最佳采收時(shí)期和產(chǎn)地差異。利用隨機(jī)引物標(biāo)記RAPD和ISSR[11]，以及特異標(biāo)記簡(jiǎn)單序列重復(fù)SSR(simple sequence repeat，SSR)[12]等分子標(biāo)記技術(shù)，探索了冬凌草遺傳多樣性規(guī)律，為篩選優(yōu)質(zhì)種質(zhì)資源提供了豐富的遺傳基礎(chǔ)，探討藥材道地性的形成機(jī)制。

高通量測(cè)序(next generation sequencing，NGS)技術(shù)是對(duì)傳統(tǒng)測(cè)序的革命性改變，轉(zhuǎn)錄組學(xué)研究方法轉(zhuǎn)錄組測(cè)序(RNA-seq)揭示藥用植物遺傳育種和基因進(jìn)化規(guī)律，其具有通量高、成本低、不受基因組序列信息限制等優(yōu)點(diǎn)[13]。對(duì)于非模式生物而言，例如冬凌草，轉(zhuǎn)錄組測(cè)序更偏重基因編碼區(qū)域，與龐大的基因組信息相比，重復(fù)元件和富含GC的區(qū)域比較少，拼接相對(duì)容易[14]，轉(zhuǎn)錄組測(cè)序針對(duì)物種轉(zhuǎn)錄本直接進(jìn)行測(cè)序，為基因組序列信息有限的非模式生物提供可能。本研究利用 Illumina Hiseq 4000對(duì)藥用植物冬凌草藥用部位葉和莖進(jìn)行了高通量轉(zhuǎn)錄組測(cè)序，建立冬凌草轉(zhuǎn)錄組數(shù)據(jù)庫(kù)，為深入研究冬凌草中冬凌草甲素等有效成分的生物合成途徑及其調(diào)控機(jī)制提供基礎(chǔ)，為藥用活性成分的合成代謝以及藥材道地性品質(zhì)的形成提供生物學(xué)依據(jù)。高效地發(fā)現(xiàn)非模式物種的新基因和SSR位點(diǎn)的信息分析，也為開(kāi)展分子標(biāo)記輔助育種、基因工程技術(shù)選育創(chuàng)制新的藥用冬凌草優(yōu)良品種奠定基礎(chǔ)，對(duì)優(yōu)質(zhì)中藥材的品種選育和規(guī)范化栽培生產(chǎn)具有重要的指導(dǎo)意義。

1 材料與方法

1.1 材料

本研究的冬凌草樣品來(lái)自于道地產(chǎn)區(qū)河南濟(jì)源[15]，采于藥材最佳采收期[8]，在同一生境下，選取冬凌草地上部分葉和莖為研究材料，迅速用液氮研磨，用于 RNA樣品提取。

1.2 總RNA的提取文庫(kù)、構(gòu)建和轉(zhuǎn)錄組測(cè)序

用 TRIzol?Reagent、試劑盒Plant RNA Purification Reaget試劑盒(Invitrogen，美國(guó))提取、純化葉與莖總RNA[16]。分別用Nanodrop 2000與 Agilent 2100 Bioanalyzer進(jìn)行總RNA質(zhì)量和純度檢測(cè)。為了獲得盡可能豐富的冬凌草采收期轉(zhuǎn)錄組信息，特別地將同珠的葉和莖的RNA混合。富集mRNA使用 Dynabeads Oligo(dT)25的磁珠。在Thermomixer中加入fragmentation buffer將mRNA打斷為短片段，添加0.1 mol·L-1Mg2+緩沖液到所得的 mRNA 中，高溫下使其片斷化后，然后將 mRNA為模板反轉(zhuǎn)合成cDNA、然后通過(guò)磁珠純化、末端修復(fù)、3’末端加堿基 A、加測(cè)序接頭等步驟后選擇片段大小，進(jìn)行PCR擴(kuò)增，adaptor構(gòu)建冬凌草轉(zhuǎn)錄組測(cè)序文庫(kù)，使用Illumina HiSeq 4000 進(jìn)行測(cè)序及分析。

1.3 測(cè)序數(shù)據(jù)質(zhì)量處理及組裝

使用SeqPrep(https：//github.com/jstjohn/SeqPrep)、Sickle(https：//github.com/najoshi/sickle)軟件進(jìn)行過(guò)濾原始數(shù)據(jù)，從而去除含有Primer/Adaptor的序列、含N比率超過(guò)10%的reads和含接頭的低質(zhì)量reads，并且使Q30(堿基測(cè)序錯(cuò)誤率低于0.1%)達(dá)到80%以上。從而獲得去掉接頭后的測(cè)序序列(cleanreads)進(jìn)行后續(xù)分析。

由于缺乏冬凌草基因組的數(shù)據(jù)，因此利用Trinity軟件對(duì)冬凌草葉與莖轉(zhuǎn)錄組重頭組裝(de novo assembly)，測(cè)序序列之間的重疊(overlap)信息組裝得到重疊群(contigs)，依據(jù)兩端信息(paired-end)和重疊群的相似性，然后使用TGICL軟件進(jìn)行聚類(lèi)整合和延長(zhǎng)得到非冗余的轉(zhuǎn)錄本(transcripts)。并從中選取最長(zhǎng)的轉(zhuǎn)錄本作為非重復(fù)序列基因(unigenes)。Trinity 參數(shù)設(shè)置為K-mer=25，當(dāng)序列延伸成重疊群時(shí)重復(fù)為 K-mer-1。拼接完畢后，將沒(méi)有匹配結(jié)果的短讀序的轉(zhuǎn)錄本刪除。處理與注釋冬凌草葉和莖的差異表達(dá)基因使用相同參考基因，得到同一時(shí)期冬凌草葉和莖所有Unigenes數(shù)據(jù)庫(kù)使用TGIGL的聚類(lèi)組裝策略。

1.4 功能基因注釋

通過(guò)Transdecoder version v 2.0.1對(duì)unigenes的CDS進(jìn)行了預(yù)測(cè)。利用Trinotate軟件包(v 2.0.0，https：//trinotate.github.io/)的標(biāo)準(zhǔn)流程對(duì)All-unigenes進(jìn)行了全面的注釋?zhuān)?Swissprot[17]、Pfam[18]、eggNOG[19]、Gene Ontology[20]、SignalP[21]、Rnammer[22]和KEGG Ontology數(shù)據(jù)庫(kù)注釋。

我們使用Trinotate軟件(版本：20140708)對(duì)預(yù)測(cè)出的編碼序列進(jìn)行功能注釋。Trinotate軟件整合了多個(gè)數(shù)據(jù)庫(kù)，可以從不同角度對(duì)轉(zhuǎn)錄本進(jìn)行功能注釋?zhuān)褂密浖驍?shù)據(jù)庫(kù)如下：使用diamond(版本：0.5.3)blastx，SwissProt、TrEMBL和kegg作為數(shù)據(jù)庫(kù)，對(duì)轉(zhuǎn)錄本進(jìn)行注釋?zhuān)煌瑫r(shí)，還使用blastp(版本：2.2.28+)，SwissProt作為數(shù)據(jù)庫(kù)，對(duì)轉(zhuǎn)錄本進(jìn)行注釋?zhuān)皇褂胔mmscan(HMMER)軟件(版本：3.1)，Pfam(Bateman et al.2002)作為數(shù)據(jù)庫(kù)，對(duì)預(yù)測(cè)的蛋白序列的結(jié)構(gòu)域進(jìn)行注釋?zhuān)挥肧ignalP(Petersen et al.2011)(版本：V4)軟件，預(yù)測(cè)信號(hào)肽；使用tmHMM(版本：2.0c)預(yù)測(cè)跨膜結(jié)構(gòu)；使用RNAmmer(版本：1.2)預(yù)測(cè)rRNA基因。

1.5 基因表達(dá)水平及表達(dá)差異量分析

基于Trinity拼接得到的轉(zhuǎn)錄本，我們進(jìn)行表達(dá)豐度計(jì)算。最先利用bowtie(版本：1.0.1；參數(shù)：mismatch=2)將每個(gè)樣品的reads分別比對(duì)到組裝出來(lái)的轉(zhuǎn)錄本上，得到bam格式的比對(duì)文件，然后使用RSEM(版本：v1.2.17)(http：//deweylab.biostat.wisc.edu/rsem/)計(jì)算比對(duì)到轉(zhuǎn)錄本上的reads的表達(dá)量。為了便于不同基因之間進(jìn)行表達(dá)量的比較，我們使用FPKM對(duì)表達(dá)量進(jìn)行標(biāo)準(zhǔn)化。FPKM(expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced)[23]是每百萬(wàn)fragments中來(lái)自某一基因每千個(gè)堿基長(zhǎng)度的fragments數(shù)目，其同時(shí)考慮了測(cè)序深度和基因長(zhǎng)度對(duì)fragments計(jì)數(shù)的影響。

差異表達(dá)基因(differentially expressed genes，DEGs)是基于基因表達(dá)水平分析中得到的read count數(shù)據(jù)，根據(jù)篩選閾值篩選出來(lái)的表達(dá)差異的基因。參考edgeR進(jìn)行差異基因分析，篩選閾值選擇為False Discovery Rate(FDR)≤0.05，|log2FoldChange|≥1。其中FC(Fold Change)表示葉與莖間表達(dá)量的比值，其絕對(duì)值越大表明差異倍數(shù)越小，F(xiàn)DR值越大則為表達(dá)差異越不顯著[24]。使用軟件GOseq(版本：3.0)對(duì)樣品間差異表達(dá)基因進(jìn)行GO庫(kù)注釋的分類(lèi)統(tǒng)計(jì)和功能富集分析，使用KOBAS(http：//ko-bas.cbi.pku.edu.cn/home.do)進(jìn)行KEGG pathway富集分析，使用Fisher精確檢驗(yàn)法對(duì)差異表達(dá)基因KEGG通路富集的顯著性進(jìn)行評(píng)價(jià)。使用超幾何檢驗(yàn)來(lái)找到差異表達(dá)基因中顯著富集的途徑。并經(jīng)過(guò)多次檢測(cè)和校正后，定義P <0.05的通路途徑為顯著富集差異表達(dá)基因的途徑。

1.6 SSR分析

以冬凌草轉(zhuǎn)錄組為研究對(duì)象，使用MISA(Microsatellite)軟件進(jìn)行SSR搜索。對(duì)Unigenes進(jìn)行SSR檢測(cè)，參數(shù)設(shè)置為“1～12、2～6、3～5、4～5、5～4、6～4”，1～12代表單堿基重復(fù)至少12次才算SSR，雙堿基6次，三堿基5次，以此類(lèi)推，重復(fù)單元最多有6個(gè)堿基。

2 結(jié)論與分析

2.1 總RNA質(zhì)量檢測(cè)和組裝結(jié)果分析

經(jīng)檢測(cè)，樣品總RNA的A260/A280均大于1.8，28 S/18 S在1.5～2.3之間，反應(yīng)RNA完整性的RIN值在8.5～9.5之間為質(zhì)量滿(mǎn)足建立和基因功能注釋、分類(lèi)及樣本間異同性分析。

通過(guò)組裝，葉與莖冬凌草樣品共獲得3 7961個(gè)Unigene序列。All-Unigene總長(zhǎng)為40366712 bp，平均長(zhǎng)度為1063 bp，反應(yīng)組裝效果的N50值為1810。其中N50是評(píng)估組裝質(zhì)量最常用的指標(biāo)之一，它反映了組裝結(jié)果的連續(xù)性和完整性。所得Unigene的數(shù)目隨長(zhǎng)度的增加而減少，其Unigene大多集中在在200～3000 bp之間，平均Unigene為1063 bp?？偟膩?lái)說(shuō)，Unigene的整體長(zhǎng)度分布均勻(圖1)。

圖1 All-Unigene長(zhǎng)度分布統(tǒng)計(jì)

2.2 Unigene蛋白數(shù)據(jù)庫(kù)比對(duì)結(jié)果

基于BLAST 算法對(duì)基因相似性分析，在37 961個(gè)Unigene 中，PFam庫(kù)中17842Unigene 同源匹配得到信息，在SignalP數(shù)據(jù)庫(kù)中有2394個(gè)，在eggNOG數(shù)據(jù)庫(kù)中有8769個(gè)，在 GO 數(shù)據(jù)庫(kù)中有18 776個(gè)，KEGG 數(shù)據(jù)庫(kù)中有23 936個(gè)(表1)。通常情況下測(cè)序的難度與基因組、轉(zhuǎn)錄組復(fù)雜程度有關(guān)。此外，由于缺乏參考基因組信息，當(dāng)前研究非模式生物的轉(zhuǎn)錄組仍具有一定挑戰(zhàn)。即使目前國(guó)內(nèi)對(duì)藥用植物的從頭轉(zhuǎn)錄組研究已有報(bào)道，但測(cè)序深度和基因檢出數(shù)目仍有待完善和提高[25-26]。

表1 All-Unigene蛋白數(shù)據(jù)庫(kù)統(tǒng)計(jì)結(jié)果

2.3 GO注釋和分類(lèi)

GO注釋和分類(lèi)是轉(zhuǎn)錄組數(shù)據(jù)研究中的一種國(guó)際標(biāo)準(zhǔn)化的基因功能分類(lèi)手段，它提供了一套動(dòng)態(tài)更新的標(biāo)準(zhǔn)詞匯表，用來(lái)全方面描述基因和基因產(chǎn)物的特性。GO功能分析由細(xì)胞組分、分子功能以及生物學(xué)過(guò)程三大部分組成。通過(guò)GO分析，37961Unigene中有18 775個(gè)獲得了GO注釋?zhuān)蛔⑨尀椤凹?xì)胞組分”有15 639個(gè)；16 455個(gè)被注釋為“分子功能”，1541515639被注釋為“生物學(xué)過(guò)程”。我們對(duì)其中的level2的GO注釋進(jìn)行了統(tǒng)計(jì)(圖2)。

2.4 KEGG 注釋和分類(lèi)

KEGG是全面分析生物體基因產(chǎn)物的代謝合成途徑以及產(chǎn)物功能的參考數(shù)據(jù)庫(kù)，根據(jù)生化途徑，提供基因的功能注釋和分類(lèi)信息，對(duì)比分析表達(dá)基因在生物學(xué)上的系統(tǒng)行為。將6456個(gè) Unigene 映射到 KEGG Pathway 數(shù)據(jù)庫(kù)中，注釋為代謝行為相關(guān)的序列共有4235個(gè)。其中涉及代謝(metabolism)、生物系統(tǒng)(organismal systems)環(huán)境信息處理(environmental information processing)、細(xì)胞進(jìn)程(cellular processes)及遺傳信息處理(genetic information processing)5個(gè)大類(lèi)及102條代謝通路，包含 Unigene 最多的是Ribosome(核糖體)共有389個(gè)，占糖酵解/糖異生(glycolysis/gluconeogenesis)共有 226個(gè)，占6%；其次是Oxidative phosphorylation(氧化磷酸化)共有244個(gè)，占3.8%，這些轉(zhuǎn)錄水平上的數(shù)據(jù)變化充分說(shuō)明了冬凌草藥用部位在特定發(fā)育時(shí)期的新陳代謝變化(圖3)。進(jìn)一步研究還發(fā)現(xiàn)，在代謝途徑中，與中藥有效活性成分生物合成相關(guān)Unigenes 序列分別是苯丙素類(lèi)合成(phenylpropanoid biosynthesis)101條、萜類(lèi)化合物骨架合成(terpenoid backbone biosynthesis)60 條、各種萜類(lèi)化合物合成(sesquiterpenoid，triterpenoid，monoterpenoid and diterpenoid biosynthesis)6條、(莨菪類(lèi)、哌啶類(lèi)、吡啶類(lèi))生物堿合成(tropane，piperidine andpyridine alkaloid biosynthesis)15條、黃酮類(lèi)成分合成(flavonoid，flavone，flavonol and isoflavonoid biosynthesis)30條(圖4)。

圖2 Level-2的GO注釋結(jié)果

圖3 KEGG 注釋和分類(lèi)結(jié)果

圖4 KEGG富集分析圖

2.5 二萜化合物合成途徑相關(guān)基因

植物二萜化合物主要合成途徑已經(jīng)明確，分為三步：最先，活性異戊烯焦磷酸(IPP)及其異構(gòu)體二甲基烯丙基焦磷酸(DMAPP)生成；之后為二萜化合物通用前體香葉基香葉基焦磷酸(GGPP)碳骨架的合成；最后在相關(guān)酶的作用下發(fā)生環(huán)化、羥基化等從而形成完整的二萜化合物分子[27]。利用本地化Blast的檢索方式在冬凌草中確定了本數(shù)據(jù)中與二萜化合物生物合成相關(guān)的候選基因。發(fā)現(xiàn)二萜合成途徑上的相關(guān)基因共有26個(gè)，其中與二萜骨架合成有關(guān)的基因有9個(gè)，其中與焦磷酸合成酶(copalyl diphosphate synthase)有關(guān)的 Unigene有6個(gè)，與貝殼杉烯合酶合成酶(KSL，kaurene synthase-like synthase)有關(guān)的unigene 有3個(gè)，具有細(xì)胞色素 P450(CYP450，cytochrome)功能的 Unigene共有5個(gè)。

2.6 差異表達(dá)基因的分析

通過(guò)冬凌草葉與莖兩個(gè)不同組織樣品的轉(zhuǎn)錄本/基因的表達(dá)量概率(FPKM scores)密度分布(圖5)，不同樣品間差異基因在KEGG 通路富集具有特別顯著性差異(P< 0.01)的通路描述，可以鑒定出一些在葉和莖兩種組織在次生物質(zhì)合成與代謝途徑中的差異表達(dá)基因。對(duì)冬凌草葉與莖兩種組織進(jìn)行差異表達(dá)分析，得到差異表達(dá)基因4565個(gè)，其中在莖中表達(dá)較高的為1668個(gè)，在葉中表達(dá)較高的有2697個(gè)，與二萜類(lèi)化合物合成相關(guān)的差異基因有15個(gè)。冬凌草葉與莖兩種組織差異基因表達(dá)圖如圖6。

圖5 冬凌草葉和莖基因表達(dá)量密度分布結(jié)果

2.7 SSR分析結(jié)果

簡(jiǎn)單重復(fù)序列又稱(chēng)微衛(wèi)星序列(Simple Sequence Repeat，SSR)由于具有高度多態(tài)性，克服一般分子標(biāo)記不確定、非特異的缺點(diǎn)，已在分子輔助育種、指紋數(shù)據(jù)庫(kù)的構(gòu)建、品種純度鑒定及遺傳多樣性分析等方面得到廣泛應(yīng)用[28-29]。在本次所測(cè)得的冬凌草轉(zhuǎn)錄組中，共檢測(cè)到8787條SSR序列，分布于7333條unigene中。SSR中二核苷酸重復(fù)的SSR標(biāo)記有4016條，占比45.70%，是主要類(lèi)型；其次是單核苷酸重復(fù)基序(2480條，占28.22%)，接著是三核苷酸重復(fù)基序(1707條，占比19.43%)。

圖6 冬凌草葉與莖差異表達(dá)基因的MA圖和火山圖

3 討論

本文通過(guò)Illumina Hiseq4000 平臺(tái)，利用RNA-seq 技術(shù)對(duì)道地產(chǎn)地冬凌草的葉與莖進(jìn)行了轉(zhuǎn)錄組測(cè)序與功能分析，初步建立了它們的基因表達(dá)圖譜，并通過(guò)比較藥用部位和非藥用部分合成途徑差異酶基因的發(fā)掘初步闡述了不同組織冬凌草成分差異的分子機(jī)制。并分析了轉(zhuǎn)錄組中的SSR，揭示了冬凌草中遺傳規(guī)律，也為藥用冬凌草的遺傳多樣性研究和分子標(biāo)記開(kāi)發(fā)奠定了分子基礎(chǔ)。

冬凌草葉和莖兩個(gè)組織樣品的高通量測(cè)序分析共獲得3 7961個(gè) Unigene序列。大量Unigene在多個(gè)數(shù)據(jù)庫(kù)中不能匹配到已知基因，可能由于數(shù)據(jù)庫(kù)中由于缺少冬凌草基因組信息；通過(guò)對(duì)轉(zhuǎn)錄組數(shù)據(jù)的差異表達(dá)分析發(fā)現(xiàn)與冬凌草合成途徑有關(guān)的差異基因主要注釋在苯丙素類(lèi)化合物、萜類(lèi)化合物骨架合成、生物堿類(lèi)化合物以及黃酮類(lèi)成分生物合成途徑等通路中。冬凌草中的活性成分為二萜類(lèi)物質(zhì)[6]，研究其生物合成途徑下游關(guān)鍵基因尤為重要?；诒疚霓D(zhuǎn)錄組數(shù)據(jù)，陳[30]等人成功克隆并初步驗(yàn)證了下游IrCYP71功能基因，利用轉(zhuǎn)錄組數(shù)據(jù)挖掘藥用植物的功能基因已經(jīng)成為重要手段[31]。對(duì)藥用植物的基因功能分化進(jìn)行深入分析，并與活性成分種類(lèi)、數(shù)量和積累變化相結(jié)合，進(jìn)而從分子生藥學(xué)角度揭示由于不同生長(zhǎng)時(shí)間和生長(zhǎng)環(huán)境基因表達(dá)的差異與藥用部位中次生代謝成分差異的機(jī)制，才能為藥材質(zhì)量控制與提高提供理論依據(jù)[32]。高通量轉(zhuǎn)錄組測(cè)序技術(shù)的應(yīng)用同時(shí)也加速了SSR 位點(diǎn)的挖掘，本研究中冬凌草的 All-Unigenes 轉(zhuǎn)錄本中 SSR 出現(xiàn)的頻率特點(diǎn)為二核苷酸的重復(fù)單元類(lèi)型最為豐富，這為后續(xù)冬凌草的種質(zhì)資源和遺傳多樣性研究提供了豐富的分子標(biāo)記數(shù)據(jù)。通過(guò)構(gòu)建藥用植物生長(zhǎng)發(fā)育的不同時(shí)期、不同器官和不同產(chǎn)地的轉(zhuǎn)錄組文庫(kù)，可以全方位地研究藥用植物中活性成份生物合成的關(guān)鍵酶基因的表達(dá)水平，從而有效地揭示次生代謝產(chǎn)物生物合成的途徑及其調(diào)控機(jī)制。