球孢白僵菌分生孢子耐低湿萌发转录组相关差异基因分析

王子夜，张海剑，路晓月，王瑞东，王志刚，阎爱华

(1河北农业大学 林学院，河北 保定 071000；2河北省林木种质资源和森林保护重点实验室，河北 保定 071000；3河北省农林科学院 植物保护研究所，河北 保定 071000；4河北省农业有害生物综合防治工程技术研究中心，河北 保定 071000；5农业部华北北部作物有害生物综合治理重点实验室，河北 保定 071000; 6河北省城市森林健康技术创新中心，河北 保定 071000)

球孢白僵菌(Beauveriabassiana)是重要的昆虫病原真菌，其宿主范围广、特异性高、安全性好，在农林害虫生物防治中具有良好的应用前景。作为生防制剂的主要成分，分生孢子能否在不良环境胁迫下大量宿存、成功萌发、高效侵入是决定球孢白僵菌(B.bassiana)防治效果的关键[1-3]。空气相对湿度是昆虫病原真菌应用的限制性环境因子之一，RH90%～100%是分生孢子萌发最适湿度，很少昆虫病原真菌菌株能在低于RH70%的湿度下萌发[4-5]。探明不同湿度对白僵菌萌发的作用机理，对于拓展白僵菌的菌种资源，增强对害虫的防治效果，扩大生防制剂的应用范围具有重要意义。

对白僵菌-昆虫-环境三者的相互作用机制已经有大量研究，随着分子生物学、高通量测序技术等广泛应用，海藻糖酶(Beauveriabassiananeutral trehalase，BbNTH1)基因、特异性生长抑制蛋白(Growth arrest specific protein，GAS1)基因、热激蛋白70(Heat shock protein，HSP70)基因、丝裂原活化蛋白激酶(Mitogen-activated protein kinases，MAPK)同源基因slt2和hog1等与白僵菌附着、侵入钉形成、定殖和产孢相关的一系列信号传导、抗逆和形态建成等功能相关的基因被克隆和鉴定[6-11]。作为真菌生长发育的起始阶段，与萌发相关基因和调控机制的报道还很少[12-13]。本研究通过构建不同空气相对湿度下白僵菌菌株萌发差异表达基因(Differently expressed genes,DEGs)转录组文库，通过基因通路富集和功能注释，挖掘可能参与白僵菌低湿萌发的差异基因及其功能和代谢过程，为球孢白僵菌的进一步广泛应用提供理论基础，也为了解其他虫生生防菌萌发功能基因研究提供了重要线索。

1 材料与方法

1.1 供试菌株

供试球孢白僵菌(B.bassiana)菌株WA采自河北省保定市，由河北农业大学林木病虫害无公害防控实验室分离自光肩星天牛(Anoplophoraglabripennis)僵虫并进行保存。在萨式培养基(葡萄糖40 g，蛋白胨10 g，酵母膏10 g，琼脂20 g/L )上28 ℃活化培养15 d，收集分生孢子。

1.2 不同相对空气湿度分生孢子萌发率测定与文库样品采集

采用小容器空气湿度调节法测定空气相对湿度对白僵菌孢子萌发的影响，在2个口径210 mm玻璃密闭干燥器中分别注入灭菌饱和NaCl溶液和灭菌水，构建RH75%和RH100%的小空间气候室。载玻片上覆盖玻璃纸，将20 μL分生孢子悬浮液(约1.03×107个/mL)滴在玻璃纸上，迅速风干后置于密闭小空间气候室，25 ℃恒温培养24 h后取下玻璃纸，镜检记录观察分生孢子萌发率和芽管长度。同时大量收集不同处理分生孢子，液氮速冻，在零下80 ℃下保存，用于文库构建。样品分别命名为WA_75和WA_100，以未经处理且未萌发分生孢子为对照，命名为WA_S。

1.3 cDNA文库构建及转录组序列测定

cDNA文库构建由诺和致源生物信息技术有限公司完成，获得的序列在 Illumina HiSeqTM2000测序仪进行转录组测序。测序所得的原始片段(Raw reads)经过过滤，得到纯净片段(Clean reads)。纯净片段经Trinity软件进行Denovo 从头合成拼接，组装去冗余后得到Unigenes进行数据分析。

1.4 差异表达Unigenes数据注释与分析

差异Unigenes的表达丰度用FPKM值表示[q value<0.05 & |log2(fold change)|>1]，序列信息分别在Nr数据库、NT数据库,Swiss-Prot数据库、KEGG数据库、COG数据库和GO数据库中比对和检索，获得与Unigenes对应的蛋白功能注释及功能分类统计。差异基因表达热图采用欧易云在线软件(https://cloud.oebiotech.com)进行。

1.5 差异表达基因实时荧光定量检测

为验证不同湿度白僵菌萌发转录组文库构建质量，以白僵菌18SrRNA基因为内参基因，随机挑选了文库中5个DEGs序列，用Probegene设计RT-qPCR引物(表1)，验证差异基因表达趋势[14]。按照材料与方法1.2收集萌发的WA孢子，经过RNA提取并反转录各样品的cDNA，每处理3个重复，按2-ΔΔCt相对定量法计算基因相对表达量。

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

1.6 数据统计

采用Excel 2003进行数据处理，采用SPSS 17.0进行单因素ANOVA检验分析，满足P<0.05,再进行方差分析和多重比较，并用Origin 2021软件作图。

2 结果与分析

2.1 不同湿度下菌株萌发情况测定

不同湿度下菌株萌发情况测定，见图1。

图1 不同空气相对湿度下分生孢子萌发率和萌发芽管平均长度

如图1所示，WA菌株具有较强的耐低湿萌发能力，在RH100%和RH75%时，菌株WA的分生孢子萌发率和萌发芽管长度没有显著差异(P<0.05)。

2.2 转录组测序数据质量分析与注释

各样品数据质量分析，见表2。

表2 样品数据质量分析Table 2 Summery of sample sequence data output quality

如表2所示，WA_S、WA_100和WA_75转录组文库数据质量评估去除不可靠数据后，3个样品各获得纯净片段为43 016 924条，44 723 918条和43 537 642条，Q30达到92%以上，错配率控制在1%以下，碱基A、T、C和G之间没有明显的偏差。组装的转录本和Unigene平均长度分别为547 bp和445 bp，测序和Unigene的拼接质量较好。序列分别在Nr、GO、KOG、Swiss-Prot、KEGG、Nt 6个数据库数据库中注释，共注释到89 903个Unigenes，其中27 786个Unigenes可以注释到至少1个文库中，数据库成功匹配率达到50.43%。

2.3 不同空气相对湿度下不同菌株DEGs分析

差异基因韦恩图分析，见图2。

图2 差异基因韦恩图分析Figure 2 Representing and Venn diagram of DEGs in response to different humidity

如图2所示，与未萌发孢子相比，不同空气相对湿度下WA菌株DEGs数量如图，WA_75与WA_S共有差异表达基因767个，其中上调342个，下调425个，WA_100与WA_S共有差异表达基因822个，上调400个，下调422个。393个DEGs在WA_75与WA_S和WA_100与WA_S中共同表达，其中上调114个，下调279个。

2.4 分生孢子萌发差异表达基因GO注释结果

差异最显著30条GO功能注释通路，见图3。

如图3所示，不同湿度下孢子萌发差异表达基因在GO文库注释到细胞组分(Cellular component，CC)，生物进程(Biological process，BP)，分子功能(Molecular function，MF)3个大类44个亚类。与未萌发的孢子相比，不同湿度条件下萌发孢子的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 分生孢子萌发差异表达基因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，确定不同基因间互作参与的生化代谢途径和信号转导途径。文库中共有23 549条DEGs注释到KEGG代谢通路中，分为细胞过程(Cellular processes)，环境信息处理(Environmental information processing)，遗传信息处理(Genetic information processing)，代谢(Metabolism)，有机系统(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)等方面数量最多(P<0.05)。

2.6 分生孢子萌发差异基因表达模式分析

2.6.1 呼吸作用相关基因 呼吸作用相关基因表达模式，见图5。

图5 呼吸作用相关基因表达模式Figure 5 Expression patterns of respiration related DEGs

如图5所示，在有氧条件下，白僵菌萌发时呼吸作用剧烈，聚类分析结果表明，氧化磷酸化相关基因在WA_75和WA_100中表达模式相似，但ATP 合成酶(ATP synthase)、NADH脱氢酶(NADH dehydrogenase，NADH-CoQ)、琥珀酸脱氢酶(Succinate dehydrogenase，SDH)、NADH-泛醌氧化还原酶(NADH-ubiquinone oxidoreductase，NDH)和细胞色素c氧化酶(Cytochromecoxidase，COX)相关序列在WA_75中上调。

2.6.2 信号转导蛋白激酶基因表达模式分析 蛋白激酶基因表达模式分析，见图6。

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

如图6所示，蛋白激酶催化蛋白磷酸化，是多种信号转导途径中的重要环节。聚类热图分析表明，低湿萌发过程中，WA_75文库中大量丝氨酸和苏氨酸蛋白激酶基因表达上调，包括BR激酶(BR Serine/threonine-protein kinase,BRSK)、DBF20激酶(Serine/threonine-protein kinase，DBF20)、GCN2激酶(Serine/threonine-protein kinase，GCN2)、钙依赖性蛋白激酶(Calcium-dependent protein kinases,CDPKs)、生发中心激酶(Germinal center kinase，GCKs)、酪蛋白激酶I(Casein kinase I，CKI)；磷酸肌醇3激酶(Phosphoinositide 3-kinase，PK3)、二酰甘油激酶(Diacylglycerol kinase)等；参与细胞周期调控的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技术验证随机挑选的5个DEGs表达量的差异，各基因FPKM值和相对表达量虽然在变化程度上有差异，但不同基因均呈上升趋势。细胞壁胁迫响应组份蛋白(WSC domain-containing protein)、假定泛素蛋白连接酶(Putative RING-H2 finger protein)和Zn (II)2Cys6锌簇蛋白[Zn (II)2Cys6 transcription factor]的FPKM值变化程度较大，锌指和BTB结构域蛋白(Zinc finger and BTB domain-containing protein)和假定泛素蛋白连接酶(Putative RING-H2 finger protein)FPKM值变化程度较小；假定外切聚磷酸酶相对表达量变化程度最大，其他4个基因变化差异较小。

3 讨论与结论

尽管温度、光照、湿度都被认为是影响孢子萌发的关键环境因素，但实践应用证明，湿度是限制白僵菌使用范围和效果的决定性因子[15]。本研究中球孢白僵菌孢子在不同湿度下萌发时都有大量差异基因表达。分生孢子生理活性状态由自身的养分状况决定，完全成熟的分生孢子在静止状态时仍存在转录活性，能够感受环境条件以进行主动的萌发决策[16]。利用转录组学技术探讨球孢白僵菌孢子在不同空气湿度条件下萌发差异基因表达模式，不仅可以揭示白僵菌相应水分胁迫的调控机制，而且为虫生真菌进一步扩大应用研究提供了良好的切入点。

侵染早期基因表达模式的差异是导致球孢白僵菌不同菌株致病性差异的原因，Wang等研究表明毒力与白僵菌的粘附、角质层降解相关功能基因的高表达密切相关[17]。萌发出芽管是真菌生长发育的开端，是菌落形态建成中第一个重大变化，在本研究中，不同湿度下孢子萌发时氧化还原进程，真菌细胞壁组份合成和蛋白质翻译进程被显著富集，为白僵菌快速萌发提供物质基础。同时，比WA_100，不利的水分条件WA_75抑制了有机氮化物、一元羧酸和氨基酸代谢进程等，但增强了对自身营养和资源的利用效率。单羧酸代谢进程、有机物分解代谢进程、有机酸代谢进程、酮酸代谢进程等相关生物学进程基因表达活跃，在不利条件下维持正常萌发。这个结果与红色毛癣菌(Trichophytonrubrum)的萌发过程相似[18]。

在这种萌发策略下，提高孢子产能利用效率极为重要，与糖酵解与糖代谢，赖氨酸代谢，果糖与甘露糖代谢，酪氨酸代谢,甘油酯代谢等途径相关的基因均被显著富集。糖类是最重要的生物大分子之一，不仅是细胞壁的组成成分，生命体重要的产能物质，还广泛参与生命体对温度、水分、盐分、紫外光和磷胁迫等非生物胁迫。对禾谷镰刀菌(Fusariumoxysporum)和马铃薯粉痂菌(Spongosporasubterranea)转录组研究都表明真菌萌发时糖代谢活性变化十分剧烈[13,19]。在本研究中，低湿萌发的孢子果糖与甘露糖代谢显著增强，糖酵解/糖异生代谢途径的关键酶pfkA、PGK和PK基因全部上调表达，通过一系列酶促反应，催化单糖产生ATP和丙酮酸，为孢子萌发提供能量。有氧条件下，ATP合酶、NADH-CoQ、SDH、NDH和COX等呼吸作用关键酶基因进一步上调表达，通过呼吸作用彻底氧化丙酮酸为二氧化碳和水，并释放能量，增强了水分与能量的供给。

WA_75文库中有大量丝氨酸和苏氨酸蛋白激酶基因表达上调，通过复杂的信号转导系统感受环境变化，对细胞的代谢和转录进行多靶向多层次调控，保证低湿萌发细胞活动的有序进行。本研究中发现的MPS1激酶是一种MAP-kinase，该基因丧失影响细胞壁完整性，CDKs、plo1激酶、DBF20激酶参与调控酵母等真核生物细胞有丝分裂等细胞周期活动，二酰甘油激酶家族(DGKs)调节脂肪代谢中甘油二酯和磷脂酸之间的平衡，在小麦抗胁迫盐中起到了重要作用[20-22]。GCN2激酶(GCN2)是蛋白质代谢关键调控因子，在细胞缺乏营养时调节氨基酸的内稳态，维持细胞存活[23]。钙依赖性蛋白激酶(CDPKs)通过对逆境早期钙信号的感应和传导，广泛参与生物对非生物胁迫抗逆性反应[24]。酪蛋白激酶I(CKI)广泛分布在各种类型的细胞中，参与膜转运、细胞周期、染色体分离、细胞凋亡等生命活动，有报道该激酶在水稻根系发育和激素调节中起到重要作用[25]。磷脂酰肌醇3-激酶和雷帕霉素靶蛋白(PI3K/mTOR)信号通路是免疫学研究的热点之一，能量缺乏、营养不足、环境胁迫均可刺激该通路的信号传导，调节细胞代谢，调控细胞凋亡[26]。

本研究采用RNA-Seq 技术分别构建了球孢白僵菌WA菌株(B.bassianaWA)未萌发分生孢子、相对湿度100%和相对湿度75%萌发分生孢子转录组文库WA_S、WA_100，WA_75，序列生物信息学分析结果表明：与WA_S相比，WA_100和WA_75孢子萌发过程中分别有822个和767个差异表达基因(DEGs)；经GO功能富集分析，WA_75与WA_100的差异DEGs集中在单羧酸代谢进程(65)，分解代谢进程(72)，有机物分解代谢进程(65)，有机酸代谢进程(107)，酮酸代谢进程(106)等方面；KEGG功能富集在糖酵解与糖代谢(31), 赖氨酸代谢(12)，果糖与甘露糖代谢(10)，酪氨酸代谢(11),甘油酯代谢(9)等途径。与WA_100相比，WA_75大量参与糖酵解、信号转导的丝苏氨酸蛋白激酶基因表达上调。本研究为阐明球孢白僵菌孢子耐低湿萌发的机制奠定了理论基础。