松材线虫诱导棘孢木霉T203抗病相关基因的转录组研究

焦 宁，尹大川，宋瑞清

（1.东北林业大学林学院，哈尔滨150040；2.沈阳农业大学林学院,沈阳110161）

松树萎蔫病（pine wilt disease，PWD）是由松材线虫（Bursaphelenchusxylophilus）引发的松树毁灭性病害[1]。该病害在我国严重流行，对森林生态系统造成了巨大破坏，主要通过松墨天牛（Monochamusalternatus）等媒介昆虫传播，寄主主要包括黑松（Pinusthunbergii）、赤松（Pinusdensiflora）、马尾松（Pinusmassoniana）等松属植物，被线虫侵染的松树，针叶呈黄褐色或红褐色，整株迅速干枯死亡[2]。目前，针对松树萎蔫病的防治还是以化学药剂为主，而绿色环保的生物防治研究还处于起步阶段，本研究分析了对松材线虫具有抑制作用的棘孢木霉（Trichodermaasperellum），揭示了线虫诱导木霉防卫反应的分子机制，为生物防治松材线虫提供了新的思路。

棘孢木霉（T.asperellum）是一种广泛应用于植物病害领域的生防菌，木霉通过代谢产生的水解酶、抗生素等活性物质对病原物产生高效抑制的作用[3]。研究表明，棘孢木霉能够有效抑制终极腐霉菌（Pythiumultimum）生长[4]。张锴等[5]在哈茨木霉（Trichodermaharzianum）体内找到了与松材线虫相关的抗性基因P450。然而，目前对松材线虫诱导棘孢木霉分子抗病机制还鲜有研究，已报道的木霉生防基因主要包括几丁质酶、外切β-1，3-葡聚糖酶以及N-乙酰氨基葡萄糖苷酶等水解酶类[6]，其中还有很多诱导抗病机制的基因有待发掘。

近年来，高通量测序技术被广泛应用于真菌的组学研究当中。为更深入地研究松材线虫诱导棘孢木霉的分子抗病机制，对线虫处理过的棘孢木霉T203菌株进行了转录组分析，根据相关抗性基因的功能注释、结构域分析以及相关的文献报道，找到了棘孢木霉T203的抗病相关基因，这不仅丰富了棘孢木霉的基因表达谱，也为松材线虫诱导木霉的分子防卫反应机制提供了基础依据。

1 材料与方法

1.1 材料

棘孢木霉（T.asperellum）T203由东北林业大学林学院森林病理实验室提供。松材线虫（B.xylophilus）和灰葡萄孢菌（Botrytiscinerea）由中国林科院森环森保所提供。

1.2 主要仪器

电子分析天平、灭菌锅、分光光度计、高速离心机、超净工作台、荧光定量PCR仪、-80℃超低温冰箱、4℃冰箱。

1.3 培养基制备与供试菌株、线虫的培养

马铃薯葡萄糖琼脂培养基（potato dextrose agar，PDA）：纱布包裹的新鲜马铃薯碎块200g煮至沸腾，加入葡萄糖和琼脂各20g，过滤后使用蒸馏水定容至1L，分装至锥形瓶中，121℃高压灭菌30min备用。

使用平板PDA培养基培养棘孢木霉T203与灰葡萄孢，置于25℃恒温培养箱中14d，待平板内长满菌丝进行后续试验。

使用灰葡萄孢培养基培养松材线虫7d，待菌丝被线虫取食干净，培养基表面近光滑油状，表明线虫长满培养基准备收集使用。

1.4 松材线虫的收集与悬浮液制备

将培养好的线虫采用贝尔曼漏斗法[7]，静置6～8h，收集至离心管中，使用3500r·min-1离心机离心5min，使用PBST缓冲液冲洗3次，配制成20条·μL-1的线虫悬浮液，置于4℃冰箱保存。

1.5 供试菌株的处理与取样

选取长势优良的棘孢木霉用于试验，处理组（T）向木霉菌板中接入500μL线虫悬浮液（约10000头），对照组（CK）加入等量的无菌水，各进行3次生物学重复，样品编号分别为T1、T2、T3与CK1、CK2、CK3，于25℃恒温培养箱培养72h，滤掉线虫，将菌丝收集至1.5mL离心管中，每管菌丝约100mg，液氮处理后使用干冰运送至哈尔滨诺禾致源生物公司进行测序。

1.6 测序样品的数据质控

利用CASAVA将高通量测序数据转化为转录组序列数据。为保证数据的准确性，从原始数据中去除带接头及平均质量分数低于Q20的数据，通过数据过滤、错误率检测、GC分布情况，获得可靠的clean reads。

1.7 差异表达基因的筛选

进行测序深度校正，将原始的read count标准化处理，利用统计学模型进行假设检验概率（pvalue）的计算，多重假设检验校正后，得到padj值，依据|log2(FoldChange)|＞0且padj＜0.05的质控标准筛选差异表达基因[8]。

1.8 功能注释

使用cluster Profiler对差异基因进行GO功能富集和KEGG途径富集聚类，根据超几何分布原理将差异基因集注释到GO数据库和KEGG数据库中，根据基因的功能注释与代谢途径对差异基因集进行富集分析。

1.9 qRT-PCR验证

为验证转录组数据的准确性，选取12个显著差异的抗病相关基因进行qRT-PCR验证。使用TaKaRa Mini BESTUniversal RNA Extrction Kit提取处理组（T）与对照组（CK）样本菌丝总RNA。使用DNaseI（Promega）消化总RNA中的DNA，测定其质量浓度。反应条件如下：42℃60min；85℃5s；16℃10min。将合成的cDNA稀释10倍，用于qRT-PCR。扩增内参基因Actin[9]以及12个差异基因的引物序列如表1，使用CFX96 Real-time PCR Detection进行qRTPCR。反应体系如下：10μL 2×SYBR premix Ex Taq酶、正向和反向引物(10μmol·L-1)各1μL，2μL稀释的cDNA模板，加去离子水补足至20μL。反应系统如下：94℃预变性30s，94℃变性12s，60℃退火45s，72℃延伸45s，81℃读板1s，45个循环，每个样品重复3次。通过标准2-ΔΔCT方法计算每个基因的相对表达水平，使用SPSS19.0进行数据误差分析与差异显著性分析，所有试验数据采用独立样本t检验（p＜0.05）[10]。

表1 试验用引物序列Table 1 Primer sequences used in the experiment

2 结果与分析

2.1 测序数据质量

棘孢木霉T203的测序原始数据经Trimmomatic质控过滤后，采用FastQC质量评价。由表2可知，每个测序样品过滤掉低质量数据后的Clean data均超过95%，错误率为0.02%～0.03%，其中，Q20百分比为97%以上，Q30百分比也达到94%以上，GC含量接近50%，说明测序原始数据质量合格。

表2 测序数据质量评估Table 2 Quality evaluation of sequencing data

2.2 基因表达分布

通过计算各样本所有基因的表达量（FPKM），得到不同样本基因表达水平的分布情况，结果如图1。基因表达量的中位数比较接近，表明各样本的基因表达水平总体相似，但处理组（T）与对照组（CK）的基因表达情况也存在明显差异。

图1 各样本基因转录分布情况Figure 1 Gene transcription distribution of each sample

2.3 样本间相关性分析

样本间基因表达水平相关性是检验差异表达基因可信度的重要指标。使用FPKM方法检测各样本间的相关性，通过皮尔逊相关系数表示，相关系数在0.8～1之间属于极高相关。由图2可知，各样本相关系数均高于0.8，表明各样本间具有较好的重复性与极高的相关性，可信度高。

图2 样品相关性系数统计图Figure 2 Statistical chart of sample correlation coefficient

2.4 差异表达基因的筛选

差异表达基因筛选结果如图3，经松材线虫处理的棘孢木霉T203总共获得964个差异表达基因，其中包含371个上调表达基因，583个下调表达基因。

图3 差异表达基因火山图Figure 3 Differential expression gene Volcano plot

2.5 差异表达基因的聚类分析

对差异表达基因进行聚类分析，横向比较为同一基因在不同样本中的表达情况，结果如图4。处理组（T）中表达量高的基因在对照组（CK）中低，处理组（T）中表达量低的基因在对照组（CK）中高，处理组（T）与对照组（CK）表达模式相反。

图4 差异表达基因聚类热图Figure 4 Clustering heatmap of differentially transcribed genes

2.6 GO功能富集分析

根据差异表达基因的功能将GO功能富集分为3个大类，分别是细胞组分（cellular component）、分子功能（molecular function）和生物学过程（biological process）。细胞组分指细胞结构的成分，这类基因主要参与细胞活动；分子功能主要涉及催化活性与信号识别相关的基因；生物学过程包括生物代谢途径与信号传导运输相关的基因。GO富集分析结果如图5，其中，参与生物学过程（BP）的差异基因最多，其次是细胞组分（CC）、分子功能（MF）。生物学过程（BP）主要包括有氨基酸转运、有机氮化合物生物合成、肽生物合成、酰胺生物合成、蛋白质代谢等过程。细胞组分（CC）主要通过内质网膜蛋白、泛素蛋白、磷酸核糖激酶、转移酶复合物等诱导细胞质与细胞膜的合成。分子功能（MF）主要包括催化酶活性、连接酶活性以及聚合酶活性。

2.7 KEGG富集分析

KEGG总共富集到75个代谢途径，包括333个差异表达基因，其中105个基因显著上调表达，筛选出前20个最显著的KEGG代谢途径进行分析，结果如图6。这些差异表达基因主要富集在以下代谢途径：过氧化物酶体、氰氨基酸代谢、脂肪酸代谢、甲烷代谢、缬氨酸、亮氨酸和异亮氨酸的降解、谷胱甘肽代谢、甘油磷脂代谢、脂肪酸降解、精氨酸和脯氨酸代谢、淀粉和蔗糖代谢、碳代谢、次生代谢产物的生物合成磷酸戊糖途径、不饱和脂肪酸的生物合成、脂肪酸生物合成、氮素代谢、抗生素的生物合成、果糖和甘露糖代谢、组氨酸代谢、酪氨酸代谢。

图6 KEGG代谢途径富集分析Figure 6 Enrichment analysis of KEGG metabolic pathway

2.8 抗病相关的MAPK信号途径及基因

在KEGG富集代谢途径中，具有诱导抗病功能的分裂原活化蛋白激酶（MAPK）信号途径中共有5个基因显著上调表达，结果如表3。这些基因分别包括G蛋白信号调节蛋白、Mcm1转录因子、Rho GTPase激活蛋白、RNA聚合酶蛋白以及过氧化氢酶。

表3 MAPK信号途径差异表达相关基因Table 3 Differential expression related genes of MAPK signal pathway

2.9 棘孢木霉抗性基因筛选与qRT-PCR验证

通过基因注释，从显著上调表达的差异基因中筛选12个参与棘孢木霉防卫反应的抗病相关基因，结果如表4。这些基因分别为短链脱氢酶、细胞色素P450、bZIP转录因子、泛素结合酶、糖苷水解酶16家族、G蛋白信号调控因子、β-1，3葡聚糖酶、枯草杆菌蛋白酶、LysM受体激酶、醛酮还原酶、乙醛脱氢酶和过氧化氢酶。这些基因在松材线虫诱导的棘孢木霉T203转录组中显著上调表达。

表4 棘孢木霉抗性功能基因Table 4 Resistance genes of Trichoderma asperellum

通过荧光定量（qRT-PCR）对松材线虫处理的棘孢木霉T203菌株中12个抗病相关基因进行验证。将处理组（T）与对照组（CK）的基因表达情况同转录组测序数据进行比较，结果如图8。尽管某些基因的差异表达倍数存在偏差，但整体的上调趋势一致，基因表达水平与转录组数据接近，因此棘孢木霉T203的转录组测序数据具备可靠性。

3 讨论与结论

木霉菌生防效果显著，应用广泛，主要通过重寄生作用、竞争作用、抗生作用、促进植物生长发育以及诱导植物抗病机制等发挥其生防功能[11]。木霉菌在代谢过程中产生的几丁质酶、葡聚糖酶等降解酶类以及活性抗生物质对病原生物具有良好的抑制作用[12]。研究表明，木霉菌可有效抑制尖刀镰孢菌（Fusariumoxysporum）、立枯丝核菌（Rhizoctoniasolanisolani）等病原真菌的孢子萌发，对根结线虫也有一定的抗性作用[13-14]。本试验通过棘孢木霉作用松材线虫，同样发现随着处理时间的延长线虫死亡率逐渐上升，因此，证明了棘孢木霉对松材线虫具有很好的抑制作用。随着当代分子科技的迅速发展，对真菌转录组研究也越来越多[15]，KAWAHARA等[16]对稻瘟病菌与感病水稻分别进行了转录组测序分析，发现了水稻的植物保护素合成基因以及稻瘟病菌的致病基因。SUN等[17]对寄生在核盘菌（Sclerotiniasclerotiorum）上的粉红螺旋聚孢霉（Clonostachysrosea）的转录组测序分析中发现了与重寄生相关的基因。本研究通过RNA-Seq分析了松材线虫处理棘孢木霉T203菌株72h的差异表达基因，结果共检测到964个差异表达基因，这些差异基因富集到不同的GO功能及KEGG代谢通路中，反映了松材线虫诱导棘孢木霉T203在抗病过程中具有重要作用。

在棘孢木霉T203转录组的差异表达基因分析中发现，差异显著的基因主要集中在细胞壁与细胞膜的降解酶、细胞色素P450、糖苷水解酶、信号传导因子、MAPK信号途径、丝氨酸苏氨酸代谢途径、谷胱甘肽代谢途径、抗生素生物合成、过氧化物酶体等调控途径上。这些基因主要包括细胞壁降解酶、细胞膜降解酶、参与细胞壁修复途径的相关基因、参与抗生素代谢产物合成的相关基因、谷胱甘肽代谢途径中关键的解毒酶基因以及参与细胞抗氧化反应的相关基因。大量研究表明，水解酶代谢途径对病原生物的细胞壁与细胞膜具有降解作用[18]，细胞色素P450是一类高效的抗性基因[19]，丝氨酸苏氨酸代谢途径对维持细胞完整性、细胞周期、调节细胞渗透压、以及对外界生物及非生物因素影响都具有重要作用[20]。MAPK信号途径对调控真菌生长发育具有重要意义[21]，谷胱甘肽代谢途径与过氧化物酶体应对外缘刺激能自身诱导解毒代谢机制以及抗氧化反应[22]。总之，棘孢木霉在松材线虫刺激下所发挥的生物防卫反应机制是由多个基因参与共同调控作用的结果。

本研究旨在探索松材线虫诱导棘孢木霉T203的分子抗病机制，首先激发棘孢木霉对病原线虫入侵的信号识别机制，通过松材线虫与棘孢木霉菌丝的相互作用，诱导木霉体内合成水解酶以及抗生素等活性代谢产物，这种接触是病原物对生防菌菌丝的纤维素结合区域以及凝集素介导蛋白质表达模式的G蛋白和MAPK信号的联合反应[23]，该信号识别传递与诱导防卫反应机制在松材线虫与棘孢木霉的互作中得以充分体现。研究表明，木霉生防机制启动是对病原物的识别，而G蛋白在信号识别与传递中起到关键作用[24]，在识别成功后，木霉菌丝会附着在线虫体表，启动防卫反应机制，通过分泌水解酶、蛋白酶溶解线虫的细胞壁与细胞膜，产生抗生作用的次生代谢产物[25]，通过自身代谢反应诱导因线虫刺激引起的抗氧化反应与解毒抗逆机制，从而发挥抗病功能。

细胞壁是线虫对外的第一层防御系统，在维系细胞形态结构稳定方面具有重要意义，同时对外界不利环境提供很好的生物屏障[26]。SYAMA等[27]研究发现粘帚霉菌产生的β-1，3-葡聚糖酶不仅可以抑制镰刀菌（Fusariumsp.）和腐霉菌（Pythiumsp.）的孢子萌发与生长，还可以降解细胞壁、形成附着孢结构；LI等[28]发现粉红螺旋聚孢霉（Clonostachysrosea）的枯草杆菌蛋白酶基因PrC能够固定线虫并降解线虫体壁；WANG等[29]从洛斯里被毛孢（Hirsutellarhossiliensis）中提取出的枯草杆菌蛋白酶基因对线虫体壁也有降解作用。本试验在棘孢木霉T203的差异表达基因中发现了显著上调表达的糖苷水解酶16家族、β-1，3-葡聚糖酶以及属于丝氨酸蛋白酶家族的枯草杆菌蛋白酶，表明这些水解酶对降解松材线虫的细胞壁及细胞膜发挥着重要的作用。

外界有害生物的刺激往往导致生物细胞内活性氧（reactive oxygen species，ROS）过量累积，加速细胞膜脂过氧化反应，生成醇类、醛类等有毒物质，从而破坏核酸及蛋白质的正常结构和功能，造成细胞过氧化损伤[30]。本研究发现，在棘孢木霉T203的差异表达基因中有一些有关抗氧化与降解毒素的关键基因，包括过氧化氢酶、乙醛脱氢酶、谷胱甘肽S转移酶（GST）、醛酮还原酶等，其中谷胱甘肽S转移酶（GST）与醛酮还原酶在谷胱甘肽代谢途径中显著上调表达。以往研究表明，谷胱甘肽代谢途径可以激活相关抗性基因、诱导抗毒素的积累[31]，谷胱甘肽S转移酶（GST）对植物抵御外源生物具有重要作用[32]。KONG等[33]在黄瓜抗灰霉病的研究中筛选出27个谷胱甘肽代谢途径中的基因。以上表明，这些基因与其调控的代谢途径在松材线虫诱导棘孢木霉抗病过程中同样发挥着解毒代谢与诱导抗逆性的功能。

在研究菌丝生长发育与应激信号识别传导方面，本试验从KEGG显著富集的代谢途径中发现了至关重要的MAPK信号途径，相关研究表明，MAPK信号途径对细胞的糖基化状态敏感，在外界生物或非生物刺激下通过信号识别与传导激发相关抗逆蛋白酶活性，调节外源因素造成的细胞壁损伤，是维系真菌孢子萌发与菌丝正常生长发育的必要信号传导途径[34]，同时MAPK信号途径在真菌对病原物产生抗性作用、诱导防卫反应机制、形成气生菌丝与分生孢子等方面也具有重要作用[35]。经松材线虫处理棘孢木霉T203菌株的MAPK信号途径中，本试验发现，G蛋白调节因子显著上调表达，该基因可以促进菌丝生长并对线虫的刺激进行信号识别与传导；研究表明转录因子Mcm1可以促进细胞完整性，对细胞周期和毒力进行调控[36]；王亚等[37]在木霉菌H6菌株防治香蕉枯萎病的研究中发现过氧化氢酶活力升高与香蕉抗病性呈正相关，这些基因在MAPK信号途径中这些基因都显著上调表达，表明棘孢木霉的MAPK信号途径及其相关基因在菌丝生长发育、细胞信号传导运输以及在有害生物引起的过氧化反应中起着关键作用从而诱导棘孢木霉抗病机制。

图7 12个抗病相关基因相对转录水平Figure 7 Relative transcription of 12 disease resistance-related genes

通过基因功能及相关代谢途径分析，本试验在差异表达基因中筛选出12个棘孢木霉抗病相关功能的基因，并分别进行注释与qRT-PCR分析。结果表明，这些基因均显著上调表达，这与转录组测序数据的表达模式相似，表明这些基因在松材线虫诱导木霉抗病反应的分子机制中发挥着重要作用，同时也证实了转录组数据的可靠性。

从松材线虫处理棘孢木霉T203菌株的转录组数据中共筛选出964个差异表达基因，其中371个上调表达，583个下调表达。GO功能富集到氨基酸转运、蛋白质合成与代谢、细胞合成以及催化酶活性等相关功能的基因，KEGG代谢途径富集到抗生素生物合成、氨基酸代谢、谷胱甘肽代谢、过氧化物酶体等相关功能的基因。在KEGG代谢途径中，本研究发现了参与菌丝生长发育及信号识别传导的MAPK信号途径，并从中分析G蛋白、过氧化氢酶等有关信号识别与诱导抗病机制的关键基因。最后，依据差异表达基因的功能注释筛选出12个显著上调表达的抗病基因，它们主要编码信号传导因子、抗氧化酶、水解酶、解毒酶、抗生素等，并对这些基因进行了qRT-PCR验证。结果发现，12个基因显著上调表达且表达模式与测序数据趋势一致，表明这些基因在诱导棘孢木霉抗病机制中发挥着关键作用。