球孢白僵菌分生孢子耐低濕萌發(fā)轉(zhuǎn)錄組相關(guān)差異基因分析

王子夜，張海劍，路曉月，王瑞東，王志剛，閻愛華

(1河北農(nóng)業(yè)大學(xué) 林學(xué)院，河北 保定 071000；2河北省林木種質(zhì)資源和森林保護重點實驗室，河北 保定 071000；3河北省農(nóng)林科學(xué)院 植物保護研究所，河北 保定 071000；4河北省農(nóng)業(yè)有害生物綜合防治工程技術(shù)研究中心，河北 保定 071000；5農(nóng)業(yè)部華北北部作物有害生物綜合治理重點實驗室，河北 保定 071000; 6河北省城市森林健康技術(shù)創(chuàng)新中心，河北 保定 071000)

球孢白僵菌(Beauveriabassiana)是重要的昆蟲病原真菌，其宿主范圍廣、特異性高、安全性好，在農(nóng)林害蟲生物防治中具有良好的應(yīng)用前景。作為生防制劑的主要成分，分生孢子能否在不良環(huán)境脅迫下大量宿存、成功萌發(fā)、高效侵入是決定球孢白僵菌(B.bassiana)防治效果的關(guān)鍵[1-3]?？諝庀鄬穸仁抢ハx病原真菌應(yīng)用的限制性環(huán)境因子之一，RH90%～100%是分生孢子萌發(fā)最適濕度，很少昆蟲病原真菌菌株能在低于RH70%的濕度下萌發(fā)[4-5]。探明不同濕度對白僵菌萌發(fā)的作用機理，對于拓展白僵菌的菌種資源，增強對害蟲的防治效果，擴大生防制劑的應(yīng)用范圍具有重要意義。

對白僵菌-昆蟲-環(huán)境三者的相互作用機制已經(jīng)有大量研究，隨著分子生物學(xué)、高通量測序技術(shù)等廣泛應(yīng)用，海藻糖酶(Beauveriabassiananeutral trehalase，BbNTH1)基因、特異性生長抑制蛋白(Growth arrest specific protein，GAS1)基因、熱激蛋白70(Heat shock protein，HSP70)基因、絲裂原活化蛋白激酶(Mitogen-activated protein kinases，MAPK)同源基因slt2和hog1等與白僵菌附著、侵入釘形成、定殖和產(chǎn)孢相關(guān)的一系列信號傳導(dǎo)、抗逆和形態(tài)建成等功能相關(guān)的基因被克隆和鑒定[6-11]。作為真菌生長發(fā)育的起始階段，與萌發(fā)相關(guān)基因和調(diào)控機制的報道還很少[12-13]。本研究通過構(gòu)建不同空氣相對濕度下白僵菌菌株萌發(fā)差異表達基因(Differently expressed genes,DEGs)轉(zhuǎn)錄組文庫，通過基因通路富集和功能注釋，挖掘可能參與白僵菌低濕萌發(fā)的差異基因及其功能和代謝過程，為球孢白僵菌的進一步廣泛應(yīng)用提供理論基礎(chǔ)，也為了解其他蟲生生防菌萌發(fā)功能基因研究提供了重要線索。

1 材料與方法

1.1 供試菌株

供試球孢白僵菌(B.bassiana)菌株WA采自河北省保定市，由河北農(nóng)業(yè)大學(xué)林木病蟲害無公害防控實驗室分離自光肩星天牛(Anoplophoraglabripennis)僵蟲并進行保存。在薩式培養(yǎng)基(葡萄糖40 g，蛋白胨10 g，酵母膏10 g，瓊脂20 g/L )上28 ℃活化培養(yǎng)15 d，收集分生孢子。

1.2 不同相對空氣濕度分生孢子萌發(fā)率測定與文庫樣品采集

采用小容器空氣濕度調(diào)節(jié)法測定空氣相對濕度對白僵菌孢子萌發(fā)的影響，在2個口徑210 mm玻璃密閉干燥器中分別注入滅菌飽和NaCl溶液和滅菌水，構(gòu)建RH75%和RH100%的小空間氣候室。載玻片上覆蓋玻璃紙，將20 μL分生孢子懸浮液(約1.03×107個/mL)滴在玻璃紙上，迅速風干后置于密閉小空間氣候室，25 ℃恒溫培養(yǎng)24 h后取下玻璃紙，鏡檢記錄觀察分生孢子萌發(fā)率和芽管長度。同時大量收集不同處理分生孢子，液氮速凍，在零下80 ℃下保存，用于文庫構(gòu)建。樣品分別命名為WA_75和WA_100，以未經(jīng)處理且未萌發(fā)分生孢子為對照，命名為WA_S。

1.3 cDNA文庫構(gòu)建及轉(zhuǎn)錄組序列測定

cDNA文庫構(gòu)建由諾和致源生物信息技術(shù)有限公司完成，獲得的序列在 Illumina HiSeqTM2000測序儀進行轉(zhuǎn)錄組測序。測序所得的原始片段(Raw reads)經(jīng)過過濾，得到純凈片段(Clean reads)。純凈片段經(jīng)Trinity軟件進行Denovo 從頭合成拼接，組裝去冗余后得到Unigenes進行數(shù)據(jù)分析。

1.4 差異表達Unigenes數(shù)據(jù)注釋與分析

差異Unigenes的表達豐度用FPKM值表示[q value<0.05 & |log2(fold change)|>1]，序列信息分別在Nr數(shù)據(jù)庫、NT數(shù)據(jù)庫,Swiss-Prot數(shù)據(jù)庫、KEGG數(shù)據(jù)庫、COG數(shù)據(jù)庫和GO數(shù)據(jù)庫中比對和檢索，獲得與Unigenes對應(yīng)的蛋白功能注釋及功能分類統(tǒng)計。差異基因表達熱圖采用歐易云在線軟件(https://cloud.oebiotech.com)進行。

1.5 差異表達基因?qū)崟r熒光定量檢測

為驗證不同濕度白僵菌萌發(fā)轉(zhuǎn)錄組文庫構(gòu)建質(zhì)量，以白僵菌18SrRNA基因為內(nèi)參基因，隨機挑選了文庫中5個DEGs序列，用Probegene設(shè)計RT-qPCR引物(表1)，驗證差異基因表達趨勢[14]。按照材料與方法1.2收集萌發(fā)的WA孢子，經(jīng)過RNA提取并反轉(zhuǎn)錄各樣品的cDNA，每處理3個重復(fù)，按2-ΔΔCt相對定量法計算基因相對表達量。

表1 RT-qPCR引物Table 1 RT-qPCR primers

1.6 數(shù)據(jù)統(tǒng)計

采用Excel 2003進行數(shù)據(jù)處理，采用SPSS 17.0進行單因素ANOVA檢驗分析，滿足P<0.05,再進行方差分析和多重比較，并用Origin 2021軟件作圖。

2 結(jié)果與分析

2.1 不同濕度下菌株萌發(fā)情況測定

不同濕度下菌株萌發(fā)情況測定，見圖1。

圖1 不同空氣相對濕度下分生孢子萌發(fā)率和萌發(fā)芽管平均長度

如圖1所示，WA菌株具有較強的耐低濕萌發(fā)能力，在RH100%和RH75%時，菌株WA的分生孢子萌發(fā)率和萌發(fā)芽管長度沒有顯著差異(P<0.05)。

2.2 轉(zhuǎn)錄組測序數(shù)據(jù)質(zhì)量分析與注釋

各樣品數(shù)據(jù)質(zhì)量分析，見表2。

表2 樣品數(shù)據(jù)質(zhì)量分析Table 2 Summery of sample sequence data output quality

如表2所示，WA_S、WA_100和WA_75轉(zhuǎn)錄組文庫數(shù)據(jù)質(zhì)量評估去除不可靠數(shù)據(jù)后，3個樣品各獲得純凈片段為43 016 924條，44 723 918條和43 537 642條，Q30達到92%以上，錯配率控制在1%以下，堿基A、T、C和G之間沒有明顯的偏差。組裝的轉(zhuǎn)錄本和Unigene平均長度分別為547 bp和445 bp，測序和Unigene的拼接質(zhì)量較好。序列分別在Nr、GO、KOG、Swiss-Prot、KEGG、Nt 6個數(shù)據(jù)庫數(shù)據(jù)庫中注釋，共注釋到89 903個Unigenes，其中27 786個Unigenes可以注釋到至少1個文庫中，數(shù)據(jù)庫成功匹配率達到50.43%。

2.3 不同空氣相對濕度下不同菌株DEGs分析

差異基因韋恩圖分析，見圖2。

圖2 差異基因韋恩圖分析Figure 2 Representing and Venn diagram of DEGs in response to different humidity

如圖2所示，與未萌發(fā)孢子相比，不同空氣相對濕度下WA菌株DEGs數(shù)量如圖，WA_75與WA_S共有差異表達基因767個，其中上調(diào)342個，下調(diào)425個，WA_100與WA_S共有差異表達基因822個，上調(diào)400個，下調(diào)422個。393個DEGs在WA_75與WA_S和WA_100與WA_S中共同表達，其中上調(diào)114個，下調(diào)279個。

2.4 分生孢子萌發(fā)差異表達基因GO注釋結(jié)果

差異最顯著30條GO功能注釋通路，見圖3。

如圖3所示，不同濕度下孢子萌發(fā)差異表達基因在GO文庫注釋到細胞組分(Cellular component，CC)，生物進程(Biological process，BP)，分子功能(Molecular function，MF)3個大類44個亞類。與未萌發(fā)的孢子相比，不同濕度條件下萌發(fā)孢子的DEGs在氧化還原進程(Oxidation-reduction process)、真菌細胞壁組分(Fungal-type cell wall)、核糖體合成(Ribosome biogenesis)等亞類中顯著富集(圖3-a、圖3-b)，而與WA_100相比，WA_75差異基因功能差異顯著(圖3-c)，主要集中在BP大類中，在單羧酸代謝進程(Monocarboxylic acid metabolic process，65)，分解代謝進程(Catabolic process，72)，有機物分解代謝進程(Organic substance catabolic process，65)，有機酸代謝進程(Organic acid metabolic process，107)，酮酸代謝進程(Oxoacid metabolic process，106)等方面差異顯著，MF大類中異檸檬酸裂解酶活性(Isocitrate lyase activity，7)差異顯著(P<0.05)。

注：A是WA_75與WA_S；B是WA_100與WA_S；C是WA_75與WA_100。

2.5 分生孢子萌發(fā)差異表達基因KEGG通路富集分析

WA_75與WA_100差異DEGs富集最多的20條KEGG通路，見圖4。

圖4 WA_75與WA_100差異DEGs富集最多的20條KEGG通路Figure 4 WA_75vsWA_100 differential DEGs enriched the 20 most KEGG pathways

如圖4所示，利用超幾何檢驗獲得差異基因顯著富集的Pathway，確定不同基因間互作參與的生化代謝途徑和信號轉(zhuǎn)導(dǎo)途徑。文庫中共有23 549條DEGs注釋到KEGG代謝通路中，分為細胞過程(Cellular processes)，環(huán)境信息處理(Environmental information processing)，遺傳信息處理(Genetic information processing)，代謝(Metabolism)，有機系統(tǒng)(Organismal systems) 5大類20個亞類297條代謝途徑。與WA_100相比，WA_75中365條差異DEGs富集，富集顯著的前20條通路，在糖酵解與糖代謝(Glycolysis / Gluconeogenesis，31)，賴氨酸降解(Lysine degradation，12)，果糖與甘露糖代謝(Fructose and mannose metabolism，10)，酪氨酸代謝(Tyrosine metabolism，11)，甘油酯代謝(Glycerolipid metabolism，9)等方面數(shù)量最多(P<0.05)。

2.6 分生孢子萌發(fā)差異基因表達模式分析

2.6.1 呼吸作用相關(guān)基因 呼吸作用相關(guān)基因表達模式，見圖5。

圖5 呼吸作用相關(guān)基因表達模式Figure 5 Expression patterns of respiration related DEGs

如圖5所示，在有氧條件下，白僵菌萌發(fā)時呼吸作用劇烈，聚類分析結(jié)果表明，氧化磷酸化相關(guān)基因在WA_75和WA_100中表達模式相似，但ATP 合成酶(ATP synthase)、NADH脫氫酶(NADH dehydrogenase，NADH-CoQ)、琥珀酸脫氫酶(Succinate dehydrogenase，SDH)、NADH-泛醌氧化還原酶(NADH-ubiquinone oxidoreductase，NDH)和細胞色素c氧化酶(Cytochromecoxidase，COX)相關(guān)序列在WA_75中上調(diào)。

2.6.2 信號轉(zhuǎn)導(dǎo)蛋白激酶基因表達模式分析 蛋白激酶基因表達模式分析，見圖6。

圖6 蛋白激酶基因表達模式分析Figure 6 Expression patterns of protein kinase

如圖6所示，蛋白激酶催化蛋白磷酸化，是多種信號轉(zhuǎn)導(dǎo)途徑中的重要環(huán)節(jié)。聚類熱圖分析表明，低濕萌發(fā)過程中，WA_75文庫中大量絲氨酸和蘇氨酸蛋白激酶基因表達上調(diào)，包括BR激酶(BR Serine/threonine-protein kinase,BRSK)、DBF20激酶(Serine/threonine-protein kinase，DBF20)、GCN2激酶(Serine/threonine-protein kinase，GCN2)、鈣依賴性蛋白激酶(Calcium-dependent protein kinases,CDPKs)、生發(fā)中心激酶(Germinal center kinase，GCKs)、酪蛋白激酶I(Casein kinase I，CKI)；磷酸肌醇3激酶(Phosphoinositide 3-kinase，PK3)、二酰甘油激酶(Diacylglycerol kinase)等；參與細胞周期調(diào)控的MPS1激酶(Serine/threonine-protein kinase，MPS1)、plo1激酶、細胞周期蛋白依賴性激酶(Cyclin-dependent kinase，CDKs)、雷帕霉素靶蛋白(TOR)激酶等。

2.7 差異表達基因的實時熒光定量驗證

DEGs表達量的差異變化，見圖7。

圖7 DEGs表達量的差異變化Figure 7 Changes of DEGs expression

如圖7所示，利用qRT-PCR技術(shù)驗證隨機挑選的5個DEGs表達量的差異，各基因FPKM值和相對表達量雖然在變化程度上有差異，但不同基因均呈上升趨勢。細胞壁脅迫響應(yīng)組份蛋白(WSC domain-containing protein)、假定泛素蛋白連接酶(Putative RING-H2 finger protein)和Zn (II)2Cys6鋅簇蛋白[Zn (II)2Cys6 transcription factor]的FPKM值變化程度較大，鋅指和BTB結(jié)構(gòu)域蛋白(Zinc finger and BTB domain-containing protein)和假定泛素蛋白連接酶(Putative RING-H2 finger protein)FPKM值變化程度較??；假定外切聚磷酸酶相對表達量變化程度最大，其他4個基因變化差異較小。

3 討論與結(jié)論

盡管溫度、光照、濕度都被認為是影響孢子萌發(fā)的關(guān)鍵環(huán)境因素，但實踐應(yīng)用證明，濕度是限制白僵菌使用范圍和效果的決定性因子[15]。本研究中球孢白僵菌孢子在不同濕度下萌發(fā)時都有大量差異基因表達。分生孢子生理活性狀態(tài)由自身的養(yǎng)分狀況決定，完全成熟的分生孢子在靜止狀態(tài)時仍存在轉(zhuǎn)錄活性，能夠感受環(huán)境條件以進行主動的萌發(fā)決策[16]。利用轉(zhuǎn)錄組學(xué)技術(shù)探討球孢白僵菌孢子在不同空氣濕度條件下萌發(fā)差異基因表達模式，不僅可以揭示白僵菌相應(yīng)水分脅迫的調(diào)控機制，而且為蟲生真菌進一步擴大應(yīng)用研究提供了良好的切入點。

侵染早期基因表達模式的差異是導(dǎo)致球孢白僵菌不同菌株致病性差異的原因，Wang等研究表明毒力與白僵菌的粘附、角質(zhì)層降解相關(guān)功能基因的高表達密切相關(guān)[17]。萌發(fā)出芽管是真菌生長發(fā)育的開端，是菌落形態(tài)建成中第一個重大變化，在本研究中，不同濕度下孢子萌發(fā)時氧化還原進程，真菌細胞壁組份合成和蛋白質(zhì)翻譯進程被顯著富集，為白僵菌快速萌發(fā)提供物質(zhì)基礎(chǔ)。同時，比WA_100，不利的水分條件WA_75抑制了有機氮化物、一元羧酸和氨基酸代謝進程等，但增強了對自身營養(yǎng)和資源的利用效率。單羧酸代謝進程、有機物分解代謝進程、有機酸代謝進程、酮酸代謝進程等相關(guān)生物學(xué)進程基因表達活躍，在不利條件下維持正常萌發(fā)。這個結(jié)果與紅色毛癬菌(Trichophytonrubrum)的萌發(fā)過程相似[18]。

在這種萌發(fā)策略下，提高孢子產(chǎn)能利用效率極為重要，與糖酵解與糖代謝，賴氨酸代謝，果糖與甘露糖代謝，酪氨酸代謝,甘油酯代謝等途徑相關(guān)的基因均被顯著富集。糖類是最重要的生物大分子之一，不僅是細胞壁的組成成分，生命體重要的產(chǎn)能物質(zhì)，還廣泛參與生命體對溫度、水分、鹽分、紫外光和磷脅迫等非生物脅迫。對禾谷鐮刀菌(Fusariumoxysporum)和馬鈴薯粉痂菌(Spongosporasubterranea)轉(zhuǎn)錄組研究都表明真菌萌發(fā)時糖代謝活性變化十分劇烈[13,19]。在本研究中，低濕萌發(fā)的孢子果糖與甘露糖代謝顯著增強，糖酵解/糖異生代謝途徑的關(guān)鍵酶pfkA、PGK和PK基因全部上調(diào)表達，通過一系列酶促反應(yīng)，催化單糖產(chǎn)生ATP和丙酮酸，為孢子萌發(fā)提供能量。有氧條件下，ATP合酶、NADH-CoQ、SDH、NDH和COX等呼吸作用關(guān)鍵酶基因進一步上調(diào)表達，通過呼吸作用徹底氧化丙酮酸為二氧化碳和水，并釋放能量，增強了水分與能量的供給。

WA_75文庫中有大量絲氨酸和蘇氨酸蛋白激酶基因表達上調(diào)，通過復(fù)雜的信號轉(zhuǎn)導(dǎo)系統(tǒng)感受環(huán)境變化，對細胞的代謝和轉(zhuǎn)錄進行多靶向多層次調(diào)控，保證低濕萌發(fā)細胞活動的有序進行。本研究中發(fā)現(xiàn)的MPS1激酶是一種MAP-kinase，該基因喪失影響細胞壁完整性，CDKs、plo1激酶、DBF20激酶參與調(diào)控酵母等真核生物細胞有絲分裂等細胞周期活動，二酰甘油激酶家族(DGKs)調(diào)節(jié)脂肪代謝中甘油二酯和磷脂酸之間的平衡，在小麥抗脅迫鹽中起到了重要作用[20-22]。GCN2激酶(GCN2)是蛋白質(zhì)代謝關(guān)鍵調(diào)控因子，在細胞缺乏營養(yǎng)時調(diào)節(jié)氨基酸的內(nèi)穩(wěn)態(tài)，維持細胞存活[23]。鈣依賴性蛋白激酶(CDPKs)通過對逆境早期鈣信號的感應(yīng)和傳導(dǎo)，廣泛參與生物對非生物脅迫抗逆性反應(yīng)[24]。酪蛋白激酶I(CKI)廣泛分布在各種類型的細胞中，參與膜轉(zhuǎn)運、細胞周期、染色體分離、細胞凋亡等生命活動，有報道該激酶在水稻根系發(fā)育和激素調(diào)節(jié)中起到重要作用[25]。磷脂酰肌醇3-激酶和雷帕霉素靶蛋白(PI3K/mTOR)信號通路是免疫學(xué)研究的熱點之一，能量缺乏、營養(yǎng)不足、環(huán)境脅迫均可刺激該通路的信號傳導(dǎo)，調(diào)節(jié)細胞代謝，調(diào)控細胞凋亡[26]。

本研究采用RNA-Seq 技術(shù)分別構(gòu)建了球孢白僵菌WA菌株(B.bassianaWA)未萌發(fā)分生孢子、相對濕度100%和相對濕度75%萌發(fā)分生孢子轉(zhuǎn)錄組文庫WA_S、WA_100，WA_75，序列生物信息學(xué)分析結(jié)果表明：與WA_S相比，WA_100和WA_75孢子萌發(fā)過程中分別有822個和767個差異表達基因(DEGs)；經(jīng)GO功能富集分析，WA_75與WA_100的差異DEGs集中在單羧酸代謝進程(65)，分解代謝進程(72)，有機物分解代謝進程(65)，有機酸代謝進程(107)，酮酸代謝進程(106)等方面；KEGG功能富集在糖酵解與糖代謝(31), 賴氨酸代謝(12)，果糖與甘露糖代謝(10)，酪氨酸代謝(11),甘油酯代謝(9)等途徑。與WA_100相比，WA_75大量參與糖酵解、信號轉(zhuǎn)導(dǎo)的絲蘇氨酸蛋白激酶基因表達上調(diào)。本研究為闡明球孢白僵菌孢子耐低濕萌發(fā)的機制奠定了理論基礎(chǔ)。