藥用資源植物山莨菪的轉(zhuǎn)錄組信息分析

張 雨 夏銘澤 張發(fā)起

(1.中國(guó)科學(xué)院高原生物適應(yīng)與進(jìn)化重點(diǎn)實(shí)驗(yàn)室，中國(guó)科學(xué)院西北高原生物研究所，西寧 810001； 2.中國(guó)科學(xué)院大學(xué)，北京 100039)

山莨菪(Anisodustanguticus)隸屬茄科(Solanaceae)山莨菪屬(Anisodus)，是多年生宿根草本植物，產(chǎn)于青海、甘肅、西藏東部、云南西北部等地區(qū)，海拔在2 800～4 200 m，生于草坡陽面。山莨菪是重要的資源植物，其根可供藥用；莨菪烷類生物堿可以從山莨菪中提取，在醫(yī)學(xué)上具有麻醉、鎮(zhèn)痛和解毒等功效；其地上部分可以摻入牛飼料中，有催膘作用[1～2]。山莨菪具有重要的藥用價(jià)值和經(jīng)濟(jì)價(jià)值，是一種重要的資源植物。但是由于山莨菪分布生境狹窄，種植量有限，需求量又較高，故被大量挖掘，其資源的可持續(xù)利用遭到了嚴(yán)重威脅，在我國(guó)《國(guó)家重點(diǎn)保護(hù)野生植物名錄》中已被列為國(guó)家Ⅱ級(jí)保護(hù)植物。目前，對(duì)山莨菪的研究主要集中在對(duì)山莨菪堿的藥理學(xué)、病理學(xué)和臨床應(yīng)用研究方面[3～4]，部分研究集中在山莨菪的光合作用、生長(zhǎng)特征、元素含量及傳粉生物學(xué)方面[5～8]，只有少數(shù)人對(duì)山莨菪進(jìn)行分子生物學(xué)方面的研究，分子標(biāo)記開發(fā)較為落后，且轉(zhuǎn)錄組和基因組信息相對(duì)缺乏[9～10]。

轉(zhuǎn)錄組是組織或細(xì)胞轉(zhuǎn)錄產(chǎn)物的總和，通過轉(zhuǎn)錄組分析，可以進(jìn)一步對(duì)基因的結(jié)構(gòu)和功能進(jìn)行研究，便于了解基因組的組成和功能，深入探索生物體表達(dá)狀況、發(fā)育的機(jī)理及調(diào)控網(wǎng)絡(luò)機(jī)制等[11～13]。近年來，高通量測(cè)序技術(shù)的迅速發(fā)展降低了測(cè)序成本，有利于轉(zhuǎn)錄組分析的發(fā)展[14]，目前已有多種資源植物完成了轉(zhuǎn)錄組測(cè)序，如李東[15]對(duì)熱脅迫下丹參(Salviamiltiorrhiza)的次生代謝比較轉(zhuǎn)錄組進(jìn)行分析，發(fā)現(xiàn)了熱脅迫對(duì)丹參轉(zhuǎn)錄組、對(duì)迷迭香酸途徑關(guān)鍵酶表達(dá)及萜類合成途徑關(guān)鍵酶表達(dá)的影響，為解釋熱脅迫與有效成分合成與積累過程提供了依據(jù)；李依民等[16]通過高通量測(cè)序技術(shù)對(duì)黃三七(Soulieavaginata)的根莖轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行分析，獲得了豐富的轉(zhuǎn)錄組數(shù)據(jù)，為黃三七代謝途徑解析、調(diào)控機(jī)制研究及基因功能鑒定奠定了基礎(chǔ)；李太強(qiáng)等[17]對(duì)長(zhǎng)梗杜鵑(Rhododendroncalvescens)進(jìn)行轉(zhuǎn)錄組分析，發(fā)現(xiàn)了176條和人類疾病相關(guān)的Unigenes，其中與內(nèi)分泌及代謝疾病相關(guān)的Unigenes有167條，與抗藥性有關(guān)的Unigenes有9條，為杜鵑屬的抗性機(jī)理研究及相關(guān)病理學(xué)研究提供了重要參考；另外還有許多研究人員[18～20]對(duì)黑果枸杞(Lyciumruthenicum)、余干子(Phyllanthusemblica)、陸地棉(Gossypiumhirsutum)等資源植物進(jìn)行了轉(zhuǎn)錄組分析，獲得了大量數(shù)據(jù)信息。

鑒于此，本研究利用高通量測(cè)序技術(shù)，對(duì)山莨菪進(jìn)行轉(zhuǎn)錄組測(cè)序分析，得到大量測(cè)序數(shù)據(jù)，結(jié)合生物信息學(xué)方法，拼接、組裝和注釋測(cè)序序列，對(duì)山莨菪的基因結(jié)構(gòu)、表達(dá)水平和差異、基因富集進(jìn)行分析，以期為山莨菪轉(zhuǎn)錄組和基因組水平的研究提供數(shù)據(jù)，同時(shí)也對(duì)山莨菪的保護(hù)及合理開發(fā)利用奠定分子生物學(xué)基礎(chǔ)。

1 材料與方法

1.1 材料

本研究的植物材料山莨菪采于青海省同仁縣(地理坐標(biāo)N35°13′58.30″，E101°51′05.50″；海拔3 532 m)，采集一年生幼葉，用純水清洗去除污漬、75%酒精殺菌消毒后迅速置于液氮中處理，再轉(zhuǎn)移到-80℃的超低溫冰箱保存?zhèn)溆?。憑證標(biāo)本(標(biāo)本號(hào)：Zhang2018050)存于中國(guó)科學(xué)院西北高原生物研究所青藏高原生物標(biāo)本館(HNWP)。

1.2 方法

1.2.1 山莨菪RNA提取與測(cè)序

提取山莨菪的RNA，在對(duì)其RNA進(jìn)行純化和質(zhì)量分析之后，構(gòu)建cDNA文庫(kù)，采用Illumina HiseqTM測(cè)序平臺(tái)進(jìn)行測(cè)序。

1.2.2 序列的拼接與組裝

通過FastQC進(jìn)行數(shù)據(jù)評(píng)估和數(shù)據(jù)質(zhì)控，通過Trimmomatic進(jìn)行質(zhì)量剪切，在質(zhì)量剪切的過程中，需要去除以下情況的序列：①N堿基序列；②reads中的接頭序列及Q值低于20的低質(zhì)量堿基；③長(zhǎng)度不高于35nt的reads及配對(duì)序列；④被污染的序列等。最終，得到相對(duì)準(zhǔn)確的有效數(shù)據(jù)，使用Trinity將有效數(shù)據(jù)從頭組裝成轉(zhuǎn)錄本，大致過程如下：

①運(yùn)行seqtk-trinity程序，將fq數(shù)據(jù)轉(zhuǎn)換成fa數(shù)據(jù)，并生成both.fa；

②運(yùn)行Jellyfish程序，生成K-mer目錄；

③運(yùn)行Inchworm程序，組裝reads初步得到contigs；

④運(yùn)行Chrysalis程序，將上一步生成的contigs進(jìn)行聚類分析，對(duì)每個(gè)類構(gòu)建Bruijn圖；

⑤運(yùn)行Butterfly程序，處理上一步得到的Bruijn圖，根據(jù)圖中reads和成對(duì)的reads尋找途徑，得到最終的isoform序列；

⑥得到Trinity.fasta文件，其haeder信息包含了轉(zhuǎn)錄本名稱、長(zhǎng)度、De Brujin圖重構(gòu)途徑等；

⑦對(duì)拼接組裝得到的轉(zhuǎn)錄本去冗余。

在此基礎(chǔ)上，統(tǒng)計(jì)所得轉(zhuǎn)錄本的各項(xiàng)信息，如GC含量、長(zhǎng)度分布狀況、所含isoform數(shù)目分布狀況等，為后續(xù)分析做準(zhǔn)備。

1.2.3 基因功能注釋

通過Blast將獲得的Unigenes與多個(gè)數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，取evalue小于等于1e-10，相似度<90%且覆蓋度<80%的比對(duì)結(jié)果，再對(duì)山莨菪進(jìn)行基因功能的注釋及分類分析?？蛇M(jìn)行對(duì)比參考的數(shù)據(jù)庫(kù)有NT(核酸序列數(shù)據(jù)庫(kù))、NR(非冗余蛋白質(zhì)序列數(shù)據(jù)庫(kù))、COG(原核生物基因直系同源關(guān)系注釋系統(tǒng))、KOG(真核生物基因直系同源關(guān)系注釋系統(tǒng))、Swiss～Prot(最常用且全面的蛋白質(zhì)數(shù)據(jù)庫(kù))、TrEMBL(蛋白質(zhì)數(shù)據(jù)庫(kù)，Swiss～Prot的增補(bǔ)本)、PFAM(最全面蛋白結(jié)構(gòu)域注釋分類系統(tǒng))、CDD(保守區(qū)域結(jié)構(gòu)數(shù)據(jù)庫(kù))、GO(國(guó)際化基因功能分類體系)、KEGG(有關(guān)生物系統(tǒng)的較完善數(shù)據(jù)庫(kù))等。通過與NR、Swiss～Prot、TrEMBL數(shù)據(jù)庫(kù)進(jìn)行對(duì)比，可以對(duì)CDS序列進(jìn)行預(yù)測(cè)；山莨菪的轉(zhuǎn)錄本序列與相近物種的近似情況，以及山莨菪同源序列的功能信息可以通過與NR庫(kù)對(duì)比得到；對(duì)山莨菪的基因進(jìn)行GO分類后，可以便于統(tǒng)計(jì)基因在生物過程、細(xì)胞成分、分子功能下的GO term；另外，對(duì)基因做KO注釋之后，可以根據(jù)KO與pathway的聯(lián)系對(duì)山莨菪進(jìn)行KEGG代謝通路的分類分析[21～23]。

2 結(jié)果與分析

2.1 轉(zhuǎn)錄本組裝與CDS預(yù)測(cè)分析

從頭測(cè)序組裝之后得到158 378個(gè)Transcripts，得到71 463個(gè)Unigenes，對(duì)編碼序列進(jìn)行預(yù)測(cè)之后得到47 685個(gè)CDS(見表1)。對(duì)Transcripts而言(圖1A)，分布在200～300 bp的數(shù)量是最多的，其次是300～400 bp(42 231條)、大于等于2 000 bp(19 215條)和300～400 bp(14 864條)。對(duì)Unigenes而言(圖1B)，按分布范圍內(nèi)數(shù)量排序的話，依次是200～300 bp(29 326條)、300～400 bp(11 003條)、400～500 bp(5 955條)、大于等于2 000 bp(4 146條)；對(duì)CDS而言，長(zhǎng)度在100～200 nt內(nèi)的序列最多，有12 928條，占27.11%，其次為200～300 nt(11 259條)、300～400 nt(4 838條)、400～500 nt(2 980條)。無論是Transcripts、Unigenes還是CDS，長(zhǎng)度最短的序列數(shù)量最多，而隨著序列長(zhǎng)度增加，所獲得的拼接數(shù)量就越少(大于等于2 000 bp的序列除外)。Unigenes的平均長(zhǎng)度為651.1 bp，最短的Unigene長(zhǎng)度為201 bp，最長(zhǎng)為8 526 bp，N50長(zhǎng)度為1 115 bp，總長(zhǎng)度為46 529 443 bp。

圖1 山莨菪Transcript與Unigenes的長(zhǎng)度分布圖Fig.1 Distribution of transcript and unigene length for A.tanguticus

表1 拼接結(jié)果統(tǒng)計(jì)

2.2 功能注釋及COG分類

以CDD、KOG、NR、NT、GO、Swissprot、PFAM、KEGG、TrEMBL這9個(gè)數(shù)據(jù)庫(kù)作為參考，將de novo組裝得到的71 463個(gè)基因序列與之進(jìn)行比對(duì)和注釋(見表2)。注釋到CDD、KOG、NR、NT、PFAM、Swissprot、TrEMBL、GO、KEGG數(shù)據(jù)庫(kù)種的Unigenes分別有19 526、13 110、39 621、38 769、14 393、24 163、39 653、29 309和3 679條，注釋到CDD數(shù)據(jù)庫(kù)種的Unigenes最多，占比27.32%，而注釋到KEGG數(shù)據(jù)庫(kù)的Unigenes最少，占5.15%。在9個(gè)數(shù)據(jù)庫(kù)中至少1個(gè)數(shù)據(jù)庫(kù)注釋成功的Unigenes為47 624條，占總Unigenes數(shù)的66.64%；共有2 415個(gè)Unigenes與9個(gè)數(shù)據(jù)庫(kù)都能匹配成功，占總Unigenes數(shù)的3.38%(見表2)；另外，尚未注釋成功的Unigenes數(shù)量較多，有23 839條，約占1/3。

表2 Unigenes注釋結(jié)果

圖2 山莨菪轉(zhuǎn)錄組Unigenes的KOG功能分布圖 A.RNA加工和修飾RNA；B.染色體結(jié)構(gòu)和動(dòng)力學(xué)；C.能源生產(chǎn)與轉(zhuǎn)化；D.細(xì)胞周期調(diào)控、細(xì)胞分裂、染色體分離；E.氨基酸轉(zhuǎn)運(yùn)和代謝；F.核酸轉(zhuǎn)運(yùn)和代謝；G.碳水化合物轉(zhuǎn)運(yùn)和代謝；H.輔酶轉(zhuǎn)運(yùn)和代謝；I.脂類轉(zhuǎn)運(yùn)和代謝；J.翻譯、核糖體結(jié)構(gòu)和生物發(fā)生；K.轉(zhuǎn)錄；L.復(fù)制、重組和修飾；M.細(xì)胞壁/細(xì)胞膜生物發(fā)生；N.細(xì)胞活性；O.翻譯后修飾、蛋白轉(zhuǎn)運(yùn)；P.無機(jī)離子轉(zhuǎn)運(yùn)和代謝；Q.次生代謝物合成、轉(zhuǎn)運(yùn)和代謝；R.只有一般功能預(yù)測(cè)；S.未知功能；T.信號(hào)傳遞機(jī)制；U.細(xì)胞間運(yùn)輸、分泌物和囊泡運(yùn)動(dòng)；V.防御機(jī)制；W.細(xì)胞外結(jié)構(gòu)；Y.核結(jié)構(gòu)；Z.細(xì)胞骨架Fig.2 KOG functional annotation distribution of unigenes of transcriptome for A.tanguticus A.RNA processing and modification; B.Chromatin structure and dynamic; C.Energy production and conversion; D.Cell cycle control,cell division,chromosome partitioning; E.Amina acid transport and metabolism; F.Nucleotide transport and metabolism; G.Carbohydrate transport and metabolism; H.Coenzyme transport and metabolism; I.Lipid transport and metabolism; J.Translation,ribosomal structure and biogenesis; K.Transcription; L.Replication,recombination and repair; M.Cell wall/membrane/envelope biogenesis; N.Cell motility; O.Posttranslational modification,protein turnover,chaperones; P.Inorganic ion transport and meabolism; Q.Secodary metabolites biosynthesis,transport and catabolism; R.General function prediction only; S.Function unknown; T.Signal transduction mechanisms; U.Intracelluar trafficking,secretion,and vesicular transport; V.Defense mechanisms; W.Extracellular structures; Y.Nuclear structure; Z.Cytoskeleton

2.3 KOG功能注釋

將Unigenes與KOG蛋白質(zhì)庫(kù)進(jìn)行比對(duì)，有13 110條Unigenes獲得注釋，占總Unigenes的18.35%，Unigenes根據(jù)其功能被分為了26個(gè)類別，都與山莨菪最基本的生命活動(dòng)相關(guān)(見圖2)。其中，在已經(jīng)明確生物學(xué)功能的Unigenes中(占94.65%)，信號(hào)傳遞機(jī)制最多，有1 873條Unigenes，占比為12.76%，其次是只有一般功能預(yù)測(cè)和翻譯后修飾，分別為1 703(11.60%)和1 580(10.76%)。與轉(zhuǎn)錄、翻譯、化合物運(yùn)輸、能量生產(chǎn)轉(zhuǎn)化、代謝等功能相關(guān)的Unigenes處于相對(duì)較多的數(shù)量，數(shù)量在500～900條，而注釋出來與染色體結(jié)構(gòu)、核結(jié)構(gòu)、細(xì)胞骨架、細(xì)胞周期等相關(guān)序列是較少的，有的甚至只有不到10條。

2.4 NR功能注釋

通過與NR庫(kù)(evalue<0.000 01)進(jìn)行比對(duì)，得到39 621個(gè)被注釋的Unigenes，未獲得注釋的基因較少，僅有22個(gè)(0.056%)。與山莨菪能比對(duì)上有8個(gè)同科(茄科)植物，這8個(gè)同科植物主要是茄屬、辣椒屬和煙草屬的植物，彼此之間具有較近的親緣關(guān)系，另外還有兩種非茄科的植物，分別是豆科和藜科植物，所占比例較小，有3 422(8.64%)的Unigenes零星分布于其他392個(gè)物種中。在所有的注釋中，注釋為與陽芋(Solanumtuberosum)相關(guān)的Unigenes最多，有8 598條(21.70%)，其次為辣椒(Capsicumannuum)和煙草(Nicotianatabacum)，分別占16.24%和13.45%，比例所占較大的前幾種同科植物均可說明與山莨菪具有較高的序列同源性。

2.5 GO功能注釋及分類

為全面描述山莨菪基因和基因產(chǎn)物的屬性，了解山莨菪表達(dá)基因的功能分布狀況及其所代表的生物學(xué)意義，將Unigenes通過與Swissprot、TrEMBL數(shù)據(jù)庫(kù)的比對(duì)，得到GO功能注釋信息，注釋后得到29 309個(gè)Unigenes，對(duì)其進(jìn)行統(tǒng)計(jì)分析后(表3)，可分為分子功能、生物學(xué)過程和細(xì)胞組分3個(gè)大類，分別有16、24、22個(gè)子類，共有62個(gè)子類，其中結(jié)合、催化活性、細(xì)胞過程、代謝過程、細(xì)胞、細(xì)胞部分及細(xì)胞器等獲得注釋較多，而注釋為金屬伴活動(dòng)分子功能調(diào)節(jié)器、結(jié)構(gòu)分子活性、傳譯調(diào)治活性、細(xì)胞、膜封閉腔、擬核的Unigenes基本是最少的。

圖3 山莨菪轉(zhuǎn)錄組Unigenes的KEGG功能注釋分布統(tǒng)計(jì)圖 A.細(xì)胞生長(zhǎng)和死亡；B.細(xì)胞運(yùn)動(dòng)；C.細(xì)胞群體；D.運(yùn)輸和分解代謝；E.膜運(yùn)輸；F.信號(hào)轉(zhuǎn)導(dǎo)；G.信號(hào)分子和互作作用；H.折疊排序與退化；I.復(fù)制和修復(fù)；J.轉(zhuǎn)錄；K.翻譯；L.氨基酸代謝；M.其他次級(jí)代謝產(chǎn)物合成；N.碳水化合物代謝；O.能量代謝；P.多糖合成與代謝；Q.類脂(化合)物代謝作用；R.輔助因子和維生素的代謝；S.其他氨基酸的代謝；T.萜類化合物和聚酮化合物的代謝；U.Nucleotide metabolism；V.概觀；W.異生素生物降解和新陳代謝Fig.3 KEGG functional annotation distribution of unigenes of transcriptome for A.tanguticus A.Cell growth and death; B.Cell motility; C.Cellular community; D.Transcript and catabolism; E.Membrane transport; F.Signal transduction; G.Signaling molecules and interaction; H.Folding sorting and degradation; I.Replication and repair; J.Transcription; K.Translation; L.Amino acid metabolism; M.Biosynthesis of other secondary metabolites; N.Carbohydrate metaboliam; O.Energy metabolism; P.Glycan biosynthesis and metabolism; Q.Lipid metabolism; R.Metabolism of cofactros and vitamins; S.Metabolism of other amino acids; T.Metabolism of terpenoids and polyketides; U.Nucleotied metabolism; V.Overview; W.Xenobiotics biodegradation and metabolism

2.6 KEGG代謝通路分析

以KEGG代謝庫(kù)(evalue<0.000 01)為參考，對(duì)代謝通路進(jìn)行統(tǒng)計(jì)和分類分析。有3 679條Unigenes被注釋，這3 679條被注釋的基因根據(jù)涉及到的代謝通路分為4大類，23個(gè)子類(見圖3)。在4種代謝大類中，代謝相關(guān)的通路獲得注釋最多，為2 640個(gè)，占比49.45%，其次分別是遺傳信息處理、細(xì)胞過程、環(huán)境信息處理，分別獲得了1 438(26.93%)、656(12.29%)、605(11.33%)個(gè)注釋。23個(gè)亞類中，與代謝相關(guān)的通路最多，有12條，包括氨基酸代謝、碳水化合物代謝、其他次生物質(zhì)代謝、能量代謝、糖生物合成與代謝、脂類物質(zhì)代謝、輔助因子和維生素代謝、其他氨基酸代謝、萜類化合物和聚酮化合物代謝、核苷酸代謝、輔助因子代謝、概觀；與遺傳信息處理、細(xì)胞過程、環(huán)境信息處理相關(guān)的通路較少，分別有4條、4條、3條。在這23個(gè)子類代謝途徑中，翻譯獲得的注釋最多，為595條，占11.14%，其次為信號(hào)轉(zhuǎn)導(dǎo)、碳水化合物代謝，折疊排序和退化，分別為576條(10.79%)、479條(8.97%)和428條(8.01%)。

表3 山莨菪轉(zhuǎn)錄組Unigenes的GO功能分類統(tǒng)計(jì)

有7 075條Unigenes歸入到209條代謝途徑中，按照基因獲得的注釋數(shù)量從高到低排列，將前13個(gè)代謝通路列于表4中，其中以核糖體代謝通路最多，占3.92%，其次為碳代謝和植物激素信號(hào)轉(zhuǎn)導(dǎo)，分別占2.64%和2.45%。

2.7 山莨菪藥用活性成分分析

山莨菪植株體內(nèi)含有生物堿、萜類、苯丙素類、黃酮類、糖苷類、醌類、聚酮類、有機(jī)酸及酚類等藥效成分(表5)，其中生物堿、萜類、黃酮類和糖苷類成分的代謝通路分別有5條(以異喹啉生物堿生物合成途徑，二苯乙烯類、二芳基庚烷類和姜辣素生物合成途徑為主)、4條(以萜類骨架生物合成途徑為主)、3條(以類黃酮生物合成為主)、2條(以泛醌和其他萜類化合物—醌生物合成途徑為主)，其余藥效成分目前基本均只有1條代謝途徑。在所有藥效成分中，萜類和苯丙素類所對(duì)應(yīng)的unigenes是最多的，均有67條，生物堿、醌類及黃酮類所對(duì)應(yīng)的unigenes數(shù)量也較多，分別有30、28、22條，而糖苷類、聚酮類、有機(jī)酸及酚類對(duì)應(yīng)的unigenes是相對(duì)較少的，分別只有2，2和11條。對(duì)山莨菪不同藥效成分的代謝通路、所對(duì)應(yīng)unigenes的數(shù)量及類型進(jìn)行分析，將便于復(fù)雜基因調(diào)控網(wǎng)絡(luò)的建立，有助于為人工合成所需藥用成分提供研究基礎(chǔ)和思路。

表4 山莨菪Unigenes數(shù)量最多的13個(gè)代謝通路

Table 4 Top thirteen metabolic pathways involved inA.tanguticusunigenes

代謝通路Pathway ID代謝通路PathwayUnigenes數(shù)量及比例Number and propotion of unigenesko03010核糖體Ribosome277(3.92%)ko01200碳代謝Carbon metabolism187(2.64%)ko04075植物激素信號(hào)轉(zhuǎn)導(dǎo)Plant hormone signal transduction173(2.45%)ko01230氨基酸合成Biosynthesis of amino acids167(2.36%)ko03040剪接Spliceosome158(2.23%)ko04141內(nèi)質(zhì)網(wǎng)蛋白加工Protein processing in endoplasmic reticulum156(2.20%)ko00230嘌呤代謝Purine metabolism137(1.94%)ko03013RNA轉(zhuǎn)運(yùn) RNA transport137(1.94%)ko00500淀粉和蔗糖代謝Starch and sucrose metabolism128(1.81%)ko00190氧化磷酸化Oxidative phosphorylation127(1.80%)ko04144胞吞作用Endocytosis113(1.60%)ko00240嘧啶代謝Pyrimidine metabolism104(1.47%)

2.8 SSR與SNP檢測(cè)分析

使用BCFtools根據(jù)Mapping結(jié)果找出可能的SNP位點(diǎn)(單核苷酸位點(diǎn)多態(tài)性)，最終一共獲得31 382個(gè)SNP位點(diǎn)，3 675個(gè)InDel，分析發(fā)現(xiàn)山莨菪Unigenes序列上SNP分布不均勻，轉(zhuǎn)換突變類型的SNP數(shù)量比顛換突變的2倍略多(見圖4)。其中轉(zhuǎn)換突變類型(A→G、C→T、G→A、T→C)有20 168個(gè)，占64.26%，顛換突變類型(A→C、A→T、C→A、C→G、G→C、G→T、T→A、T→G)有11 219個(gè)，占35.74%。采用MISA基于拼接所得轉(zhuǎn)錄本序列信息進(jìn)行SSR(簡(jiǎn)單序列重復(fù)標(biāo)記)分析，圖5的SSR密度分布圖顯示，有6種SSR重復(fù)類型，其中單堿基重復(fù)、二堿基重復(fù)和三堿基重復(fù)類型占絕大多數(shù)，每百萬堿基中出現(xiàn)的單堿基重復(fù)的SSR個(gè)數(shù)有56.52個(gè)，占45.30%，二堿基重復(fù)和三堿基重復(fù)的個(gè)數(shù)分別為28.22和37.29，分別占22.62%、29.89%，而四堿基重復(fù)、五堿基重復(fù)和六堿基重復(fù)所占比例僅占2.19%。豐富的SSR和SNP位點(diǎn)對(duì)山莨菪遺傳圖譜構(gòu)建、遺傳多樣性分析和親緣關(guān)系研究具有重要的意義和價(jià)值。

表5 山莨菪藥用活性成分代謝通路及基因統(tǒng)計(jì)

圖4 山莨菪轉(zhuǎn)錄組Unigenes的突變譜系圖Fig.4 Mutation pedigree of unigenes of transcriptome for A.tanguticus

圖5 山莨菪轉(zhuǎn)錄組Unigenes的SSR密度分布圖Fig.5 SSR Density Distribution Map of unigenes of transcriptome for A.tanguticus

3 討論

本研究采用Illumina測(cè)序平臺(tái)對(duì)山莨菪進(jìn)行高通量測(cè)序、從頭組裝、拼接、去冗余，結(jié)合多個(gè)數(shù)據(jù)庫(kù)對(duì)基因進(jìn)行注釋和分析，一共獲得了158 378條Transcripts，71 463條Unigenes，對(duì)編碼序列進(jìn)行預(yù)測(cè)后獲得，47 685條CDS。最長(zhǎng)的Unigenes為8 526 bp，最短的為201 bp，N50為1 115 bp，平均長(zhǎng)度是651.1 bp，結(jié)合一些分析較為成熟的常見資源植物的拼接結(jié)果，如芝麻Sesamumindicum(N50為1 006 bp，平均長(zhǎng)度683 bp)[24]、樟樹Cinnamomumbodinieri(N50為1 023 bp，平均長(zhǎng)度584 bp)[25]、百合Liliumbrownii(平均長(zhǎng)度673 bp)[26]、鐵皮石斛Dendrobiumofficinale(平均長(zhǎng)度660 bp)[27]，可知山莨菪的序列拼接質(zhì)量與上述植物相近，而仿刺參Apostichopusjaponicus與山莨菪雖然在基因組和轉(zhuǎn)錄組發(fā)展方面具有類似的情況，但其拼接后序列的平均長(zhǎng)度卻僅為344 bp[28]，因此，在山莨菪分子生物學(xué)方面的研究較少、轉(zhuǎn)錄組和基因組信息相對(duì)缺乏的情況下，還能獲得與轉(zhuǎn)錄組分析較為成熟的植物比較相近的拼接結(jié)果，可知其拼接質(zhì)量是較高的。

將Unigenes與9大數(shù)據(jù)庫(kù)(CDD、KOG、NR、NT、GO、Swissprot、PFAM、KEGG、TrEMBL)進(jìn)行對(duì)比，注釋成功的Unigenes有47 624(66.64%)條，但仍有23 839(33.36%)條序列未能獲得注釋，數(shù)量較多，究其原因，可能有以下3點(diǎn)：①一般而言，在測(cè)序質(zhì)量較高的前提下，序列片段的長(zhǎng)度越長(zhǎng)，注釋的成功率和可靠性就越高，因而山莨菪的拼接質(zhì)量較優(yōu)，但結(jié)果卻有33.36%的Unigenes無法獲得注釋，這可能是因?yàn)槠湫蛄斜旧黹L(zhǎng)度就較短，因此很難在同源性比對(duì)時(shí)獲得結(jié)果；②目前山莨菪的分子生物學(xué)研究仍處于初級(jí)階段，基因組學(xué)研究的匱乏導(dǎo)致基因功能注釋信息的缺失，因此會(huì)有部分序列很難得到注釋信息；③山莨菪本身轉(zhuǎn)錄組可能具有一定的特異性，由于之前尚未有山莨菪轉(zhuǎn)錄組方面的研究作為參考，相關(guān)功能基因的研究也較少，所以可能使得部分特異性基因得不到識(shí)別和注釋。因此，測(cè)序技術(shù)及質(zhì)量的提高，基礎(chǔ)功能數(shù)據(jù)庫(kù)的不斷完善，以及對(duì)山莨菪基因組學(xué)方面研究的進(jìn)一步探索，很大程度上可有效減少未被注釋序列的數(shù)量。

成功注釋到KOG數(shù)據(jù)庫(kù)中的Unigenes有131 100條，根據(jù)表1中Unigenes的注釋數(shù)量及比例可知，注釋到該庫(kù)的比例與注釋到NR數(shù)據(jù)庫(kù)的55.44%、注釋到NT數(shù)據(jù)庫(kù)的54.25%、注釋到GO的41.01%以及注釋到TrEMBL的55.49%相比而言，注釋的比例相對(duì)較低，僅有18.35%，排除這幾個(gè)數(shù)據(jù)庫(kù)本身注釋難易程度差異的影響，這可能還與數(shù)據(jù)庫(kù)內(nèi)部數(shù)據(jù)的豐富程度、山莨菪本身轉(zhuǎn)錄組信息是否完善等因素有關(guān)。從KOG功能注釋圖中可以了解到，在獲得注釋的26類Unigenes中，信號(hào)傳遞、轉(zhuǎn)錄翻譯、能量生產(chǎn)轉(zhuǎn)化、代謝等相關(guān)基因被注釋的數(shù)量較多，而與細(xì)胞結(jié)構(gòu)、周期相關(guān)的基因數(shù)量卻較少，由被注釋功能基因的數(shù)量可以推測(cè)不同基因表達(dá)豐富度的差異，從而判斷不同生命活動(dòng)對(duì)山莨菪生長(zhǎng)發(fā)育的重要程度。

在與NR數(shù)據(jù)庫(kù)比對(duì)之后，注釋為陽芋同源種的Unigenes有8 598條，遠(yuǎn)高于其他幾個(gè)被注釋物種，出現(xiàn)這種情況，可能是由于山莨菪與陽芋的進(jìn)化史和生活史較為接近，可以說明山莨菪與陽芋的親緣關(guān)系更緊密一點(diǎn)，而與其他物種的親緣關(guān)系相對(duì)較遠(yuǎn)。對(duì)拼接組裝獲得的Unigenes進(jìn)行了代謝途徑分析和功能分類，共獲得3 679條Unigenes注釋，注釋可細(xì)分為4個(gè)代謝通路大類，其中定位到代謝通路相關(guān)的基因最多，占49.45%，證明山莨菪具有較強(qiáng)的代謝活動(dòng)能力。進(jìn)一步對(duì)山莨菪的藥用活性成分(物堿、萜類、苯丙素類、黃酮類、糖苷類、醌類、聚酮類、有機(jī)酸及酚類等)的代謝通路及相關(guān)Unigenes數(shù)量和類型進(jìn)行統(tǒng)計(jì)分析，發(fā)現(xiàn)與生物堿相關(guān)的代謝通路最多，而萜類和苯丙素類所對(duì)應(yīng)的Unigenes數(shù)量最多，該分析將有助于推進(jìn)山莨菪藥用成分合成機(jī)理的探索，為所需成分的提取和復(fù)合成分的合成提供研究基礎(chǔ)。

轉(zhuǎn)錄組測(cè)序分析后，共檢測(cè)出31 382個(gè)SNP多態(tài)位點(diǎn)，兩種突變類型中以轉(zhuǎn)換突變?yōu)橹鳎?/3；在SSR分析中，獲得六種SSR重復(fù)類型，其中以單堿基重復(fù)、二堿基重復(fù)和三堿基重復(fù)為主。通過SNP分析與SSR分析，可以構(gòu)建較高分辨率的山莨菪遺傳圖譜，有助于提高目的基因定位的準(zhǔn)確性，同時(shí)促進(jìn)山莨菪群體遺傳學(xué)、比較基因組學(xué)和分子系統(tǒng)學(xué)的發(fā)展。