宋 鑫，蘇 浩，蔣詩雨，李 力*

1.河南工業(yè)大學 糧油食品學院，河南 鄭州 450001 2.四子王旗市場監(jiān)督管理局，內蒙古 烏蘭察布 011800

來源廣泛且具有抗氧化、調節(jié)腸道菌群等功能活性的麥麩膳食纖維屬于木質纖維素，由纖維素、半纖維素和木質素構成，復雜的物化構造限制了其作為生態(tài)學重要資源的發(fā)展和使用[1]。白腐菌中蘊含完整的細胞外木質素降解酶系統，是目前最高效的木質及膳食纖維素降解微生物之一[2]。實驗室前期通過篩選得到1株對小麥麩皮具有高改性效率的白腐菌菌株A.polytricha5.584，通過優(yōu)化后對小麥麩皮木質素的降解率達到68.72%[3-4]，需要在分子層面上對其降解麥麩纖維的機制進行深入研究。

轉錄組是某個群體或細菌在某一生長發(fā)育階段或各種特殊功能境況下轉錄得到的全部RNA的集合體，可以從總體高度深入研究基因組結構各種功用和基因組構造，闡明單個生物過程和生物過程中的基因分子機制。轉錄組測序法(RNA-seq)是研究生物功能基因組的重要方式，同時也是基因功能與構造研究的重要基礎與出發(fā)點[5]。劉家豪[6]對紅側耳菌的轉錄組進行測序，發(fā)現了與玉米秸稈纖維素、半纖維素和木質素降解有關的基因89個，其中70個為上調基因，證明了溶解性多糖單氧酶在玉米秸稈木質纖維素的分解中起著重要作用。

目前，使用野生型白腐真菌對麥麩纖維進行生物改性已經取得了較好的效果，但其降解麥麩木質纖維素的分子機理并不明晰，且白腐真菌以麥麩纖維作為單一碳源時，其胞外纖維素酶系的表達調控機制也鮮有報道，因此難以使用分子生物學方法對其進行定向突變和選育[7]。作者通過RNA-seq對白腐菌A.polytricha5.584在不同碳源下的差異表達基因開展功能注釋、GO與KEGG功能分類和蛋白網絡互作的分析，篩選出與麥麩纖維降解相關的基因，為在分子水平上揭示白腐真菌降解麥麩纖維的機制奠定理論基礎。

1 材料與方法

1.1 材料

麥麩：河南中鶴有限公司；菌株A.Polytricha5.584：中國國家普通微生物學菌種保存管理中心福建省三明真菌研究院，試驗前將菌株儲存于PDA(Potato Dextrose Agar)培養(yǎng)基上或25 ℃的培養(yǎng)基盒中。

Trizol、RNase-free：賽默飛世爾科技公司；異戊醇、乙醇(75%)：天津市科密歐化學試劑有限公司。

PDA培養(yǎng)基(g/100 mL)：PDA 3.7，KH2PO40.25，MgSO40.15，VB10.05，土霉素 0.02。

麥麩纖維液態(tài)培養(yǎng)基(g/100 mL)：麥麩纖維5.0，K2HPO40.1，FeSO40.005，KH2PO40.1，MgSO40.02，CaCl2·2H2O 0.01，(NH4)2SO40.2，MnSO40.002。

葡萄糖液態(tài)培養(yǎng)基(g/100 mL)：葡萄糖2.0，FeSO40.005，K2HPO40.1，KH2PO40.1，MgSO40.02，CaCl2·2H2O 0.01，(NH4)2SO40.2，MnSO40.002。

1.2 儀器與設備

HWS智能恒溫恒濕箱：寧波市東南儀表公司；JP-1500B高速萬能粉碎機：永康市久品工貿公司；JM-F50膠體磨：浙江省麥特隆機器有限公司；數顯恒溫式水浴鍋：金壇市佳美儀表公司；DHG-9011A電熱恒溫干燥盒：上海市精宏實驗儀器機械設備公司。

1.3 試驗方法

1.3.1 麥麩膳食纖維的制備

麥麩磨碎后過40目篩，將麥麩和蒸餾水按質量比1∶ 10配制，用膠體磨處理25 min；干燥后調節(jié)加水量(1∶ 10)，95 ℃蒸30 min；降溫至50 ℃以下，加入HCl調pH值為5.6，升溫至95 ℃；加入淀粉酶(1.5%)，攪動30 min；降溫至50 ℃以下，加入NaOH調pH值至9.0，攪拌均勻；加入堿性蛋白酶(3%)，攪動2 h，升溫至100 ℃滅酶；酶解的麥麩用自來水洗滌數遍，直到洗滌液不渾濁；置于烘箱中，60 ℃直至烘干。

1.3.2 菌絲培養(yǎng)及取樣

將保藏菌株注射在PDA平板培養(yǎng)基(9 cm)上，25 ℃培養(yǎng)14 d使菌絲完全長滿平板，取中央菌絲，分別注射至裝有100 mL麥麩纖維液態(tài)培養(yǎng)基、葡萄糖液態(tài)培養(yǎng)基的250 mL三角瓶中，封口后置于27 ℃、180 r/min的振蕩培養(yǎng)箱中。培養(yǎng)7 d后依次取培養(yǎng)液6 mL，室溫條件下4 000 r/min離心10 min，去除上清液，菌體用無菌水充分重懸3次，洗滌步驟重復5次以上，直至完全去除培養(yǎng)基雜質，菌體洗滌液干凈透明。得到不少于200 mg的活菌體，轉入冷凍存管封口，儲存于液體氮氣罐中。

1.3.3 RNA 提取與轉錄組測序

用液氮把菌株細胞研磨成粉末，轉移至2 mL試管中。加1.5 mL的Trizol，充分搖勻后，室溫靜置5 min，10 000 r/min離心5 min。向細胞上清液中加入200 μL異戊醇(上清液與異戊醇體積比24∶ 1)和1 mL裂解試劑，10 000 r/min離心10 min，把細胞上清液轉至等體積異丙醇的試管中，于-20 ℃的冰箱中冷卻1 h。13 600 r/min離心20 min后，用1 mL 75%乙醇結晶上清液，并靜置3～5 min。RNA用30～100 μL RNase-free水溶解。采集后的RNA通過NanoDrop系統檢測濃度，RNA的結構完整性通過Agilent2100生物分析儀測試，檢驗無誤后，開展菌株的轉錄組測序。

用mRNA富集法對總RNA進行加工。用含有OligodT的磁珠富集含有polyA尾巴的mRNA，然后再用打斷buffer將獲得的RNA片段化，用隨機的N6引物進行反轉錄，將人工合成的cDNA二鏈形成雙鏈DNA。將制備的雙鏈DNA末端補平并5′末端磷酸化，3′端成為凸起1個“A”的粘末端，再接上一條3′端有凸起“T”的鼓泡形連接器。把鏈接物質經過特異的引物實現PCR擴增。然后將PCR物質經過熱改性為單鏈結構，再用一種橋型引物使單鏈結構DNA環(huán)化，獲得單鏈結構環(huán)狀DNA文庫并上機檢測。

轉錄組測序委托深圳市華大基因技術服務公司，利用高通量測序平臺BGI-500，開展paired-end測序。為使后續(xù)的組裝和分析結果更為可信，使用華大自主開發(fā)的過濾軟件SOAPnuke(版本：v1.4.0；參數：-l 5 -q 0.5 -n 0.1)進行統計，并使用trimmomatic(版本：v0.36；參數：ILLUMINACLIP:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:50)進行過濾，過濾后的Clean reads用于后續(xù)的研究和分析。

1.3.4 基因信息的基礎分析

數據過濾方式：采用SOAPNuke進行統計，trimmomatic進行過濾，過程包括：除去含有連接的reads(連接污染)；除去未知堿基N濃度超過5%的reads；除去較低品質的reads(質量值小于10的堿基占該reads總堿基數的比率超過20%的reads)。

參照基因組的比對方式：首先采用Bowtie2 (版本：v2.2.5;參數：-q-p hred64-sensitive-dp ad 0-gbar 99999999-mp 1，1-np 1-score-min L，0，-0.1-p 16-k 200)[8]將clean reads比對到參考基因組順序，然后再采用RSEM(版本：v1.2.8)[9]估計基因組和轉錄本的發(fā)表水平。

基因信息De novo組裝：首先使用Trinity(版本：v2.0.6；參數：-min_contig_length 150-CPU 8-min_kmer_co v 3-min_glue 3-bfly_opts′-V 5-edge-thr=0.1-stderr′)對clean reads(去除PCR重復以提高組裝效率)進行De novo裝配，之后再通過Tgicl將已組裝的標本進行聚類分析除去數據的冗余，從而得到Unigene[10-11]。Trinity包括3個獨立模塊：Inchworm、Chrysalis、Butterfly，分別處理了大量reads。Trinity將reads建立為大量獨立的de Bruijn圖，然后對每個圖分別提取全長的轉錄本水平剪切亞型。De novo的組裝品質評價：基于組裝結果，通過BUSCO數據庫對Trinity組裝的轉錄本進行品質評價，通過與保守基因進行比較，進而說明轉錄組的組裝效果[12]。

編碼序列CDS(Coding sequence)的預測：使用Trinity推薦的軟件Transdecoder(版本：v3.0.1)識別Unigene中的候選編碼區(qū)，并首先獲得最長的開放閱讀框，之后再使用Blast比對Swiss-Prot和Hmmscan搜索Pfam蛋白同源順序，進而預測編碼區(qū)。具體步驟：(1)使用 TransDecoder.LongOrfs通過對Unigene.fa進行開放的讀碼框查找，以獲得更長讀碼框順序；(2)通過序列相似性用Diamond Blastp把Unigene.fa結果比對在Swiss-Prot數據庫系統上，然后再使用Hmmscan對Blast結果在pfam數據庫檢索；(3)利用TransDecoder.Predict估計Unigene.fa的編碼順序。

微衛(wèi)星序列標記方法：使用MISA(版本：v1.0;參數：1-12，2-6，3-5，4-5，5-4，6-4 100 150)對Unigene進行檢測，使用Primer3 (版本：v2.2.2)對檢測到的SSR 進行引物設計[13-14]。

1.3.5 基因表達量與表達量分布的分析

通過Bowtie2 (修訂版：v2.2.5；參數：-q-p hred64-sensitive-dpad 0-gbar 99999999-mp 1，1-np 1-score-min L，0，-0.1-p 16-k 200)將clean reads比對到基因組順序上，并且利用RSEM(修訂版：v1.2.8)測算各個樣本中的基因組表達水平。根據各個樣本的表達量FPKM信息，繪制密度圖。

1.3.6 差異基因的篩選與功能注釋

DEGseq(參數：差異倍數2倍以上且Q-value≤0.001)方法基于泊松分布，根據Wang等[15]描述的方式，進行了差異基因測定。RNA測序也可被視為隨機取樣的步驟，每次read在樣本中都獨立而均勻地采樣。在這種假設下，來自整個基因組(轉錄本)的read數目遵循二項分布(并且近似由泊松分布代替)。通過上述的數據模型，DEGseq給出了一種基于MA-plot的新方式，MA-plot是應用于芯片數據的統計分析工具，讓C1和C2顯示由2個樣品中得到的特定基因的read總量，按照Ci～binomial(ni，pi)，i=1，2分布，其中ni代表所有比對上的read總量，pi代表該基因的概率。DEGseq定義M=log2C1-log2C2和A=(log2C1+log2C2)/2，并證明在隨機抽樣假設下，給定A=a (a是A的觀察值)的條件下分布，M遵循近似正態(tài)分布。對MA圖上的每個基因都進行假設并檢驗H0： p1=p2與H1： p1p2。然后根據正態(tài)分布計算P, 參考Scheid等[16]的方法把P-values修正為Q-values。為增加DEGs的準確率，界定了差異倍數為2倍以上且Q-value≤0.001的基因，并遴選為顯著差異表達基因。

使用hmmscan(版本：v3.0)、Blast(版本：v2.2.23)和Blast2GO(版本：v2.5.0)對Trinity組裝得到的Unigene與七大功能數據庫(NR(Non-Redundant Protein Sequence Database)、NT(Non-Redundant Protein Sequence Database)、Swiss-Prot、KOG(Clusters of orthologous groups for eukaryotic complete genomes)、GO(Gene Ontology)、KEGG(Encyclopedia of Genes and Genomes)、Pfam(Protein Families Database))進行比對注釋。當基因功能注釋結果不同時，為保證生物學意義，按照NR、Swiss-Prot和KOG的順序將比對結果作為此基因的注釋。

2 結果與分析

2.1 轉錄組數據與質量控制評估結果

采用BGISEQ-500平臺一共測了38.18 Gb數據，重組和去冗余后共得出31 531個Unigene，總長度、平均總長度、n50及GC濃度依次為57 226 025、1 814、2 537 bp和52.42%。通過把Unigene比對到七大功能數據庫并加以注解，結果依次有26 646 (NR：84.51%)、5 590(NT：17.73%)、16 359(Swiss-Prot:51.88%)、15 404(KOG:48.85%)、18 475(KEGG:58.59%)、8 638(GO:27.40%)以及20 686(Pfam:65.61%)個Unigene獲得功能注釋。利用Transdecoder檢查出24 404個CDS，同樣還檢查出1 989個SSR散布于1 536個Unigene中。

轉錄組測序數據及與參考基因組的比對結果統計見表1。各樣品堿基質量值Q30均不小于86.65%。各樣品的Clean Reads與參考基因組的比對效率為 86.29%～89.49%，數據利用率較高，表明所選參考基因組能夠滿足后續(xù)分析的需求。

表1 測序數據及比對參考基因組結果統計Table 1 Results statistics of sequencing data and reference genome

2.2 麥麩纖維處理組與葡萄糖對照組差異表達基因與表達量分布分析

使用DEGseq軟件進行麥麩纖維組與葡萄糖組之間的差異性基因表達分析，以差異倍數≥2且Q-value≤0.001的基因體作為檢驗標準，得出了2組樣品中間的不同表達基因體集，DEGs(differentially expressed genes)共12 096條，當中上調基因8 214條，下調基因3 882條。使用hmmscan、Blast和Blast2GO軟件對所有DEGs在各個數據庫系統中執(zhí)行功能注釋的結果如表2所示，可知在各個數據庫系統中被注釋的DEGs總量有所不同，最后共計31 531個DEGs獲得注解，共計2 073個DEGs被所有數據庫系統同時注釋。

表2 注釋的差異表達基因數量統計Table 2 Statistics of annotated differentially expressed genes

對于各種注釋的差異表達基因數據，選取了KEGG、GO、NR、NT、Swiss-Prot基因注釋數據庫作出5維的差異表達基因的Venn圖(圖1)，由各樣本FPKM密度分布可以看出麥麩纖維液態(tài)培養(yǎng)下與葡萄糖液態(tài)培養(yǎng)下RNA的表達量密度曲線分布不同，但組間表達密度分布曲線相互接近，表明分組合理且組間存在差異，樣本表達量密度如圖2所示。

圖1 KEGG、GO、NR、NT、Swiss-Prot差異表達基因注釋維恩圖Fig.1 Venn diagram of differentially expressedgenes by KEGG, GO, NR, NT, and Swiss-Prot

圖2 樣本表達量密度圖Fig.2 Density map of sample expression

GO是一種國際標準化的生物基因活性分類系統，用來描述生命體中基因及其所翻譯蛋白質的多種屬性[17]，并把基因表達形成的不同功能歸納為與生物學過程、細胞組成和分子功能有密切聯系的類型。本研究對處理組和對照組通過GO注釋進行了統計分析，認為所獲得的基因轉錄組共計8 638個基因得到GO注釋，其中有3 415個與得到注釋的基因是不同基因，并且它們都被標記在了三大主要功能分類目的42個二級條目中(20+12+10)。在3 415個不同基因中，注釋在“生物過程”功能的二級功能“細胞新陳代謝步驟(GO:0008152)”中的不同基因表達數共有1 088個；注釋在“細胞組分”功能的二級功能“膜(GO:0016020)”中的差異基因數量有1 111個；注釋在“基因分子結構和功能”功能的二級功能“催化活性(GO:0003824)”中的差異基因數量最多，有1 818個。

對全部的3 415個功能差異基因進行GO注釋富集解析，篩選出了在三大主要功能分類中校正P<0.01的二級條目共26條，包括“生物過程”類別中的10個類目；“細菌組成”類別中的5個類目；“分子物質構成和功用”類別中的11個類目。綜合GO數據庫中提出的主要類目功能描述，可將這些類目分為兩組：一組是直接參與細胞能量新陳代謝流程，是白腐菌在麥麩纖維等基質環(huán)境生長過程中所必需的能量代謝過程；另一組則是參與木質纖維素等麥麩纖維基質氧化分解與代謝過程，包括糖酵解、三羧酸循環(huán)、蛋白酶的氧化與還原活性等4個基因本體功能類目，如表3所示。

表3 GO富集麥麩纖維降解相關類目統計Table 3 Statistics of GO enriched wheat bran fiber degradation related categories

KEGG是理解高等各種特殊功能和生物學體系(如細菌、生物學和人類)，從分子水平數據資料信息，尤其是大規(guī)模分子水平數據信息集形成的基因測序和一些高通量科學實驗技術的實用程序數據庫系統資源[18]，是全球最常見的生物信息數據庫管理系統之一。

經過比對，在KEGG信息庫中共有4 310條Unigene信息進行了標注，涉及5條主通道：細胞進程(cellular processes)、環(huán)境信息處理(environmental information processing)、遺傳信息處理(genetic information processing)、新陳代謝(metabolism)、有機合成信息系統(organismal systems)；21條子通路[19]。如圖3所示，參與新陳代謝總體通路 (global and overview maps)的基因有2 132個；參與新陳代謝下的碳水化合物代謝(carbohydrate metabolism)和細胞進程下的運送和分解代謝物(transport and catabolism)通路的基因分別有1 087和701個；最少的是參與新陳代謝下的外源性生物蛋白與新陳代謝物(xenobiotics biodegradation and metabolism)通路的基因，僅有12個。

圖3 KEGG通路分類圖Fig.3 KEGG pathway classification map

2.3 蛋白功能網絡互作圖

為了尋找與麥麩纖維降解過程相關的重要基因，通過分析降解麥麩纖維的白腐菌中不同表達基因，獲得相關的調控交互作用網絡，并利用在線網站STRING[20]分析不同表達基因的蛋白交互作用網絡和Hub基因。使用DIAMOND[21]將所有基因比對至STRING數據庫，并使用與所有已知蛋白的同源關系獲得所有基因之間的互作關聯[22]。基因間的互作關系有評分(150～1 000)，得分越高則互作關系越可靠，默認篩選出評分≥300的網絡關系作圖。

在GO數據庫中，將生物調節(jié)、碳利用、新陳代謝流程、負控制生物學流程、正控制生物學流程、控制生物學流程、抗氧化活性、結合、催化活性條目下FC(Fold Change)前500的差異基因進行蛋白功能網絡互作分析，獲得了一個含有不同基因的66條互作關聯的蛋白相互作用網絡，將互作數據中的Hub基因確認并展示出來，如圖4所示。結果顯示Hub基因分別為10條互作關系的CL16.Contig1_All和10條互作關系的CL1024.Contig1_All。

圖4 GO分類下的蛋白互作網絡圖Fig.4 Protein interaction network diagram with GO classification

如圖5所示，在KEGG條目下，可以從碳水化合物新陳代謝、能源新陳代謝、糖類生物學組成和新陳代謝、輔助因子和維生素代謝通路下FOC值前500的差異基因進行蛋白功能網絡互作，獲得了一個含有113條互作關聯不同基因的蛋白相互作用網絡，并將互作數據中的Hub基因確認并展示出來。結果顯示，Hub基因分別為10條互作關系的CL16.Contig1_All、CL710.Contig10_All、CL3860.Contig3_All和15條互作關系的CL2915.Contig1_All。它們之間存在較強的相互作用關系，可能是白腐真菌降解麥麩纖維的關鍵基因，為進一步探索麥麩纖維降解的菌種培養(yǎng)、菌種優(yōu)化和效率提升提供基礎。

圖5 KEGG分類下的蛋白互作網絡圖Fig. 5 Protein interaction network diagramwith KEGG classification

3 結論

通過對麥麩纖維或葡萄糖為單一碳源培養(yǎng)的白腐菌A.polytricha5.584的基因轉錄組進行高通量測序，獲得了12 096個差異表達基因，GO富集結果顯示差異表達基因主要參與了白腐菌在麥麩纖維等基質環(huán)境生長過程中所必需的能量代謝過程，以及麥麩纖維基質的酶促氧化分解與代謝過程。經GO和KEGG功能注釋與富集解析后通過蛋白網絡互作分析，得到5個與A.polytricha5.584降解麥麩木質素機制相關的關鍵基因。其中有1個注釋到GO數據庫中“氧化還原酶活性”的過氧體水合酶脫氫酶基因(CL16.Contig1_All)表達調高，有2個注釋到Pfam數據庫中“醛脫氫酶家族”的醛脫氫酶基因(CL1024.Contig1_All、CL2915.Contig1_All)，1個注釋到Pfam數據庫中“短鏈脫氫酶”的短鏈脫氫酶基因(CL710.Contig10_All)，1個注釋到KEGG數據庫中“碳代謝”的三羥基萘還原酶基因(CL3860.Contig3_All)。這些基因與白腐菌A.polytricha5.584降解麥麩纖維過程高度相關。其中注釋到GO數據庫中“氧化還原酶活性”的過氧體水合酶脫氫酶基因(CL16.Contig1_All)同時出現在了經篩選GO條目和經篩選KEGG條目后得到的蛋白網絡互作圖中，且其在兩個蛋白網絡互作圖中的熱度之和最高，值得進一步研究關注。