張 雨 夏銘澤 張發(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áng)面。山莨菪是重要的資源植物，其根可供藥用；莨菪烷類生物堿可以從山莨菪中提取，在醫(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)記開(kāi)發(fā)較為落后，且轉(zhuǎn)錄組和基因組信息相對(duì)缺乏[9～10]。

轉(zhuǎn)錄組是組織或細(xì)胞轉(zhuǎn)錄產(chǎn)物的總和，通過(guò)轉(zhuǎn)錄組分析，可以進(jìn)一步對(duì)基因的結(jié)構(gòu)和功能進(jìn)行研究，便于了解基因組的組成和功能，深入探索生物體表達(dá)狀況、發(fā)育的機(jī)理及調(diào)控網(wǎng)絡(luò)機(jī)制等[11～13]。近年來(lái)，高通量測(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á)的影響，為解釋熱脅迫與有效成分合成與積累過(guò)程提供了依據(jù)；李依民等[16]通過(guò)高通量測(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ù)及合理開(kāi)發(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 序列的拼接與組裝

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

①運(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 基因功能注釋

通過(guò)Blast將獲得的Unigenes與多個(gè)數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，取evalue小于等于1e-10，相似度<90%且覆蓋度<80%的比對(duì)結(jié)果，再對(duì)山莨菪進(jìn)行基因功能的注釋及分類分析。可進(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ù))等。通過(guò)與NR、Swiss～Prot、TrEMBL數(shù)據(jù)庫(kù)進(jìn)行對(duì)比，可以對(duì)CDS序列進(jìn)行預(yù)測(cè)；山莨菪的轉(zhuǎn)錄本序列與相近物種的近似情況，以及山莨菪同源序列的功能信息可以通過(guò)與NR庫(kù)對(duì)比得到；對(duì)山莨菪的基因進(jìn)行GO分類后，可以便于統(tǒng)計(jì)基因在生物過(guò)程、細(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(見(jiàn)表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條)。無(wú)論是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ì)和注釋(見(jiàn)表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%(見(jiàn)表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.無(wú)機(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)(見(jià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條，而注釋出來(lái)與染色體結(jié)構(gòu)、核結(jié)構(gòu)、細(xì)胞骨架、細(xì)胞周期等相關(guān)序列是較少的，有的甚至只有不到10條。

2.4 NR功能注釋

通過(guò)與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è)物種中。在所有的注釋中，注釋為與陽(yáng)芋(Solanumtuberosum)相關(guān)的Unigenes最多，有8 598條(21.70%)，其次為辣椒(Capsicumannuum)和煙草(Nicotianatabacum)，分別占16.24%和13.45%，比例所占較大的前幾種同科植物均可說(shuō)明與山莨菪具有較高的序列同源性。

2.5 GO功能注釋及分類

為全面描述山莨菪基因和基因產(chǎn)物的屬性，了解山莨菪表達(dá)基因的功能分布狀況及其所代表的生物學(xué)意義，將Unigenes通過(guò)與Swissprot、TrEMBL數(shù)據(jù)庫(kù)的比對(duì)，得到GO功能注釋信息，注釋后得到29 309個(gè)Unigenes，對(duì)其進(jìn)行統(tǒng)計(jì)分析后(表3)，可分為分子功能、生物學(xué)過(guò)程和細(xì)胞組分3個(gè)大類，分別有16、24、22個(gè)子類，共有62個(gè)子類，其中結(jié)合、催化活性、細(xì)胞過(guò)程、代謝過(guò)程、細(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è)子類(見(jiàn)圖3)。在4種代謝大類中，代謝相關(guān)的通路獲得注釋最多，為2 640個(gè)，占比49.45%，其次分別是遺傳信息處理、細(xì)胞過(guò)程、環(huán)境信息處理，分別獲得了1 438(26.93%)、656(12.29%)、605(11.33%)個(gè)注釋。23個(gè)亞類中，與代謝相關(guān)的通路最多，有12條，包括氨基酸代謝、碳水化合物代謝、其他次生物質(zhì)代謝、能量代謝、糖生物合成與代謝、脂類物質(zhì)代謝、輔助因子和維生素代謝、其他氨基酸代謝、萜類化合物和聚酮化合物代謝、核苷酸代謝、輔助因子代謝、概觀；與遺傳信息處理、細(xì)胞過(guò)程、環(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倍略多(見(jiàn)圖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ù)，每百萬(wàn)堿基中出現(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àn)資源植物的拼接結(jié)果，如芝麻Sesamumindicum(N50為1 006 bp，平均長(zhǎng)度683 bp)[24]、樟樹(shù)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無(wú)法獲得注釋，這可能是因?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ì)之后，注釋為陽(yáng)芋同源種的Unigenes有8 598條，遠(yuǎn)高于其他幾個(gè)被注釋物種，出現(xiàn)這種情況，可能是由于山莨菪與陽(yáng)芋的進(jìn)化史和生活史較為接近，可以說(shuō)明山莨菪與陽(yáng)芋的親緣關(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ù)為主。通過(guò)SNP分析與SSR分析，可以構(gòu)建較高分辨率的山莨菪遺傳圖譜，有助于提高目的基因定位的準(zhǔn)確性，同時(shí)促進(jìn)山莨菪群體遺傳學(xué)、比較基因組學(xué)和分子系統(tǒng)學(xué)的發(fā)展。