基于轉錄組測序的異常漢遜酵母菌不同發酵時期差異表達基因功能分析

唐朝，張寒玉，王婷，馮光文，錢衛東，蔡長龍，毛培宏*

1(新疆大學 物理科學與技術學院，放射生態與離子束生物技術中心，新疆 烏魯木齊，830046) 2(陜西科技大學 食品與生物工程學院，陜西 西安，710032) 3(西安工業大學，離子束生物工程與生物多樣性研究中心，陜西 西安，710032)

異常漢遜酵母菌具有較強的發酵力和酯化力,是釀造食品香氣成分的主要貢獻者之一[1-3],在其發酵過程中還能積累色氨酸[4]。色氨酸是動物的必需氨基酸之一,其生理作用不可替代[5]。此外,在離子注入介導外源DNA大分子的遺傳轉化研究中,異常漢遜酵母菌常作為外源基因的受體菌[6-7]。

轉錄組對于解釋基因組的功能元件、揭示細胞和組織的分子成分以及對發育和進化的理解具有重要作用。近年來,轉錄組測序(RNA-seq)技術成熟,測序數據具有高度重現性,在技術重復方面幾乎沒有系統差異,同時具有信噪比高、分辨率高、應用范圍廣等優勢,正成為研究基因表達水平、新基因的檢測及注釋的重要手段[8]。

本研究基于RNA-seq,利用生物信息學方法分析異常漢遜酵母菌不同發酵時期的差異表達基因的信息及其功能,以期為人們深入認識異常漢遜酵母菌的基因表達與調控及其色氨酸代謝機制提供線索。

1 材料與方法

1.1 菌株

異常漢遜酵母(Hansenulaanomala),來源于中國普通微生物菌種保藏管理中心(編號2.340,產色氨酸)，由陜西科技大學食品與生物工程學院功能微生物研究室保存。

1.2 培養基

YPD斜面培養基[9](g)：酵母膏10,蛋白胨20,葡萄糖20,瓊脂粉10,蒸餾水1 000 mL。

液體培養基(g)：一水葡萄糖100,NaNO310,酵母膏粉 5，K2HPO4·3H2O 1，MgSO4·7H2O 0.5,蒸餾水1 000 mL,用濃度100 g/L的NaOH調pH至7.0。

1.3 試驗方法

將YPD斜面培養基上培養48 h的新鮮菌種接種至液體培養基中,置搖床230 r/min、28～30 ℃發酵96 h。每間隔24 h取樣(含0 h),于4 ℃、8 000×g離心5 min,充分棄上清,迅速將菌體置于液氮速凍保存。本研究設置3個生物學重復,每個生物學重復的0、24、48、72、96 h的樣品均設置6個技術重復。

分別制備各發酵時間點菌體樣本的mRNA，并應用Illumina進行RNA-seq測序，每個樣本的測序數據量5G。將獲得的每個樣本的轉錄組文庫經過CASAVA堿基識別分析轉化為原始測序序列(raw data)，再經過去除接頭和過濾低質量reads處理，所獲得的RNA序列(clean data)，作為本研究的基本數據。同時，制備異常漢遜酵母菌基因組DNA，并應用PacBio單分子測序技術進行全基因組De novo測序，所獲得的全基因組DNA序列，作為RNA-seq的參考基因組。

根據RNA-seq規程[10],選取HISAT2方法[11]和HTSeq方法[12]分別作為異常漢遜酵母不同發酵時間點差異表達基因的表達量的比對分析方法和FPKM(fragments per kilobase of transcript per million fragments mapped)定量方法。

本研究應用HTSeq方法對差異表達基因的表達量進行FPKM定量和count定量,對表達基因通過BLAST2GO[13]和KAAS(KEGG automatic annotation server)[14]分別比對到GO(gene orthology)數據庫和KEGG(kyoto encyclopedia of genes and genomes)數據庫，利用DESeq2方法[15]對定量結果進行差異表達分析并獲得差異表達基因列表,再通過GOSeq方法[16]采取超幾何分布樣本抽取方法進行差異表達的基因功能的GO富集和KEGG代謝通路(KEGG pathway)富集。

2 結果與分析

2.1 差異表達基因(DEGs)篩選

利用HISAT2方法將RNA序列比對到異常漢遜酵母的基因組后得到sam數據,比對結果均高于90%,表明相關實驗不存在污染。

以FC(fold change,差異表達倍數)和p-adjusted(校正后的概率P值)作為衡量各發酵時間點基因表達的差異程度及顯著性的尺度,以FC>2(|log2FC|>1)、Padj<0.05判定基因是否在樣本間存在差異表達。

通過各發酵時間點的表達基因的相互對比,獲得的DEGs如圖1所示。各發酵時間點與其相鄰的上一個時間點相比,發酵24、48、72、96 h的上調表達的基因數量分別為585、487、154、615,下調表達的基因數量分別為1 112、725、5、245；與0 h相比,發酵48、72、96 h的上調表達基因數量分別為701、639、639,下調表達的基因數量分別為1 315、1 211、1 023；與24 h相比,發酵72 h和96 h的上調表達基因數量分別為581和689,下調表達基因數量分別為587和471(圖1)。

圖1 異常漢遜酵母液體發酵不同時間點間的差異表達基因Fig.1 Differential expressed genes of different culture time at Hansenula anomala liquid culture注：t1/t2表示以t1時間樣本為處理組、t2為對照組。

從各發酵時間點與其相鄰的上一個時間點的比較來看,無論是上調還是下調的差異表達基因,其數量從24～72 h呈下降趨勢,然后增加。發酵96 h時,其上調的差異表達基因數量高于24 h,而下調的差異表達基因數量只有24 h的22.03%。

從各發酵時間點與0 h的比較來看,下調表達的基因數量均超過1 000,其數量在24～48 h上升,然后下降；而上調表達的基因數量在24～48 h呈上升趨勢,48～72 h下降,然后趨于平穩。這些分析結果表明,異常漢遜酵母從新鮮斜面接入液體發酵時(0 h),其基因表達水平已進入一次峰值,在液體發酵96 h時再進入新一輪基因表達峰值。

2.2 差異表達基因功能的GO富集分析

通過FPKM表達矩陣得到異常漢遜酵母菌表達的基因數目為6 086個,與GO數據庫的比對結果顯示,共有3 455個基因對應到3 686種GO功能。

將異常漢遜酵母菌差異表達基因通過GO功能富集(P<0.05)得到的功能與GO數據庫(http://www.geneonthology.org)進行比對，尋找其二級分類,發現其共涉及43個生物學過程(BP,biological process)、13個細胞組分(CC,cellular component)和13個分子功能(MF,molecular function)。在生物學過程的功能分類中,差異表達基因的功能主要涉及細胞過程(cellular process)、代謝過程(metabolic process)、細胞成分(cellular component of organization or biogenesis)；在細胞組分的功能分類中,差異表達基因的功能主要涉及細胞(cell)、細胞成分(cell part)、蛋白質復合物(protein-containing complex)；在分子功能的分類中,差異表達基因功能主要涉及催化活性(catalytic activity)、結合(binding)、轉錄調節活性(transcription regulator activity)。

對異常漢遜酵母菌各個發酵時間點均含有GO富集的差異基因的功能進行分析,結果表明,rDNA沉默(GO: 0000183)、過氧化物酶(GO: 0002215)、ADP磷酸化(GO: 0006119)、蛋白質接合作用(GO: 0045116)、質子運輸(GO: 0015992)、對熱刺激反應(GO: 0009408)、胞質小核糖體亞基(GO: 0022627)、層黏連蛋白-3復合物(GO: 005608)、蛋白質結合(GO: 0030674)均在24 h和48 h下調富集、72 h和96 h上調富集,其變化趨勢均為從0～48 h下降,48～96 h上升(圖2)。

圖2 異常漢遜酵母液體發酵的各時間點差異表達基因的GO功能類別Fig.2 Differential gene GO enrichment in different culture time at Hansenula anomala liquid culture注：t1/t2表示以t1時間樣本為處理組、t2為對照組。

除去發酵時間連續的GO富集結果之后,對剩余GO功能類別進行分析,上調富集的GO功能大部分與細胞分裂有關,其中凋亡DNA片段化(GO: 0006309)、中心體循環(GO: 0007098)、DNA雙鏈斷裂處理(GO: 0000729)、DNA雙鏈斷裂修復(GO: 0006302)、微管解聚負調控(GO: 0007027)、正向調節彈性蛋白(GO: 0051544)、蛋白質泛素化的調節(GO: 0031396)、內膜的組成成分(GO: 0031356)在96 h對比0 h或24 h時均上調富集。將其余富集的GO Term與GO數據庫比對尋找其二級分類發現功能主要分布在催化活性和代謝過程，涉及催化活性相關的有α-1,2-甘露糖基轉移酶(GO:000026)、ADP-L-甘油基-d-甘露-庚糖-6-差向異構酶(GO:0008712)、谷氨酸合酶(GO:0016040)、ATP酶(GO:0016887)、3-羥基-N-甲基-(S)-丙氨酸4-O-甲基轉移酶(GO:0030784)、[核酮糖二磷酸羧化酶]-賴氨酸N-甲基轉移酶(GO:0030785)、質子轉運ATP酶(GO:0046961)、鄰二羥基香豆素7-O-葡糖基轉移酶(GO:0047208)；而和代謝相關的有三環化合物rRNA轉錄物的內切核酸切割(GO:0000449)、ITS1內核裂解性裂解(GO:0000464)、外切核酸修剪(GO:0000465)、產生成熟的LSU-rRNA(GO:0000468)、來自四順反子rRNA轉錄物的SSU-rRNA成熟(GO:0000474)、5S rRNA的成熟(GO:0000481)、來自四順反子rRNA轉錄物的5S rRNA成熟(GO:0000482)、腺嘌呤分解代謝(GO:0006146)、鳥嘌呤分解代謝(GO: 0006147)、dUMP生物合成(GO: 0006226)、蛋白水解(GO:0006508)、谷氨酸生物合成(GO:0006537)、異亮氨酸分解代謝(GO:0006550)、亮氨酸分解代謝(GO:0006552)、賴氨酸分解代謝(GO:0006554)、脯氨酸代謝(GO:0006560)、色氨酸代謝(GO:0006568)、纈氨酸分解代謝(GO:0006574)、脂肪酸生物合成(GO:0006633)、谷胱甘肽代謝(GO:0006749)、脂多糖生物合成(GO:0009103)、生物堿生物合成(GO:0009821)、苯甲酸鹽代謝(GO:0018874)、去甲腎上腺素生物合成(GO:0042421)、去甲腎上腺素分解代謝(GO:0045422)、轉錄調節(GO:0045449)。

2.3 KEGG富集的差異表達基因的功能

將異常漢遜酵母菌表達的6 086個基因與KEGG數據庫進行比對,結果顯示,共有3 275個基因注釋到2 865種功能。

KEGG富集(FDR<0.05)分析結果顯示,DEGs共涉及50條通路,其中與氨基酸代謝相關的通路有11條,包括苯丙氨酸代謝(Phenylalanine metabolism,ko00340)、色氨酸代謝(Tryptophan metabolism,ko00360)、芳香族氨基酸(苯丙氨酸、酪氨酸和色氨酸)的生物合成(Phenylalanine,tyrosine and tryptophan biosynthesis,ko00380)。

異常漢遜酵母菌具有積累色氨酸的代謝特性,對其色氨酸的代謝通路的分析表明,參與色氨酸合成的為色氨酸合成酶(EC: 4.2.1.20)、參與其分解的為色氨酸氨基轉移酶(EC: 2.6.1.27)和吲哚胺2,3-二氧化酶(EC: 1.13.11.52)。由酶的表達圖(圖3)可知,色氨酸的合成速率在0～24 h上升,隨后下降；而分解速率從0～48 h上升,隨后下降。

圖3 異常漢遜酵母液體的各時間點色氨酸代謝相關相關酶Fig.3 The synthesis for metallic of tryptophan at Hansenula anomala liquid culture

對異常漢遜酵母菌各個發酵時間點均含有KEGG富集的差異基因的功能進行分析,結果表明,核糖體(ribosome,ko03010)的代謝趨勢為0 h至 48 h 下降,48 h至96 h上升,這與各個發酵時間點均含有的GO富集的差異基因的功能分析結果一致,說明了異常漢遜酵母菌的合成代謝強度在0 h至48 h下降,48 h至96 h上升。

3 討論

本研究通過Illumina測序平臺首次對異常漢遜酵母5個發酵時間點共15個生物學樣本進行轉錄組測序,共獲得6 086個表達基因。對5個時間點的異常漢遜酵母表達譜的分析顯示,其在生長過程中的基因表達存在顯著差異。與0 h相比,其他發酵時間點均存在大量下調基因,這是由于接入液體培養基的新鮮斜面菌種的代謝旺盛,導致0 h之后的表達基因出現大量下調。另外,本研究在選取差異表達的基因時,差異表達倍數設置為2,使得時間上連續的基因功能富集的種類較少。

差異表達基因的變化趨勢與GO功能富集結果的關聯分析結果表明，下調表達基因數量24～48 h上升,然后下降；而上調表達的基因數量在24～48 h呈上升趨勢,48～72 h下降,然后趨于平穩。0～48 h下調的基因富集功能涉及多種代謝過程和催化活性，而上調基因功能富集的只有蛋白水解和信號識別粒子；48～96 h上調的基因功能也涉及多種代謝過程和催化活性，而其下調基因的功能富集只有谷胱甘肽代謝和脂肪調節酶。可能由于基因表達功能較多，導致即使存在大量差異表達基因，功能富集項也較少，而未富集的基因功能可能與細胞繁殖與衰老相關。異常漢遜酵母菌在細胞整體轉錄組變化和細胞數量變化的共同作用下，導致差異基因表達結果與酵母菌細胞的生長曲線存在差別，其內在機制有待進一步研究。

在異常漢遜酵母菌的差異表達基因功能的GO富集結果中,ADP磷酸化和質子運輸的變化趨勢代表其代謝強度的變化[17],即菌體代謝強度0～48 h下降,48～96 h上升。rDNA沉默的變化趨勢和代謝強度變化趨勢一致,說明菌體在48 h時,細胞已儲存足夠多的核糖體蛋白用于后續蛋白質合成,因此在48 h之后rDNA沉默也呈上升趨勢。除此之外，其他與合成代謝產物前體相關的富集的GO功能類別,其變化趨勢與核糖體蛋白作用速度的變化一致,這表明48 h后菌體代謝產物的合成速度上升。

在96 h上調的微管解聚基因的負調控與酵母菌的有絲分裂相關[18],彈性蛋白基因的正向調節與細胞有絲分裂末期形成細胞膜的組成成分相關[19]；而凋亡DNA片段化的基因表達上調表明96 h時細胞死亡進入一個峰值,這與細胞生長曲線上的細胞進入衰老期相符；DNA雙鏈斷裂處理、修復上調說明大部分細胞進入分裂間期,表明其將進入第二輪基因表達。KEGG富集結果中的核糖體變化趨勢與此一致,可為異常漢遜酵母菌代謝相關研究提供依據。

在異常漢遜酵母菌的差異表達基因功能的KEGG富集的通路中,合成色氨酸的酶為色氨酸合成酶(EC:4.2.1.20)；而分解色氨酸的酶分別為色氨酸氨基轉移酶(EC: 2.6.1.27)和吲哚胺2,3-二氧化酶(EC: 1.13.11.52)。基因表達矩陣的分析結果說明了色氨酸的合成速率在0～72 h下降,隨后上升；而分解速率從0～48 h上升,隨后下降。本研究從基因表達角度印證了此前關于酵母菌發酵色氨酸的相關實驗[20]。

差異表達基因功能的GO和KEGG的富集結果顯示,在本研究條件下,異常漢遜酵母的代謝強度在0～48 h下降,48～96 h上升,并且在96 h時其轉錄、復制和翻譯等基因表達豐度較高,具有較強的代謝活動和遺傳信息處理能力。本研究結果將為異常漢遜酵母菌的分子育種及其代謝調控研究提供理論依據。