西藏黑木耳高產菌種不同生長發育時期轉錄組測序分析

薛 嬌, 于 洋, 張彥龍, 雷 虹, 曾偉民

(1.黑龍江大學 農業微生物技術教育部工程研究中心, 哈爾濱 150500; 2.黑龍江大學 生命科學學院, 哈爾濱 150080)

0 引 言

黑木耳屬于真菌界,擔子菌亞門,層菌綱,木耳目,木耳科,木耳屬[1]。黑木耳口感脆嫩,味道鮮美可口,我國對其的食用已有上千年的歷史。黑木耳子實體營養豐富,其多糖具有降血脂、抗血栓和抗腫瘤等作用[2],黑色素能夠明顯改善四氯化碳造成的急性肝損傷,多酚類和黃酮類物質具有抗氧化和清除自由基等功能。黑木耳在我國種植廣泛,在云南、貴州、新疆、華北、內蒙古、西藏及東北地區均有分布。西藏黑木耳由于生長位置海拔較高和紫外線輻射較強的獨特地理條件而具有抗逆性強、耐低溫、生長周期短和營養價值高的優點[3]。

轉錄組在廣義上是指整個細胞在某一特定條件或狀態下所轉錄出的編碼RNA和非編碼RNA的總和,轉錄組學是在整體水平上研究生物個體所轉錄出RNA的差異情況[4]。轉錄組測序是通過將RNA進行測序,獲得表達差異基因,注釋到各個數據庫從而進行分析。目前,用于轉錄組研究的主要有SAGE基因表達序列分析技術、大規模平行測序技術、轉錄組測序技術和基因芯片技術[5]。轉錄組測序技術并不局限于對已知轉錄本信息的樣品進行檢測,還可以利用高通量測序對無參考序列信息的樣品進行檢測。近年來隨著測序成本的降低和其靈敏度高、檢測范圍廣和重復性好的特點被廣泛應用到各個領域。王思雨等對成菇期和變色期雞腿菇子實體進行轉錄組測序分析,發現變色期DNA復制、轉錄和翻譯能力下降,這可能與雞腿菇自溶現象有關[6]。劉璐等通過對羊肚菌正常菌絲和退化菌絲進行轉錄組測序,篩選出一些差異表達基因,注釋后得出羊肚菌菌絲退化與代謝過程有關,為后續解決菌絲退化問題提供參考[7]。趙震宇等通過對不同濃度單寧處理過的玉木耳菌絲進行轉錄組測序,得到2 012個差異基因,為進一步研究單寧對玉木耳菌絲的影響機制提供參考[8]。本試驗采用課題組分離選育出的西藏黑木耳高產菌種,對其子實體3個不同發育階段轉錄組進行分析研究。對黑木耳子實體的unigene使用生物信息學分析工具,進行功能注釋和分析。該研究為下一步篩選和挖掘有關黑木耳子實體生長發育的功能基因提供重要數據參考,以及為黑木耳品種改良提供了科學依據。

1 材料與方法

1.1 材料與儀器

1.1.1供試菌株

西藏黑木耳高產菌種(專利號:ZL 2013 1 0278146.8)由本實驗室前期對比選育分離得到,命名為西藏6號[3]。

1.1.2樣品的采集

黑木耳子實體采摘自黑龍江省勃利縣黑木耳種植基地,在黑木耳子實體生長的三個發育階段(幼耳期、旺盛期和成熟期),采取無雜菌、無污染、品質優良的黑木耳置于放有冰袋的泡沫盒中,立即運至實驗室,轉移至凍存管后,使用液氮將樣品速凍0.5 h,置于-80 ℃超低溫冰箱中用于RNA提取,委托上海美吉生物醫藥科技有限公司測序。

1.1.3儀器

DU530紫外可見分光光度計(美國Beckman公司);DYY-12電泳儀(北京六一生物科技有限公司);DYCZ-24DN水平電泳槽(北京六一生物科技有限公司);Agilent 2100生物芯片分析系統(美國Agilent科技有限公司);Illumina NovaSeq 6000測序系統(美國Illumina公司);Elite超低溫冰箱(美國Revco公司)。

1.2 試驗方法

1.2.1總RNA樣品提取

分別在西藏6號黑木耳幼耳期、旺盛期和成熟期取樣,每個樣品設置3個平行實驗,共9個樣品。每個樣品取50 mg,分別加入250 μL NucleoZOL裂解細胞組織,再加入200 μL RNase-free的水沉淀污染物,分別用250 μL異丙醇和250 μL 75%乙醇沉淀和洗滌RNA。用NanoDrop[9]檢測RNA吸收峰是否正常,利用瓊脂糖凝膠電泳對提取到的待測樣品RNA進行檢測,使用Agilent 2100[10]精確檢測RNA的RIN值、28S/18S和5S峰。所有指標均需符合測序的質量要求。

1.2.2RNA測序

提取后的RNA為總RNA,需要對其中的mRNA進行富集。將堿基末端帶有T的磁珠與測序RNA 3′末端的ployA結合,加入buffer將富集后的mRNA打成短片段。以mRNA為模板,加入六堿基隨機引物和逆轉錄酶,合成單鏈cDNA,再加入dDNTPs、DNA聚合酶和緩沖液合成第二條鏈,最后加A尾,PCR擴增15個循環,對文庫進行測序[11]。RNA測序服務由上海美吉生物醫藥科技有限公司提供。

1.2.3轉錄組測序數據處理

測量后的數據中包含測序接頭序列、不確定堿基信息N率較高的序列、一些低質量的讀段片段和一些影響分析的長度過短的片段,這些因素都會影響后續分析結果的準確性。對于需要去除接頭序列,并且序列中含有質量小于10的堿基的reads,整條序列都需要剔除。為了得到高質量的有效讀數片段,需要剔除不確定堿基信息N率大于10%的片段和修剪后片段長度小于30 bp的序列。本研究無參考基因組,利用Trinity[12](https://github.com/trinityrnaseq/trinityrnaseq)將所有質量控制后的測序數據進行組裝,組裝后利用TransRate[13](http://hibberdlab.com/transrate/)和CD-HIT[14](http://weizhongli-lab.org/cd-hit/)對組裝結果進行優化過濾,利用BUSCO[15](Benchmarking universal single-copy orthologs, http://busco.ezlab.org)對組裝結果進行評估,評估基因組或轉錄組的組裝完整性。取測序數據最長的轉錄本作為unigene進行后續分析。

1.2.4差異表達基因篩選

使用DESeq2軟件[16]對Unigene進行分析,以篩選差異表達基因,FDR<0.05的基因為顯著差異基因。

1.2.5Unigene功能注釋

將所有的Unigene與六大數據庫(GO、KEGG、COG、Nr、Swiss-Prot和Pfam)進行比對,獲得Unigene和Transcript的注釋情況,并對各注釋信息進行統計分析。

1.2.6GO和KEGG功能富集分析

利用Goatools軟件[17](https://github.com/tanghaibao/GOatools)和R語言編寫腳本對基因集中的轉錄本進行GO和KEGG富集分析。GO和KEGG通路富集均默認,P值(Pvalue_corrected)<0.05時,GO和KEGG通路存在顯著富集情況。

2 結果與分析

2.1 測序數據質量控制

黑木耳子實體轉錄組測序數據如表1所示。由表1可以看出,從西藏黑木耳的9個樣品中共獲得51 482 040～77 672 734條統計過濾后測序數據,Clean bases為7 685 784 278～11 565 940 347 G。樣品測序的Q20含量超過98%,GC含量為60.54%～61.48%,Q30堿基含量均不小于94.84%。根據結果對比統計得出,各樣品的Reads測序錯誤率為0.023 9～0.024 3。熊宇晴等對無患子葉片進行高通量測序后,得出數據質量情況,GC占總堿基數的42.73%,Q20和Q30含量分別為97.79%和94.00%,符合下一步的組裝要求[18]。從西藏6號黑木耳質控數據結果得出,樣品測序結果較好,符合轉錄組進行組裝和注釋的要求。

表1 質控數據表

2.2 轉錄組從頭組裝結果

對Transcript有效讀取片段進行組裝拼接,去掉其中冗雜多余的片段,共生成28 817個Transcript(表2),組裝得到的所有Transcript堿基為36 006 442個,平均長度為1 249.49 bp。Unigene共有13 765個,將所有的Unigene組裝得到的堿基為15 530 526個,平均長度為1 128.3 bp。西藏6號黑木耳Unigene的N50數值為1 954。在物種差異、拼接算法和測序平臺都正確的前提下,N50數值能夠作為評估序列測度的重要參數,反映拼接所得序列的完整程度。E90N50長度分別為2 891和2 640。所有樣本質控后的有效數據合并后與組裝Unigene/Transcript進行比較,獲得的mapped率分別為59.606%和66.171%?？倝A基數量中的58.19%和58.91%是GC堿基。李方東利用BUSCO對茶樹中直系同源單拷貝基因完整性進行考察,以評價轉錄本組裝的準確性和完整性,相似性評分越高,組裝結果越好[19]。用BUSCO評估Unigene和Transcript,得分分別為83.5%和90.0%,表明組裝結果完整性和準確性較高。

表2 優化組裝結果評估表

如表3所示,過濾后測序數據為25 741 020～38 836 367條,可以對比到組裝Transcript上的過濾后測序數據為17 330 951～27 367 399條,可以追蹤到組裝Transcript上的過濾后測序數據所占百分比為67.03%～70.47%。

表3 測序數據與組裝結果對比

2.3 差異表達基因分析

表4為差異基因數目統計表,從獲得的差異表達基因集的基因數目統計得出,西藏黑木耳子實體旺盛期和成熟期差異基因為1 225個,上調基因為661個,下調基因為564個。幼耳期和旺盛期差異基因為2 096個,比旺盛期和成熟期的差異基因多,同時上調基因也比下調基因表達得更多,上調基因為1 225個,下調基因為871個。

表4 差異基因數目統計表

圖1 差異基因數目統計圖(a)和維恩圖(b)

圖1為差異基因數目的統計圖和信息圖。維恩圖顯示,西藏6號黑木耳幼耳期、旺盛期和成熟期3個不同發育時期的共同差異基因有13 281個。其中,成熟期和旺盛期的差異基因有14 335個,幼耳期和成熟期的差異基因有13 717個,幼耳期和旺盛期的差異基因有14 093個。

2.4 功能注釋及統計

將組裝后的黑木耳子實體轉錄組序列與六大數據庫進行對比,所有的Unigene和Transcript數量分別為27 529個和57 634個,表達的Unigene和Transcript數量分別為27 444個和57 406個,注釋到數據庫的表達Unigene和Transcript數量分別為18 258個和40 469個。在六大數據庫上表達的Unigene分別注釋了1 1967個(GO)、7 658個(KEGG)、5 252個(COG)、17 132個(Nr)、9 240個(Swiss-Prot)和11 328個(Pfam),表達的Transcript分別注釋到27 232個(GO)、15 425個(KEGG)、9 090個(COG)、38 968個(NR)、19 555個(Swiss-Prot)和24 214個(Pfam)。

表5 轉錄本功能注釋統計表

2.5 GO功能富集分析

對西藏6號黑木耳子實體幼耳期和旺盛期差異表達基因進行GO功能富集分析,結果如圖2所示?；蚣械幕蚩梢苑殖葿P(Biological Process, 生物過程)、MF(Molecular Function, 分子功能)及CC(Cellular Component, 細胞組分)三個類別。在“生物過程”類別中,參與類異戊二烯生物合成過程(GO:0008299)的基因有6個,類異戊二烯代謝過程(GO:0006720)的基因有6個,谷氨酸脫羧成琥珀酸過程(GO:0006540)基因有3個,GO term顯著富集,表明黑木耳子實體在生長發育和代謝調節等方面活動頻繁。在“細胞組分”類別中,差異表達基因主要富集在膜的組成部分(GO:0016021)和膜的固有組分(GO:0031224),分別有415個。在“分子功能”類別中,最為富集的Unigene具有催化活性和結合活性。表明黑木耳子實體中大量的Unigene可以翻譯為具有催化活性的酶和具有結合功能的蛋白。

圖2 黑木耳子實體幼耳期和旺盛期差異表達基因GO功能富集分析圖

將黑木耳子實體旺盛期和成熟期的差異表達基因進行GO功能富集分析,結果如圖3所示。在“細胞組分”類別中,參與膜的固有組分(GO:0031224)和膜的組成部分(GO:0016021)的差異表達基因分別有168個。在“分子功能”類別中,參與金屬末端肽酶活性(GO:0004222)的基因有12個,參與肽鏈內切酶活性(GO:0004175)的基因有23個,表明黑木耳子實體中具有催化肽類活性的Unigene最為富集。從GO功能富集分析結果可以看出,在生物學過程類別中,差異表達基因主要富集在異類戊二烯的合成和代謝過程,說明差異表達基因在“生物學過程”主要參與細胞內進行的各種異類戊二烯的合成和代謝活動。在“細胞組分”類別中,富集到膜的固有組分和膜的組成部分的差異表達基因最多,表明細胞膜成分隨著西藏6號黑木耳不斷生長而發生改變。在“分子功能”類別中,在催化活性通路的差異表達基因最多。曾倩倩等對小麥葉片差異表達基因進行GO富集分析,檢測到的1 207個DEGs富集到3個類別的34個條目下,主要參與催化、代謝和水解酶活性等,表明小麥受干旱環境的脅迫作用是多方面的[19]。

圖3 黑木耳子實體旺盛期和成熟期差異表達基因GO功能富集分析圖

2.6 KEGG功能富集分析

將西藏6號子實體幼耳期和旺盛期差異表達基因進行KEGG功能富集分析,結果如圖4所示。富集到的Unigene有6個類別,20個亞類。在子類別中,p<0.05的有8個亞類。Unigene最為富集的途徑是碳水化合物代謝(16%)、氨基酸代謝(15%)、脂質代謝(12%)、輔助因子和維生素代謝(9%)、其他氨基酸代謝(7%)以及運輸和分解代謝(7%),說明西藏6號黑木耳子實體生長發育過程中碳代謝和氮代謝較為旺盛,為黑木耳子實體生長發育提供營養與能量。

圖4 黑木耳子實體幼耳期和旺盛期差異表達基因KEGG功能富集分析圖

曾倩倩等對小麥葉片DEGs的KEGG通路富集發現,存在淀粉和蔗糖代謝、葉綠素代謝、光合作用、苯丙烷生物合成、甘油磷脂代謝、氨基酸代謝及光合生物中的碳固定等通路。表明植物受干旱環境的脅迫,光合作用-天線蛋白、葉綠體組分、氨基酸代謝和碳代謝等通路受到影響,光合作用和物質積累受到抑制,這可能與小麥葉片相對含水量降低及MDA含量增加有關[20]。盧園萍等通過對巴氏蘑菇子實體3個階段的差異表達基因進行KEGG富集分析發現,原基期上調DEGs主要富集在核糖體蛋白和DNA復制,表明原基期細胞代謝旺盛,其中核糖體蛋白基因上調為后期蛋白質合成提供重要場所;采收期和開傘期差異表達基因主要富集在碳水化合物代謝、脂肪酸降解和氨基酸代謝等途徑,為巴氏蘑菇子實體的生長發育與成熟提供營養和能量[21]。吳曉梅等對雙孢蘑菇子實體3個發育時期進行KEGG富集分析,發現DEGs參與了氨基酸代謝、碳水化合物代謝、核苷酸代謝、脂類代謝和能量代謝這五大代謝通路,其中差異基因主要富集在氨基酸代謝通路中,氨基酸合成相關的多數基因上調表達,表明雙孢蘑菇子實體發育形成需要一系列代謝反應協同調控,氨基酸代謝相關基因可能在雙孢蘑菇子實體發育過程中起重要作用[22]。

將西藏6號黑木耳子實體旺盛期和成熟期差異表達基因進行KEGG功能富集分析,結果如圖5所示。富集到的Unigene主要有6個類別,19個亞類。有5個亞類在子類別中p<0.05,Unigene最為富集的途徑是氨基酸代謝(19%)、碳水化合物代謝(13%)、脂質代謝(12%)、聚糖的生物合成和代謝(9%)、運輸和分解代謝(7%)、其他氨基酸代謝(6%)及能量代謝(4%)。表明黑木耳生長發育過程中代謝過程極其豐富,氨基酸代謝在黑木耳子實體生長中提供能量和營養。與吳曉梅等對雙孢蘑菇子實體不同發育階段轉錄組分析結果[21]類似,氨基酸代謝和碳水化合物代謝也有大量的差異基因富集。

圖5 黑木耳子實體旺盛期和成熟期差異表達基因KEGG功能富集分析圖

3 結 論

對黑木耳子實體發育過程中3個不同階段樣品的轉錄組進行測序,并對轉錄組進行注釋以及GO和KEGG功能富集分析。對黑木耳子實體9個樣品的測序共獲得 81.93 Gb Clean data,組裝后得到Unigene為13 765個,Transcript為28 817個,平均長度為1 128.30 bp,N50長度為1 954 bp。以上結果表明,黑木耳子實體轉錄組測序和組裝質量都較高。采用生物信息學進行分析,將組裝后黑木耳子實體轉錄組序列與六大數據庫進行對比,注釋到數據庫的Unigene共有18 258個,占全部Unigene的66.53%。注釋到數據庫的Transcript共有40 469個,占全部Transcript的70.5%。

從GO功能富集分析的結果可以看出,在“生物學過程”類別中,差異表達基因主要富集在異類戊二烯的合成和代謝過程,說明差異表達基因在生物學過程主要參與細胞內進行的各種異類戊二烯的合成和代謝活動。在“細胞組分”類別中,富集到膜的固有組分和膜的組成部分的差異表達基因最多,表明隨著西藏6號黑木耳不斷生長,細胞膜成分發生改變,細胞膜能夠調節物質運輸,維持細胞內環境的穩定。在“分子功能”類別中,在催化活性通路的差異表達基因最多。KEGG富集分析有利于了解基因的生物學功能,本研究對西藏6號黑木耳子實體3個不同時期發育階段差異基因進行KEGG分析,篩選出多條與黑木耳子實體不同生長發育階段有關的代謝途徑,包括氨基酸代謝、碳水化合物代謝、脂質代謝和聚糖的生物合成和代謝等。表明黑木耳生長發育過程中代謝過程極其豐富,氨基酸代謝在黑木耳子實體生長中提供營養和能量。本研究結果為下一步篩選和挖掘有關西藏黑木耳子實體生長發育的功能基因提供重要數據參考,為黑木耳品種改良提供了科學依據。