基于轉錄組測序的兩種豬苓菌絲體多糖代謝途徑比較分析

任 潔， 李太元， 李艷茹， 梁運江， 許廣波

(延邊大學農學院，吉林 延吉 133000)

豬苓[Polyporusumbellatus(Pers.) Fries.]是非褶菌目、多孔菌科、多孔菌屬的一種藥用真菌[1-4]。我國的豬苓資源主要分布于陜西、甘肅等西南省區以及東北的長白山區[5-7]。長白山區的豬苓在生物學特性上與陜西等地的豬苓有很大差異，俗稱雞爪苓[8-9]。長白山豬苓菌絲體的生長速度比較快，菌絲潔白，能夠形成菌核狀組織體，而陜西豬苓菌絲體長勢緩慢，易老化，不能形成菌核狀組織體。長白山豬苓是一種特異的豬苓種質資源[10-14]。

豬苓多糖和甾體是豬苓的主要藥用成分，且豬苓多糖是最早從豬苓菌核中分離出來的活性成分。現代醫學研究發現，豬苓多糖還具有抑制腫瘤生長，增強機體免疫力，降血糖等藥理作用[15-16]。目前對豬苓多糖的研究還是以多糖的分離純化及藥理作用等方面研究居多，在分子水平上對豬苓多糖合成機制的研究相對較少[17-22]。該研究利用Illumina高通量測序技術對長白山豬苓菌絲體以及陜西豬苓菌絲體進行轉錄組測序及分析，挖掘氨基糖與核苷酸糖代謝以及果糖與甘露糖代謝相關的功能基因，分析長白山豬苓菌絲體與陜西豬苓菌絲體的基因表達差異特性，為今后深入研究豬苓多糖的生物合成機制和關鍵性功能基因的克隆提供依據。

1 材料與方法

分別從長白山豬苓菌核和陜西豬苓菌核中通過組織分離法獲得長白山豬苓菌株(C)和陜西豬苓菌株(S)，將獲得的2種豬苓菌株進行活化擴繁，接種于裝有液體培養基的三角瓶中，置于25 ℃恒溫培養箱中進行培養并得到2種豬苓的菌絲體。每個菌株分別選取3份生長狀況最佳的菌絲體(標記成C1、C2、C3、S1、S2、S3，共6個樣品)作為該研究的試驗材料，送至上海生工生物工程技術服務有限公司進行建庫、測序及相關生物信息學分析。

2 結果與分析

2.1 測序數據質量評估及質控

對2種豬苓菌絲體共6個樣品進行轉錄組測序，共獲得原始測序序列330 834 156個。對測序得到的原始數據進行過濾后，得到Clean數據325 077 392個(表1)。由表1可知，S3、S2、S1、C3、C2、C1 6個樣品的Q20值分別為99.36%、99.36%、99.35%、99.20%、99.40%和99.38%。Q30值分別為97.33%、97.32%、97.29%、96.79%、97.45%和97.41%。GC相對含量分別為55.53%、55.71%、55.67%、55.77%、55.57%和55.45%。經過轉錄組測序得到的數據質量較高，可用于后續的生物信息學分析。

表1 長白山豬苓和陜西豬苓菌絲體各樣本QC數據統計總表Table 1 Statistics of QC data of each sample

2.2 Unigene功能注釋

2.2.1 NR庫比對注釋

通過與NR庫的比對，可以查看物種轉錄本序列與相近物種的近似情況以及同源序列的功能信息。將得到的20 999條Unigene與NR數據庫進行同源性比對注釋，有10 815條Unigene可以找到注釋結果(圖1)。其中，有8 670條Unigene屬于多孔菌科，排名前3的最為相近物種是污叉絲孔菌(Dichomitus squalens LYAD-421 SS1)、云芝栓孔菌(Trametes versicolor FP-101664 SS1)、朱紅栓菌(Trametes cinnabarina)，分別占41.51%、18.03%和15.12%。尚有2 145條(19.82%)Unigene比對到其他物種當中。

圖1 同源物種分布餅圖Fig.1 Pie chart of homologous species distribution

2.2.2 GO功能分類注釋

將拼接得到的所有Unigene與基因本體論數據庫(gene ontology，GO)進行比對和GO分類，統計基因在生物學過程(Biological Process)，細胞組分(Cellular Component)，分子功能(Molecular Function)3個類別的GO條目，結果如圖2所示。共有9 389個Unigene被注釋，占所有Unigene總數的44.71%，涉及分子功能、生物學過程和細胞組分3大類68個分支。其中，參與分子功能的基因有6 848個，涉及到20個分支，其中結合(binding)和催化活性(catalytic activity)包含的基因數量最多，分別有5 428和5 038條；與細胞組分有關的基因有7 494個，涉及22個分支，與細胞(cell)及細胞組分(cell part)相關的基因分別達到6 089、6 079條；參與生物學過程的基因有6 848個，涉及到26個分支，其中包括的細胞過程(cellular process)和代謝過程(metabolic process)涉及的基因數量最多，分別有5 739和5 090條。

圖2 GO功能分類圖Fig.2 GO function classification

2.2.3 KOG功能分類注釋

KOG分為26個group，將KOG注釋成功的基因按KOG的group進行分類，注釋結果如圖3所示。在把所有的Unigene注釋到蛋白質真核同源數據庫(eukaryotic orthologous groups，KOG)時，發現有4 976條Unigene對應25個KOG類別找到了對應的功能信息。根據統計分析，在25個KOG類別當中，一般功能預測類(General function prediction only)要比其他類更多，共有620條，占12.46%；其次為翻譯后修飾、蛋白質轉換、伴侶(Posttranslational modification, protein turnover, chaperones)以及翻譯、核糖體結構和生物合成(Translation, ribosomal structure and biogenesis)相關的基因，分別有564個和512個，占11.33%和10.29%；與細胞運動(Cell motility)有關的只有8條，比例僅占0.16%。

圖3 KOG功能分類圖Fig.3 KOG function classification

2.2.4 KEGG功能分類注釋

對基因做KO注釋后，根據KO與Pathway的聯系對其進行KEGG代謝通路分類。將拼接得到的所有Unigene KEGG進行比對分析(圖4)，顯示共有1 959條Unigene被注釋到，分屬于細胞過程(Cellular Processes)、環境信息處理(Environmental Information Processing)、遺傳信息處理(Genetic Information Processing)、代謝(Metabolism)和生物體系統(Organismal Systems)5大類33小類，共280條代謝通路。其中，細胞過程涉及的Unigene可分為細胞生長與死亡(Cell growth and death)、細胞運動(Cell motility)、細胞群體(Cellular community)、運輸和分解代謝(Transport and catabolism)4類；環境信息處理相關的Unigene可分為膜運輸(Membrane transport)、信號轉導(Signal transduction)、信號分子與相互作用(Signaling molecules and interaction)3類；遺傳信息處理涉及的Unigene可分為折疊分類和降級(Folding sorting and degradation)、復制和修復(Replication and repair)、轉錄(Transcription)、翻譯(Translation)4類；代謝相關的Unigene可分為氨基酸代謝(Amino acid metabolism)、其他次生代謝產物的生物合成(Biosynthesis of other secondary metabolites)、碳水化合物代謝(Carbohydrate metabolism)、能量代謝(Energy metabolism)等12類；生物體系統相關的Unigene可分為免疫系統(Immune system)、神經系統(Nervous system)、內分泌系統(Endocrine system)、循環系統(Circulatory system)等10類。

圖4 KEGG功能分類圖Fig.4 KEGG function classification

在KEGG中注釋最多的前15位代謝通路中(表2)，涉及核糖體(Ribosome的Unigene)最多，有156條；其次是氨基酸的生物合成(Biosynthesis of amino acids)，有97條；與氧化磷酸化(Oxidative phosphorylation)相關的Unigene達到90條；與碳代謝(Carbon metabolism)相關的Unigene達到85條；與RNA轉運(RNA transport)相關的Unigene有70條。

表2 KEGG注釋最多的前15位代謝通路Table 2 Top 15 metabolic pathways with the most annotations by KEGG

表2 KEGG注釋最多的前15位代謝通路Continue to table 2 Top 15 metabolic pathways with the most annotations by KEGG

2.3 表達差異分析

根據組間比較得到的差異基因繪制差異基因火山圖(圖5)。與陜西豬苓菌絲體相比，長白山豬苓菌絲體共有5 881個表達差異基因出現，其中,上調基因有2 616個，下調基因有3 265個。

圖5 比較組表達差異火山圖Fig.5 Volcano diagram of the difference in expression of the comparison group

2.4 差異表達基因pathway富集分析

長白山豬苓菌絲體和陜西豬苓菌絲體的Pathway顯著性富集分析結果(表4)，氨基糖和核苷酸糖代謝顯著富集。在具有明顯差異的12個pathway中，核苷酸切除修復(Nucleotide excision repair)、基礎切除修復途徑(Base excision repair)、植物-病原互作(Plant-pathogen interaction)代謝途徑中的基因在長白山豬苓菌絲體中全部下調表達；另外氨基糖和核苷酸糖代謝(Amino sugar and nucleotide sugar metabolism)、甘氨酸、絲氨酸和蘇氨酸的代謝(Glycine, serine and threonine metabolism)、淀粉和蔗糖代謝(Starch and sucrose metabolism)、溶酶體(Lysosome)等相關代謝途徑中的基因在長白山豬苓菌絲體中基本下調表達，表明長白山豬苓菌絲體生長發育初期細胞代謝水平不如陜西豬苓菌絲體旺盛。

表4 差異基因pathway顯著富集分析表Table 4 Significant enrichment analysis of differential gene pathway

2.5 基于KEGG注釋的差異基因代謝通路

通過對氨基糖和核苷酸糖代謝通路進行分析(表5)，結果顯示與氨基糖和核苷酸糖代謝相關的Unigene共有49個，涉及到28個基因產物，其中，有13條差異基因參與到8個基因產物的代謝過程中。其中，編碼幾丁質酶(chitinase)的Unigene最多，有8條(包含上調基因1個，下調基因5個)；其次是編碼幾丁質合酶(chitin synthase)的Unigene，有7條(包含下調基因1個)；另外，有2條Unigene(包含上調基因1個)編碼磷酸葡萄糖突變酶(phosphoglucomutase)、3條Unigene(包含下調基因1個)編碼己糖胺酶(hexosaminidase)，其余4條差異基因都為下調基因，分別編碼N-乙酰氨基葡萄糖-6-磷酸脫乙酰酶(N-acetylglucosamine-6-phosphate deacetylase)、UTP-葡萄糖-1-磷酸尿酸轉移酶(UTP-glucose-1-phosphate uridylyltransferase)、磷酸甘露糖突變酶(phosphomannomutase)和GDP-L-巖藻糖合酶(GDP-L-fucose synthase)。

表5 氨基糖和核苷酸糖代謝相關UnigeneTable 5 Unigene related to amino sugar and nucleotide sugar metabolism

對果糖和甘露糖代謝通路進行分析的結果表明(表6)，與果糖和甘露糖代謝相關的Unigene共有28個，涉及到20個基因產物，其中，有4條差異基因參與到8個基因產物的代謝過程中。編碼果糖二磷酸醛縮酶I類(fructose-bisphosphatealdolase, class I)的Unigene最多，有5條；其次是編碼醛還原酶(aldehyde reductase)和果糖二磷酸醛縮酶II類(fructose-bisphosphatealdolase, class II)的Unigene，有2條(分別包含上調基因1個)；其余2條差異基因都為下調基因，分別編碼磷酸甘露糖突變酶(phosphomannomutase)和GDP-L-巖藻糖合酶(GDP-L-fucose synthase)。

表6 果糖和甘露糖代謝相關UnigeneTable 6 Unigene related to fructose and mannose metabolism

3 討論

隨著轉錄組測序技術和相應分析軟件的成熟與完善，通過比較轉錄組研究來揭示食藥用真菌的發育相關基因和活性物質代謝調控解析得到了真菌研究者的青睞。王東等對銀耳(Tremellafuciformis)菌絲體進行轉錄組測序分析[23]，組裝得到17 008條unigene序列，比白蟻蘑菇和雙孢蘑菇得到的序列數量要多，對果糖與甘露糖代謝途徑相關基因分析，得到74條unigene，編碼21個酶。聶文強等對灰樹花(Grifolafrondosa)進行轉錄組測序分析[24]，得到115條unigene與多糖合成代謝有關。Huang等對桑黃(Phellinuslinteus)菌絲體進行轉錄組測序[25]，獲得25 811個Unigene，探討了桑黃甾醇的生物合成基因。該研究發現氨基糖和核苷酸糖代謝過程顯著富集，與其相關的差異基因共有13條，參與了8個基因產物的代謝過程。其中有1個上調基因和5個下調基因參與了幾丁質酶的合成，有1個下調基因參與編碼幾丁質合酶，有1個上調基因參與編碼磷酸葡萄糖突變酶，1個下調基因參與編碼己糖胺酶，其余4條差異基因都為下調基因，分別參與編碼N-乙酰氨基葡萄糖-6-磷酸脫乙酰酶、UTP-葡萄糖-1-磷酸尿酸轉移酶、磷酸甘露糖突變酶和GDP-L-巖藻糖合酶。與果糖和甘露糖代謝相關的差異基因有4條，參與到8個基因產物的代謝過程中。其中有2條上調差異基因分別編碼醛還原酶和果糖二磷酸醛縮酶II類，2條下調差異基因分別編碼磷酸甘露糖突變酶和GDP-L-巖藻糖合酶。

4 結論

該研究篩選出了2種豬苓菌絲體的差異表達基因，并對差異表達基因進行pathway富集分析以及豬苓多糖的相關代謝通路分析，結論如下。

1) 通過對轉錄組原始數據進行數據質控共得到20 999個Unigene。

2) 差異表達分析結果表明，長白山豬苓與陜西豬苓菌絲體相比共有5 881個表達差異基因出現，在長白山豬苓菌絲體中上調表達的基因有2 616個，下調表達的基因有3 265個。

3) 通過對氨基糖與核苷酸糖和果糖與甘露糖的代謝通路進行分析，發現在氨基糖與核苷酸糖代謝通路中共涉及到28個基因產物，與其相關的Unigene共有49條，其中有13條差異基因參與到8個基因產物的代謝過程中。在果糖與甘露糖代謝相關的Unigene共有28個，涉及到20個基因產物，其中有4條差異基因參與到8個基因產物的代謝過程中。