唐賢 丁祥 董明明 朱淼 宋志強(qiáng) 侯怡鈴

（1. 西華師范大學(xué)生命科學(xué)學(xué)院 西南野生動(dòng)植物資源保護(hù)教育部重點(diǎn)實(shí)驗(yàn)室，南充 637009；2. 西華師范大學(xué)環(huán)境科學(xué)與工程學(xué)院，南充 637009）

食用菌子實(shí)體的形成和發(fā)育的分子機(jī)理對食用菌的產(chǎn)量和品質(zhì)具有關(guān)鍵的影響。近年來，我國食用菌產(chǎn)業(yè)發(fā)展迅速，香菇、杏鮑菇等已實(shí)現(xiàn)機(jī)械化生產(chǎn)，但目前一些重要食用菌仍主要依賴于野生資源，野生資源過度開發(fā)及環(huán)境的變化導(dǎo)致野生食用菌類資源量急劇下降，甚至一些種類瀕臨滅絕［1-2］。轉(zhuǎn)錄組測序能夠反應(yīng)不同生命階段的基因的表達(dá)、功能，能對基因進(jìn)行注釋及標(biāo)記［3］，已廣泛應(yīng)用于探索大型真菌發(fā)育相關(guān)基因及關(guān)鍵分子機(jī)制。榮成博等［4］通過轉(zhuǎn)錄組測序發(fā)現(xiàn)編碼苯丙氨酸裂解酶的基因（CL4509）和編碼乙醇氧化酶的基因與白靈側(cè)耳的子實(shí)體發(fā)育相關(guān)。吳小梅等［5］分析了雙孢蘑菇3個(gè)發(fā)育時(shí)期的轉(zhuǎn)錄組測序結(jié)果，發(fā)現(xiàn)雙孢蘑菇子實(shí)體發(fā)育與氨基酸代謝相關(guān)，能通過MAPK 信號(hào)途徑傳遞環(huán)境脅迫信號(hào)。施肖堃等［6］分析雙孢蘑菇As2796子實(shí)體4個(gè)發(fā)育階段的差異基因，發(fā)現(xiàn)疏水蛋白在其子實(shí)體發(fā)育與后熟過程中有重要的作用。李帆等［7］分析了刺芹側(cè)耳3個(gè)發(fā)育時(shí)期的轉(zhuǎn)錄組測序結(jié)果，發(fā)現(xiàn)菌絲期的碳代謝和氨基酸代謝、原基期的脂肪酸代謝和抗氧化酶、子實(shí)體時(shí)期的環(huán)境因子影響刺芹側(cè)耳發(fā)育。探索食用菌子實(shí)體的形成和發(fā)育的分子機(jī)制有助于實(shí)現(xiàn)野生食用菌的人工栽培生產(chǎn)。

棒瑚菌（Clavariadelphus pistillaris（Fr.）Donk）屬于擔(dān)子菌門，釘菇目，棒瑚菌科，棒瑚菌屬，子實(shí)體呈棒狀，不分枝，頂部鈍圓。主要分布于我國四川、云南等地［8］。棒瑚菌屬真菌能與殼斗科櫟屬植物形成菌根，保證樹木的正常生長，維持生態(tài)平衡［9］，在生態(tài)學(xué)上有很大的應(yīng)用潛力。棒瑚菌屬真菌含有豐富的微量元素，味道鮮美，可食用，因此它們適合作為營養(yǎng)食品和功能食品［10］。研究發(fā)現(xiàn)，從平截棒瑚菌中（Clavariadelphus truncatus.（Quel.）Donk.）提取的棒瑚菌酸（Clavaric acid）能作為法呢基轉(zhuǎn)移酶（FPTase）的抑制劑［11］，其菌絲體培養(yǎng)物具有廣泛的抗菌性能［12］。棒瑚菌屬是具有巨大發(fā)展?jié)摿Φ奶烊凰幬镔Y源。棒瑚菌屬共發(fā)現(xiàn)了22種，在我國僅發(fā)現(xiàn)9種，國內(nèi)外學(xué)者對棒瑚菌屬的研究主要在形態(tài)學(xué)、解剖學(xué)、顯微化學(xué)和開發(fā)應(yīng)用等方面［13-14］。目前，棒瑚菌尚未實(shí)現(xiàn)人工栽培，野生產(chǎn)量稀少，其獨(dú)特的形態(tài)在大型真菌中獨(dú)樹一幟。因此，分析棒瑚菌子實(shí)體不同發(fā)育時(shí)期的基因表達(dá)情況、差異表達(dá)基因，探索其子實(shí)體不同發(fā)育時(shí)期關(guān)鍵的功能基因和代謝通路，對棒瑚菌資源的保護(hù)與開發(fā)應(yīng)用具有一定的指導(dǎo)意義。

本文通過RNA-Seq測序技術(shù)，比較分析采自四川省阿壩藏族羌族自治州小金縣棒瑚菌子實(shí)體幼菇期和成熟期差異表達(dá)的基因，探索棒瑚菌發(fā)育過程中與子實(shí)體生長發(fā)育有關(guān)的功能基因及代謝機(jī)制，為棒瑚菌子實(shí)體發(fā)育的進(jìn)一步研究及食藥用資源開發(fā)提供一定的理論基礎(chǔ)。

1 材料與方法

1.1 材料

本實(shí)驗(yàn)的棒瑚菌采自四川省阿壩藏族羌族自治州小金縣，選取幼菇期（CP-1）與成熟期（CP-2）的棒瑚菌子實(shí)體作為研究材料。用無菌水將新鮮完整子實(shí)體清洗（重復(fù)3次），用無菌剪刀將子實(shí)體剪碎后，立即轉(zhuǎn)移至液氮中保存供下一步實(shí)驗(yàn)。

1.2 方法

1.2.1 樣品總RNA的提取及檢測 按照OMEGA（OMEGA's fungal total RNA extraction kit） 真 菌 總RNA提取試劑盒的使用說明提取棒瑚菌子實(shí)體中總RNA。稱取50 mg經(jīng)液氮研磨加入裂解液，離心；轉(zhuǎn)移上清加入0.5倍體積無水乙醇，離心棄濾液；加 Buffer I離心，棄去濾液；加入Buffer II離心棄去濾液（重復(fù)2次）；加入 DEPC Water離心。通過1%的瓊脂糖電泳和Agilent 2100檢測RNA樣品的完整性，用NanoDrop檢測RNA的純度，OD260/OD280應(yīng)在 1.8-2.2 之間，OD260/ OD230＞2［15］。

1.2.2 生物信息學(xué)分析 樣品總RNA檢測合格后送至北京諾禾致源生物信息科技有限公司進(jìn)行cDNA文庫構(gòu)建，庫檢合格后進(jìn)行RNA-Seq測序。

將原始測序數(shù)據(jù)（Raw data）進(jìn)行過濾篩選，得到分析數(shù)據(jù)Clean reads。利用Trinity［16］對分析數(shù)據(jù)Clean reads進(jìn)行拼接，獲得最終轉(zhuǎn)錄本序列，選取每個(gè)基因中最長的轉(zhuǎn)錄本作為非重復(fù)獨(dú)立基因Unigene。將非重復(fù)獨(dú)立基因Unigene在KOG/COG數(shù)據(jù)庫進(jìn)行功能注釋。利用RSEM［17］軟件進(jìn)行基因表達(dá)定量分析（FPKM＞0.3）。采用 DESeq［18］進(jìn)行基因差異性表達(dá)分析，閾值為padj＜0.05，并對差異基因進(jìn)行GO［19］富集分析和KEGG［20］富集分析。

2 結(jié)果

2.1 轉(zhuǎn)錄組的數(shù)據(jù)質(zhì)量分析

CP-1與CP-2進(jìn)行RNA-Seq測序后分別獲得原始數(shù)據(jù)（Raw reads），去除帶接頭的reads和低質(zhì)量的reads后分別獲得8.78 G（CP-1）和8.28 G（CP-2）分析數(shù)據(jù)（Clean reads）。綜合3次測序結(jié)果的平均值，其 中 CP-1的 Q20（%） 為 97.60%，Q30（%） 為94.12%，GC含量為57.75%，CP-2的Q20（%）為97.65%，Q30（%）為94.03%，GC含量為57.70%，堿基錯(cuò)誤率均為0.01%。拼接分析數(shù)據(jù)Clean reads后一共獲得151 367個(gè)轉(zhuǎn)錄本（Transcript），101 291個(gè)不同長度的非重復(fù)獨(dú)立基因 Unigene，占總轉(zhuǎn)錄本的66.92%。以上結(jié)果說明該測序結(jié)果可靠性較高，可用于下一步分析研究（表1，圖1）。

2.2 轉(zhuǎn)錄組基因KOG功能注釋及差異基因GO富集分析

轉(zhuǎn)錄組基因KOG功能注釋結(jié)果（表2）顯示，共有19 204個(gè)Unigene在KOG數(shù)據(jù)庫成功注釋，共涉及26類。非重復(fù)獨(dú)立基因Unigene數(shù)量最多的集中在僅通用功能預(yù)測（2 375條），翻譯后修飾、蛋白質(zhì)周轉(zhuǎn)、伴侶蛋白（2 194條）和翻譯、核糖體結(jié)構(gòu)、生物發(fā)生（2 030條）等，最少的是細(xì)胞運(yùn)動(dòng)（10條）。該結(jié)果顯示，CP-1和CP-2中基因所涉及到的KOG功能類別比較全面復(fù)雜，包括蛋白質(zhì)合成相關(guān)過程。值得注意的是，有800條Unigene的功能是未知的，且存在未命名蛋白（1條）。由此可見，棒瑚菌還有大量基因有待于進(jìn)一步研究。

表1 棒瑚菌測序轉(zhuǎn)錄組的數(shù)據(jù)質(zhì)量統(tǒng)計(jì)

用Goseq方法對CP-1和CP-2的差異基因進(jìn)行了GO富集分析，統(tǒng)計(jì)后結(jié)果（圖2）顯示，共產(chǎn)生 2 169個(gè)生物過程（Biological process），517個(gè)細(xì)胞成分（Cellula component）和1 076個(gè)分子功能（Molecular function）。差異基因GO富集結(jié)果提示，棒瑚菌在不同生長階段（CP-1和CP-2）的生物代謝存在一定的差異，同時(shí)提示，生物過程（Biological process）中的新陳代謝過程、細(xì)胞成分（Cellula component）中的細(xì)胞和細(xì)胞組分和分子功能（Molecular function）中的氧化還原酶活性是棒瑚菌子實(shí)體發(fā)育過程中的關(guān)鍵基因門類。

2.3 不同發(fā)育時(shí)期基因表達(dá)水平分析

采用RSEM軟件對基因表達(dá)進(jìn)行定量分析，F(xiàn)PKM超過60的基因高水平表達(dá)，而在0-0.3 FPKM的基因?yàn)榈退奖磉_(dá)或根本不表達(dá)?；虮磉_(dá)水平分析結(jié)果顯示，棒瑚菌在不同生長階段（CP-1和CP-2）基因表達(dá)水平的分布總體相似，但仍存在一定的差異。在棒瑚菌子實(shí)體幼菇期（CP-1）有39 140個(gè)基因表達(dá)（FPKM ≥ 0.3），約占基因總數(shù)的76.38%，其中，2 458個(gè)基因高表達(dá)（FPKM ＞60），約占基因總數(shù)的4.96%；而在棒瑚菌子實(shí)體成熟期（CP-2），有35 188個(gè)基因表達(dá)（FPKM ≥ 0.3），約占基因總數(shù)的68.5%，其中，2495個(gè)基因高表達(dá)（FPKM ＞ 60），約占基因總數(shù)的4.86%。

表2 棒瑚菌Unigene的KOG功能注釋統(tǒng)計(jì)

圖2 差異基因GO富集分析

進(jìn)一步分析發(fā)現(xiàn)，棒瑚菌子實(shí)體幼菇期（CP-1）和成熟期（CP-2）分別均有8個(gè)基因?yàn)闃O高表達(dá)（FPKM ＞ 5 000）。棒瑚菌子實(shí)體幼菇期（CP-1）的8個(gè)極高表達(dá)基因?yàn)椋篊cmD，HFB，RP-L41，PPO，MhfA，CuBPs，Hed，Hyd（表 3），在棒瑚菌子實(shí)體成熟期（CP-2）的8個(gè)基因極高表達(dá)基因?yàn)椋篊cmD，GH2，Ph-zid，PPO，CuBPs，MhfA，Hed，Hyd（表 4）。

編碼多酚氧化酶基因PPO和編碼疏水蛋白基因Hyd在CP-1和CP-2看均極高表達(dá)，hfb和RP-L41僅在CP-1中高表達(dá)，hfb為編碼真菌疏水蛋白的基因，該基因僅在CP-1中極高表達(dá)，說明在棒瑚菌的發(fā)育過程中至少有兩種疏水蛋白參與細(xì)胞壁的生物合成；RP-L41為編碼核糖體大亞基L41的基因，參與核糖體的生物合成及轉(zhuǎn)錄翻譯，提示棒瑚菌在成熟期（CP-2）核糖體的代謝發(fā)生了變化。而GH2和Ph-zid僅在CP-2 中高表達(dá)，編碼糖苷水解酶家族2的基因GH2僅在CP-2中極高表達(dá)，而編碼糖苷水解酶家族2的基因GH2和編碼苯丙氨酸拉鏈的基因Ph-zid僅在CP-2中極高表達(dá)，以上篩選出的高表達(dá)基因?qū)沂景艉骶诓煌L發(fā)育期的關(guān)鍵基因提供了一定的科學(xué)依據(jù)。

2.4 不同發(fā)育時(shí)期差異表達(dá)基因（DEGs）的分析

使用DESeq對CP-1和CP-2的差異表達(dá)基因（DEGs）進(jìn)行分析，并設(shè)定顯著差異表達(dá)的閾值Padj ＜ 0.05。棒瑚菌不同發(fā)育時(shí)期差異表達(dá)基因（DEGs）的分析結(jié)果顯示，棒瑚菌子實(shí)體幼菇期（CP-1）和成熟期（CP-2）的差異基因有6 625個(gè)，包括3 949個(gè)上調(diào)基因和2 676個(gè)下調(diào)基因（圖3），其中有2 290個(gè)基因只在幼菇期（CP-1）中表達(dá)，1 231個(gè)基因只在成熟期（CP-2）中表達(dá)。

表3 CP-1組中基因的定量表達(dá)（FPKM＞5 000）

在| log2fold change|＞ 6.5閾值內(nèi)的差異基因中，27個(gè)基因被上調(diào)，包括COX3、CYTB、wrbA、ND1、ARD1、EF3、ND4、SDHA、RP-S3e、RP-L5、RP-S5、EEF1A、ATP1、PDIA1和ARF1等；14個(gè)基因被下調(diào)，包 括 E1.11.1.5、PDIA1、htpG、UBE2D-E、CaM、RP-S9e、RP-L19e、ALDO和FDH等（表5-6）。

進(jìn)一步分析顯示，上調(diào)基因中，涉及一些氧化呼吸電子傳遞鏈的輔基基因，例如，編碼電子傳遞鏈中復(fù)合體Ⅰ的輔基的基因ND1和ND4，復(fù)合體Ⅱ的輔基的基因SDHA，復(fù)合體Ⅲ的輔基的基因CYTB，復(fù)合體Ⅳ的輔基的基因COX3，ATP合酶輔基的基因ATP1，提示氧化磷酸化途徑的基因出現(xiàn)上調(diào)。在下調(diào)基因中，包括編碼鈣調(diào)蛋白基因CaM和編碼甲酸脫氫酶基因FDH。

表4 CP-2組中基因的定量表達(dá)（FPKM＞5 000）

圖3 CP-1與CP-2基因差異表達(dá)分析火山圖

2.6 不同發(fā)育時(shí)期代謝通路KEGG分析

KEGG（Kyoto encyclopedia of genes and genomes）是有關(guān)代謝通路的主要公共數(shù)據(jù)庫［28］。棒瑚菌子實(shí)體幼菇期（CP-1）和成熟期（CP-2）差異基因KEGG Pathway富集結(jié)果（表7）顯示，有2 678個(gè)差異表達(dá)的基因成功注釋在104條通路中。其中，在核糖體（438個(gè)差異基因，16.36%）、糖酵解/糖異生途徑（83個(gè)差異基因，3.10%）、三羧酸循環(huán)（52個(gè)差異基因，1.94%）、氧化磷酸化（111個(gè)差異基因，4.14%）、內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工（118個(gè)差異基因，4.41%）、蛋白酶體（49個(gè)差異基因，1.83%）和MAPK信號(hào)通路（52個(gè)差異基因，1.94%）中，差異基因富集顯著。

本研究中三羧酸循環(huán)（Citrate cycle，TCA cycle）途徑（圖4）上調(diào)基因，包括pdhc（K00627）、DLST（K00658）、pdhD（K00382）、IDH3（K00030）、CS（K01647）、E4.2.1.2（K01676）、SDHA（K00234）、MDH1（K00025）和ACO（K01681）。MAPK信號(hào)通路（MAPK signaling pathway）（圖5）上調(diào)基因，包括Cdc42（k04393）、PAK1（k04409）、MKK（K11226）、PKC1（k02677）、MKK1_2（K08294）、MAPK7（K04464）、YPD1（K11232）和GNG（K07973）；下調(diào)基因包括STE12（K11215）和SSK2（K11230）。

此外，幼菇期（CP-1）與成熟期（CP-2）差異基因還顯著富集在核糖體、內(nèi)質(zhì)網(wǎng)中的蛋白質(zhì)加工、

蛋白酶體途徑，顯示棒瑚菌子實(shí)體幼菇期（CP-1）與成熟期（CP-2）蛋白質(zhì)合成途徑亦存在差異。

表6 差異基因的表達(dá)：下調(diào)（|log2fold change| ＞6.5）

表7 差異基因KEGG pathway富集分析

圖4 三羧酸循環(huán)

3 討論

大型真菌子實(shí)體發(fā)育過程涉及物質(zhì)能量代謝、細(xì)胞分化以及信號(hào)轉(zhuǎn)導(dǎo)等生命活動(dòng)。本研究選用四川省小金縣棒瑚菌幼菇期（CP-1）和成熟期（CP-2）的子實(shí)體進(jìn)行RNA-Seq測序并進(jìn)行生物信息學(xué)分析，不同發(fā)育時(shí)期基因表達(dá)水平分析發(fā)現(xiàn)，編碼多酚氧化酶的基因PPO在CP-1和CP-2中均極高表達(dá)，可能與棒瑚菌子實(shí)體的褐變及抗病性有關(guān)［21-22］。編碼疏水蛋白的基因Hyd在CP-1和CP-2中均極高表達(dá)，有利于孢子的擴(kuò)散和營養(yǎng)菌絲的分散，還能保護(hù)細(xì)胞免受惡劣環(huán)境的侵害［23-24］。而編碼糖苷水解酶家族2的基因GH2僅在CP-2中極高表達(dá)，說明GH2作為糖代謝的關(guān)鍵基因參與寡糖和多糖的合成，修飾和降解［25］。編碼苯丙氨酸拉鏈的基因Ph-zid僅在CP-2中極高表達(dá)，說明CP-2能在沒有細(xì)胞外配體的情況下激活和調(diào)節(jié)酪氨酸激酶的活性［26］。比較分析CP-1和CP-2差異表達(dá)的基因，發(fā)現(xiàn)上調(diào)基因主要富集在氧化呼吸電子傳遞鏈，說明棒瑚菌在幼菇期（CP-1）需要大量能量維持正常的生長發(fā)育，同時(shí)大量的耗氧代謝過程提示，環(huán)境含氧量對棒瑚菌的生長發(fā)育有重要影響。在下調(diào)基因中，編碼鈣調(diào)蛋白基因CaM下調(diào)，說明棒湖菌在成熟期（CP-2），細(xì)胞內(nèi)多種生理活動(dòng)如細(xì)胞增殖和細(xì)胞周期正在減弱［27］；編碼甲酸脫氨酶基因FDH下調(diào)，提示棒瑚菌在發(fā)育過程中甲酸含量增加，為維持一碳化合物穩(wěn)態(tài)，成熟期（CP-2）甲酸代謝過程增強(qiáng)［28］，初步篩選出了與棒瑚菌生長發(fā)育有關(guān)的代謝途徑，包括三羧酸循環(huán)、氧化磷酸化、核糖體途徑、蛋白酶體途徑以及MAPK信號(hào)通路等。

圖5 MAPK信號(hào)通路

三羧酸循環(huán)不僅是細(xì)胞能量的來源，而且能產(chǎn)生其他生物大分子合成的前體，更是生物體內(nèi)3大營養(yǎng)物質(zhì)代謝的樞紐［29］。本研究中三羧酸循環(huán)中的編碼丙酮酸脫氫酶的基因PDH、蘋果酸脫氫酶的基因MDH、檸檬酸合酶的基因CS、異檸檬酸脫氫酶的基因IDH、琥珀酸脫氫酶的基因SDH均在CP-1中表達(dá)上調(diào)，說明三羧酸循環(huán)在棒瑚菌幼菇期（CP-1）能量代謝的中起關(guān)鍵作用；而磷酸烯醇丙酮酸羧化激酶基因PCKA在棒瑚菌子實(shí)體成熟期（CP-2）中表達(dá)表達(dá)上調(diào)，說明棒瑚菌子實(shí)體成熟期（CP-2）的糖異生過程比幼菇期（CP-1）更為活躍，該過程與寡糖和多糖的合成、修飾等生理生化過程有關(guān)。MAPK信號(hào)通路是細(xì)胞內(nèi)一種重要的信號(hào)轉(zhuǎn)導(dǎo)系統(tǒng)，通過MAPK激酶激酶（MAP kinase kinase kinase，MKKK）、MAPK 激酶（MAP kinase kinase，MKK）和MAPK依次激活，共同調(diào)節(jié)細(xì)胞的生長、分化以及對環(huán)境的應(yīng)激等多種重要的細(xì)胞生理過程［30-31］。我們發(fā)現(xiàn)在MAPK信號(hào)通路中，細(xì)胞分裂控制蛋白42基因（Cdc42）表達(dá)上調(diào)，提示Cdc42可能參調(diào)控與菌絲的極性生長和細(xì)胞周期，說明真菌細(xì)胞形態(tài)轉(zhuǎn)換與Rho家族蛋白基因有關(guān)［32-33］，以上結(jié)果顯示，MAPK信號(hào)通路在棒瑚菌子實(shí)體的生長發(fā)育過程中有重要作用。

4 結(jié)論

在棒瑚菌的生長發(fā)育過程中，新陳代謝過程、細(xì)胞和細(xì)胞組分以及氧化還原酶活性是棒瑚菌子實(shí)體發(fā)育過程中的關(guān)鍵基因門類。多酚氧化酶，核糖體，糖苷水解酶家族等基因在棒瑚菌子實(shí)體的蛋白質(zhì)和細(xì)胞壁合成中起重要作用。氧化磷酸化和三羧酸循環(huán)是棒瑚菌子實(shí)體幼菇期（CP-1）生長發(fā)育的關(guān)鍵代謝途徑，且環(huán)境含氧量對棒瑚菌的生長發(fā)育有重要影響；同時(shí)MAPK信號(hào)通路是參與調(diào)節(jié)棒瑚菌子實(shí)體生長發(fā)育的關(guān)鍵信號(hào)通路，在調(diào)控菌絲的極性生長，細(xì)胞周期和真菌細(xì)胞形態(tài)轉(zhuǎn)換中起重要作用。