藥用資源植物山莨菪的轉錄組信息分析

張 雨 夏銘澤 張發起

(1.中國科學院高原生物適應與進化重點實驗室，中國科學院西北高原生物研究所，西寧 810001； 2.中國科學院大學，北京 100039)

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

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

鑒于此，本研究利用高通量測序技術，對山莨菪進行轉錄組測序分析，得到大量測序數據，結合生物信息學方法，拼接、組裝和注釋測序序列，對山莨菪的基因結構、表達水平和差異、基因富集進行分析，以期為山莨菪轉錄組和基因組水平的研究提供數據，同時也對山莨菪的保護及合理開發利用奠定分子生物學基礎。

1 材料與方法

1.1 材料

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

1.2 方法

1.2.1 山莨菪RNA提取與測序

提取山莨菪的RNA，在對其RNA進行純化和質量分析之后，構建cDNA文庫，采用Illumina HiseqTM測序平臺進行測序。

1.2.2 序列的拼接與組裝

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

①運行seqtk-trinity程序，將fq數據轉換成fa數據，并生成both.fa；

②運行Jellyfish程序，生成K-mer目錄；

③運行Inchworm程序，組裝reads初步得到contigs；

④運行Chrysalis程序，將上一步生成的contigs進行聚類分析，對每個類構建Bruijn圖；

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

⑥得到Trinity.fasta文件，其haeder信息包含了轉錄本名稱、長度、De Brujin圖重構途徑等；

⑦對拼接組裝得到的轉錄本去冗余。

在此基礎上，統計所得轉錄本的各項信息，如GC含量、長度分布狀況、所含isoform數目分布狀況等，為后續分析做準備。

1.2.3 基因功能注釋

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

2 結果與分析

2.1 轉錄本組裝與CDS預測分析

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

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

表1 拼接結果統計

2.2 功能注釋及COG分類

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

表2 Unigenes注釋結果

圖2 山莨菪轉錄組Unigenes的KOG功能分布圖 A.RNA加工和修飾RNA；B.染色體結構和動力學；C.能源生產與轉化；D.細胞周期調控、細胞分裂、染色體分離；E.氨基酸轉運和代謝；F.核酸轉運和代謝；G.碳水化合物轉運和代謝；H.輔酶轉運和代謝；I.脂類轉運和代謝；J.翻譯、核糖體結構和生物發生；K.轉錄；L.復制、重組和修飾；M.細胞壁/細胞膜生物發生；N.細胞活性；O.翻譯后修飾、蛋白轉運；P.無機離子轉運和代謝；Q.次生代謝物合成、轉運和代謝；R.只有一般功能預測；S.未知功能；T.信號傳遞機制；U.細胞間運輸、分泌物和囊泡運動；V.防御機制；W.細胞外結構；Y.核結構；Z.細胞骨架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蛋白質庫進行比對，有13 110條Unigenes獲得注釋，占總Unigenes的18.35%，Unigenes根據其功能被分為了26個類別，都與山莨菪最基本的生命活動相關(見圖2)。其中，在已經明確生物學功能的Unigenes中(占94.65%)，信號傳遞機制最多，有1 873條Unigenes，占比為12.76%，其次是只有一般功能預測和翻譯后修飾，分別為1 703(11.60%)和1 580(10.76%)。與轉錄、翻譯、化合物運輸、能量生產轉化、代謝等功能相關的Unigenes處于相對較多的數量，數量在500～900條，而注釋出來與染色體結構、核結構、細胞骨架、細胞周期等相關序列是較少的，有的甚至只有不到10條。

2.4 NR功能注釋

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

2.5 GO功能注釋及分類

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

圖3 山莨菪轉錄組Unigenes的KEGG功能注釋分布統計圖 A.細胞生長和死亡；B.細胞運動；C.細胞群體；D.運輸和分解代謝；E.膜運輸；F.信號轉導；G.信號分子和互作作用；H.折疊排序與退化；I.復制和修復；J.轉錄；K.翻譯；L.氨基酸代謝；M.其他次級代謝產物合成；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代謝庫(evalue<0.000 01)為參考，對代謝通路進行統計和分類分析。有3 679條Unigenes被注釋，這3 679條被注釋的基因根據涉及到的代謝通路分為4大類，23個子類(見圖3)。在4種代謝大類中，代謝相關的通路獲得注釋最多，為2 640個，占比49.45%，其次分別是遺傳信息處理、細胞過程、環境信息處理，分別獲得了1 438(26.93%)、656(12.29%)、605(11.33%)個注釋。23個亞類中，與代謝相關的通路最多，有12條，包括氨基酸代謝、碳水化合物代謝、其他次生物質代謝、能量代謝、糖生物合成與代謝、脂類物質代謝、輔助因子和維生素代謝、其他氨基酸代謝、萜類化合物和聚酮化合物代謝、核苷酸代謝、輔助因子代謝、概觀；與遺傳信息處理、細胞過程、環境信息處理相關的通路較少，分別有4條、4條、3條。在這23個子類代謝途徑中，翻譯獲得的注釋最多，為595條，占11.14%，其次為信號轉導、碳水化合物代謝，折疊排序和退化，分別為576條(10.79%)、479條(8.97%)和428條(8.01%)。

表3 山莨菪轉錄組Unigenes的GO功能分類統計

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

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

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

表4 山莨菪Unigenes數量最多的13個代謝通路

Table 4 Top thirteen metabolic pathways involved inA.tanguticusunigenes

代謝通路Pathway ID代謝通路PathwayUnigenes數量及比例Number and propotion of unigenesko03010核糖體Ribosome277(3.92%)ko01200碳代謝Carbon metabolism187(2.64%)ko04075植物激素信號轉導Plant hormone signal transduction173(2.45%)ko01230氨基酸合成Biosynthesis of amino acids167(2.36%)ko03040剪接Spliceosome158(2.23%)ko04141內質網蛋白加工Protein processing in endoplasmic reticulum156(2.20%)ko00230嘌呤代謝Purine metabolism137(1.94%)ko03013RNA轉運 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檢測分析

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

表5 山莨菪藥用活性成分代謝通路及基因統計

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

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

3 討論

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

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

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

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

轉錄組測序分析后，共檢測出31 382個SNP多態位點，兩種突變類型中以轉換突變為主，占2/3；在SSR分析中，獲得六種SSR重復類型，其中以單堿基重復、二堿基重復和三堿基重復為主。通過SNP分析與SSR分析，可以構建較高分辨率的山莨菪遺傳圖譜，有助于提高目的基因定位的準確性，同時促進山莨菪群體遺傳學、比較基因組學和分子系統學的發展。