基于轉(zhuǎn)錄組測序的異常漢遜酵母菌不同發(fā)酵時期差異表達(dá)基因功能分析

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

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

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

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

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

1 材料與方法

1.1 菌株

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

1.2 培養(yǎng)基

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

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

1.3 試驗方法

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

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

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

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

2 結(jié)果與分析

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

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

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

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

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

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

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

2.2 差異表達(dá)基因功能的GO富集分析

通過FPKM表達(dá)矩陣得到異常漢遜酵母菌表達(dá)的基因數(shù)目為6 086個,與GO數(shù)據(jù)庫的比對結(jié)果顯示,共有3 455個基因?qū)?yīng)到3 686種GO功能。

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

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

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

除去發(fā)酵時間連續(xù)的GO富集結(jié)果之后,對剩余GO功能類別進(jìn)行分析,上調(diào)富集的GO功能大部分與細(xì)胞分裂有關(guān),其中凋亡DNA片段化(GO: 0006309)、中心體循環(huán)(GO: 0007098)、DNA雙鏈斷裂處理(GO: 0000729)、DNA雙鏈斷裂修復(fù)(GO: 0006302)、微管解聚負(fù)調(diào)控(GO: 0007027)、正向調(diào)節(jié)彈性蛋白(GO: 0051544)、蛋白質(zhì)泛素化的調(diào)節(jié)(GO: 0031396)、內(nèi)膜的組成成分(GO: 0031356)在96 h對比0 h或24 h時均上調(diào)富集。將其余富集的GO Term與GO數(shù)據(jù)庫比對尋找其二級分類發(fā)現(xiàn)功能主要分布在催化活性和代謝過程，涉及催化活性相關(guān)的有α-1,2-甘露糖基轉(zhuǎn)移酶(GO:000026)、ADP-L-甘油基-d-甘露-庚糖-6-差向異構(gòu)酶(GO:0008712)、谷氨酸合酶(GO:0016040)、ATP酶(GO:0016887)、3-羥基-N-甲基-(S)-丙氨酸4-O-甲基轉(zhuǎn)移酶(GO:0030784)、[核酮糖二磷酸羧化酶]-賴氨酸N-甲基轉(zhuǎn)移酶(GO:0030785)、質(zhì)子轉(zhuǎn)運ATP酶(GO:0046961)、鄰二羥基香豆素7-O-葡糖基轉(zhuǎn)移酶(GO:0047208)；而和代謝相關(guān)的有三環(huán)化合物rRNA轉(zhuǎn)錄物的內(nèi)切核酸切割(GO:0000449)、ITS1內(nèi)核裂解性裂解(GO:0000464)、外切核酸修剪(GO:0000465)、產(chǎn)生成熟的LSU-rRNA(GO:0000468)、來自四順反子rRNA轉(zhuǎn)錄物的SSU-rRNA成熟(GO:0000474)、5S rRNA的成熟(GO:0000481)、來自四順反子rRNA轉(zhuǎn)錄物的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)、轉(zhuǎn)錄調(diào)節(jié)(GO:0045449)。

2.3 KEGG富集的差異表達(dá)基因的功能

將異常漢遜酵母菌表達(dá)的6 086個基因與KEGG數(shù)據(jù)庫進(jìn)行比對,結(jié)果顯示,共有3 275個基因注釋到2 865種功能。

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

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

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

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

3 討論

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

差異表達(dá)基因的變化趨勢與GO功能富集結(jié)果的關(guān)聯(lián)分析結(jié)果表明，下調(diào)表達(dá)基因數(shù)量24～48 h上升,然后下降；而上調(diào)表達(dá)的基因數(shù)量在24～48 h呈上升趨勢,48～72 h下降,然后趨于平穩(wěn)。0～48 h下調(diào)的基因富集功能涉及多種代謝過程和催化活性，而上調(diào)基因功能富集的只有蛋白水解和信號識別粒子；48～96 h上調(diào)的基因功能也涉及多種代謝過程和催化活性，而其下調(diào)基因的功能富集只有谷胱甘肽代謝和脂肪調(diào)節(jié)酶?？赡苡捎诨虮磉_(dá)功能較多，導(dǎo)致即使存在大量差異表達(dá)基因，功能富集項也較少，而未富集的基因功能可能與細(xì)胞繁殖與衰老相關(guān)。異常漢遜酵母菌在細(xì)胞整體轉(zhuǎn)錄組變化和細(xì)胞數(shù)量變化的共同作用下，導(dǎo)致差異基因表達(dá)結(jié)果與酵母菌細(xì)胞的生長曲線存在差別，其內(nèi)在機(jī)制有待進(jìn)一步研究。

在異常漢遜酵母菌的差異表達(dá)基因功能的GO富集結(jié)果中,ADP磷酸化和質(zhì)子運輸?shù)淖兓厔荽砥浯x強度的變化[17],即菌體代謝強度0～48 h下降,48～96 h上升。rDNA沉默的變化趨勢和代謝強度變化趨勢一致,說明菌體在48 h時,細(xì)胞已儲存足夠多的核糖體蛋白用于后續(xù)蛋白質(zhì)合成,因此在48 h之后rDNA沉默也呈上升趨勢。除此之外，其他與合成代謝產(chǎn)物前體相關(guān)的富集的GO功能類別,其變化趨勢與核糖體蛋白作用速度的變化一致,這表明48 h后菌體代謝產(chǎn)物的合成速度上升。

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

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

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