杨凌糖丝菌Hhs.015拮抗苹果树腐烂病菌转录组分析及抗菌机制探究

王建勋，康 凤，孙 玲，颜 霞，黄丽丽

(1.西北农林科技大学 生命科学学院，陕西杨凌 712100；2.国家旱区作物逆境生物学国家重点实验室，陕西杨凌 712100; 3.西北农林科技大学 植物保护学院，陕西杨凌 712100)

苹果是中国重要的经济作物，由苹果树腐烂病菌Valsamali引起的苹果树腐烂病严重制约了中国苹果产业的发展[1-3]。目前对于苹果树腐烂病主要依赖化学药剂防治，但是化学药剂残留给环境和人类带来了巨大的危害[4]。因此，利用环保型和可持续的生物菌剂保护果园免受苹果树腐烂病菌的侵害势在必行[5]。目前为止，已有一些用于防治Valsamali的有益微生物被发现，但防治仅处于初步阶段[5-7]。

放线菌是一类革兰氏阳性细菌，可以产生大量的代谢物和抗生素，能够对植物病原菌起到抑制的作用[8]。放线菌通过降解目标真菌细胞壁、产生抗真菌化合物来阻止病原菌的侵染[9-10]。杨凌糖丝菌(SaccharothrixyanglingeniesHhs.015，Hhs.015)是从黄瓜根部分离出来的一株稀有放线菌[11]，该菌在田间试验显示出对番茄叶霉病和苹果树腐烂病具有良好的防治效果[12-13]。有研究从Hhs.015发酵液中提取了活性物质异黄酮和戊霉素，以及2个戊烯大环类脂类物质WH01和WH02[14-15]。转录组是指特定组织或细胞在某个时间或某个状态下转录出来的所有RNA的总和，主要包括mRNA和编码RNA。转录组可以研究生物在不同环境、宿主等情况下的基因表达变化，阐述宿主或环境对生物的影响[16]。Kang等通过转录组分析，揭示了生防菌芽孢杆菌CC09在小麦体内抗病的分子机制[17]。Kröber等[18]利用转录组分析，发现生防菌BacillusamyloliquefaciensFZB42生物膜形成过程中，抗菌肽LCI基因和多种次级代谢物高效表达，揭示了生物膜的形成对解淀粉芽孢杆菌发挥生防作用的重要作用。

放线菌作为生防菌已被广泛研究，但是对于稀有放线菌属糖丝菌的抗菌机制研究相对较少。本研究对杨凌糖丝菌Hhs.015拮抗苹果树腐烂病菌Valsamali的过程进行转录组分析，试图探究Hhs.015在抗真菌过程中的分子机制，对阐明杨凌糖丝菌Hhs.015拮抗苹果树腐烂病菌的分子机理奠定一定的基础。

1 材料与方法

1.1 材 料

1.1.1 供试菌株 杨凌糖丝菌Hhs.015和苹果树腐烂病菌Valsamali均来自西北农林科技大学植物保护学院果树病害病原生物学及综合防治实验室。

1.1.2 主要培养基 马铃薯葡萄糖琼脂培养基(PDA)：200 g去皮马铃薯，切成小块加水煮沸30 min，然后用4层纱布过滤，在滤液中加入20 g葡萄糖和15 g琼脂粉，定容至1 L，高压灭菌后 备用。

黄豆培养基：称取20 g黄豆，用粉碎机打碎后煮沸30 min，4层纱布过滤，然后在滤液中加入10 g葡萄糖，5 g可溶性淀粉，2 g蛋白胨，2 g酵母膏，2 g NaCl，1 g CaCO3，0.5 g MgSO4·7H2O，0.5 g KH2PO4，15 g琼脂粉，定容至1 L，高温灭菌备用。所用试剂均为国产分析纯。

1.2 杨凌糖丝菌Hhs.015样品制备

将保存的杨凌糖丝菌Hhs.015采用划线法在黄豆培养基上进行活化培养，培养温度为 28 ℃。Valsamali则利用直径为10 mm的打孔器打取菌饼在PDA培养基上进行活化，培养温度为25 ℃。将活化好的Hhs.015直线转接至新鲜黄豆培养基中央28 ℃培养3 d，然后将活化好的Valsamali用无菌手术刀切成2～3 mm的长条置于Hhs.015两侧，共培养3 d(图1)，以空白PDA为对照，设置3个生物学重复。

图1 Hhs.015拮抗Valsa mali试验Fig.1 Hhs.015 antagonizes Valsa mali

1.3 杨凌糖丝菌RNA的提取、建库、数据质控

使用无菌牙签分别收集Valsamali拮抗培养和空白PDA处理的Hhs.015菌体，用锡箔纸包裹后立即置于液氮中，利用装有干冰的冰盒将收集的样品送至北京诺禾致源生物公司提取RNA、质检、建库以及转录组测序。

将质量合格的RNA用来构建文库，使用Qubit2.0进行初步定量，再用qRT-PCR对文库有效浓度进行准确定量，以保证文库质量。库检合格后，把不同文库按照有效浓度及目标下机数据量的需求pooling后进行Illumina测序。为了保证数据分析的质量及可靠性，对原始数据进行过滤。同时对Clean data进行Q20、Q30和GC含量计算，后续所有分析均是基于Clean data进行的高质量分析。使用R语言对样品间基因表达水平相关性进行Pearson相关系数计算，以评估RNA-seq的整体质量。相关系数越接近1，表明样品之间表达模式的相似度越高[19]。

1.4 差异表达基因(DEGs)分析

1.5 差异表达基因GO、KEGG注释

为了明确差异表达基因的生物学功能，阐明Hhs.015的抗菌机制。将差异表达基因进行Gene Ontology(GO)和Kyoto Encyclopedia of Genes and Genomes(KEGG)富集分析。即以Hhs.015基因组中所有基因的GO和KEGG注释为背景，使用超几何分布检验对DEG进行GO和KEGG途径富集分析。GO数据库分为三大类：细胞学组件(Cellular Component，CC)，分子功能(Molecular Function，MF)，生物学途径(Biological Process，BP)。通过R包软件GOseq实现差异表达基因的GO富集分析[21]，以校正后P<0.05的GO term作为差异表达基因显著富集，并绘制成散点图。

KEGG是一个重要的通路数据库，可以系统的分析基因功能和基因组信息，是生物代谢分析及代谢网络研究的重要工具[22]。通过KOBAS v2.0软件分析差异表达基因KEGG富集分析，结合GraphPad Prism软件绘制柱状图。

1.6 差异表达基因COG注释

Cluster of Orthologous Groups of Proteins(COG)，是由NCBI创建并维护的蛋白数据库，通过比对可以将某个蛋白序列注释到某一个COG中，从而推测该序列的功能。将差异表达基因进行COG注释，可明确差异表达基因所调节的生物学功能，从而分析Hhs.015在拮抗Valsamali过程中的差异基因功能分布。

1.7 DEGs中分泌蛋白及碳水化合物活性酶 分析

将差异表达基因利用SignalP-5.0 Server(http://www.cbs.dtu.dk/services/SignalP/)和TMHMM Server v2.0(http://www.cbs.dtu.dk/services/TMHMM/)在线工具预测其中的分泌蛋白。

Carbohydrate-Active enzymes Database是碳水化合物酶相关的专业数据库[23]，包括能催化碳水化合物降解、修饰、以及生物合成的相关酶系家族。包含五个主要分类：糖苷水解酶(Glycoside Hydrolases, GHs)、糖基转移酶(Glycol Transferases, GTs)、多糖裂解酶(Polysaccharide Lyases, PLs)、糖酯酶(Carbohydrate Esterases, CEs)和氧化还原酶(Auxiliary Activities, AAs)。此外，还包含与碳水化合物结合结构域(Carbohydrate-Binding Modules，CBMs)。使用Diamond软件，将差异表达基因与CAZy数据库进行比对，再结合分泌蛋白得到注释结果。

1.8 抗菌过程的次级代谢产物基因簇分析

利用antiSMASH(5.2.0)(https://antismash.secondarymetabolites.org/)预测Hhs.015中的次级代谢产物，并分析差异表达基因在次级代谢产物中的分布，明确Hhs.015在拮抗Valsamali过程中关键的次级代谢产物。

1.9 qRT-PCR验证DEGs

为了充分证实测序结果的可靠性，选取13个基因进行qRT-PCR验证，并以rpoA[24]作为内参基因。使用软件Primer Premier 5.0设计定量引物(表1)，以RevertAid RT逆转录试剂盒反转录得到cDNA作为模板，使用RealStar Green Mixture 进行实时荧光定量PCR检测。总体系为20 μL：cDNA 1 μL、Primer F/R 1 μL、RealStar Green Mixture 10 μL、ddH2O 7 μL。使用Roche lightcycler 96进行反应，反应程序为：95 ℃，10 min；45个循环(95 ℃，15 s；60 ℃，30 s；72 ℃，30 s)；溶解：95 ℃，15 s ；65 ℃，60 s；97 ℃，1 s。每组试验设置3个生物学重复，使用2-ΔΔCt法计算基因的相对表达量[25]。

表1 qRT-PCR检测基因及引物序列Table 1 Gene and primer sequence detected by qRT-PCR

2 结果与分析

2.1 测序数据质控分析

通过测序获得raw reads总数为108 587 756条，过滤后的clean reads总数为107 410 364条，clean reads的碱基总数为16.15 Gb(表 2)。样品误差率小于0.03%，Q20>97.50%，Q30> 92.88%，GC含量为70.56%～70.79%。从样本相关性热图(图2)可以看出，样本间相关性较好，取样合理。因此，本次测序数据质量优良，可以用于后续分析研究。

表2 测序结果与映射比率统计Table 2 Statistics of sequencing production and mapping ratio

横纵坐标为各样本相关系数的平方

2.2 差异表达基因分析

图3 差异表达基因火山图Fig.3 Volcano map of different expression genes

2.3 GO和KEGG富集分析

将GO分析的Term绘制成散点图(图 4)，可以看到上调表达基因主要集中在辅因子结合、水解酶活性、反转运蛋白活性、碳水化合物代谢等途径。下调表达基因主要集中在阳离子结合、草酸和羧酸代谢、类异戊二烯生物合成和代谢过程、脂质生物合成过程、谷氨酰胺家族氨基酸代谢过程、脂质代谢等过程。氧化还原过程及酶活性所涉及的基因一部分上调一部分下调。虽然最终表现为Hhs.015成功抑制了Valsamali的生长，但是在此过程中Hhs.015的代谢过程受到了Valsamali极大的影响。

图4 差异基因GO显著性富集分析Fig.4 Enrichment analysis of GO in differentially expressed genes

KEGG富集分析发现(图5)，ABC转运蛋白、次生代谢产物的生物合成、氨基酸的生物合成、代谢途径、抗生素的生物合成、核糖体等通路基因显著性上调表达。同时，ABC转运蛋白、淀粉和蔗糖代谢、次生代谢产物的生物合成、代谢途径等通路基因部分下调表达。可以看到，Hhs.015的淀粉和蔗糖代谢相关基因下调表达，阻碍了能量的供应，但是通过增强部分氨基酸代谢补充一部分能量的供应，满足Hhs.015生长的基本需求。

图5 差异基因KEGG分析Fig.5 KEGG analysis of differentially expressed genes

2.4 COG注释差异基因功能分析

将差异表达基因进行COG(Cluster of Orthologous Groups of proteins)注释，533个差异基因共有360个注释到COG数据库，包含241个上调表达基因，119个下调表达基因(图6)。在上调表达基因中“[Q]次级代谢产物生物合成，运输和分解代谢(42)”“[G]碳水化合物运输和代谢(42)”“[R]仅通用功能预测(35)”“[E]氨基酸运输和代谢(33)”“[I]脂质转运和新陈代谢(19)”这5类最多，占到70.95%。而类别“[E]氨基酸运输和代谢(16)”“ [T]信号转导机制(16)”“[R]仅通用功能预测(14)”“[K]转录(11)”“[M]细胞壁、膜包膜生物发生(10)”在下调表达基因中所占比例较高，达到56.30%。

图6 差异表达基因的COG注释分析Fig.6 Statistics of COG annotation for DEGs

可以进一步发现，在面对Valsamali的胁迫时，Hhs.015次级代谢产物合成、运输基因高效表达，起到关键的抗菌作用。同时碳水化合物、氨基酸和脂类的代谢增强，为Hhs.015的生长及次级代谢物的合成提供必要的能量。部分信号转导机制基因显著性下调，Valsamali分泌的物质破坏了Hhs.015信号的传递，还破坏了Hhs.015细胞壁、细胞膜的合成。

基因Hhs.015_007438，在[G]、[E]、[P]、[R]中均有注释，且为显著性上调表达，进一步分析发现，Hhs.015_007438是一个MFS(major facilitator superfamily)转运蛋白家族基因，还有另外2个显著性上调的MFS家族基因，Hhs.015_000830和Hhs.015_002439。

2.5 分泌蛋白及碳水化合物活性酶差异表达基因分析

差异表达基因结合分泌蛋白分析发现有16个上调表达的分泌蛋白编码基因，10个下调表达基因(表 3)。在上调表达的分泌蛋白中，包含α-甘露糖苷酶、ErfK/YbiS/YcfS/YnhG家族蛋白、铁配合物转运系统底物结合蛋白、ABC转运蛋白底物结合蛋白、α-葡萄糖苷酶等等。在下调表达分泌蛋白中，只有一个含有SH3结构域的蛋白，其他均为未知功能的假想蛋白。

表3 差异表达基因中的分泌蛋白及功能Table 3 Secretory proteins and functions in differentially expressed genes

结合碳水化合物酶数据库分析发现，差异表达基因中，有36个基因具有碳水化合物注释，27个基因上调表达，9个基因下调表达(图7)。分泌型的碳水化合物活性酶全部上调表达，主要集中在碳水化合物结合结构域和糖苷水解酶类。在非分泌型的碳水化合物活性酶中大部分为上调表达，同样也主要集中在碳水化合物结合结构域和糖苷水解酶类，部分还有氧化还原酶。下调表达的基因主要是糖基转移酶，还有部分碳水化合物结合结构域和糖苷水解酶类。

图7 分泌型碳水化合物活性酶(A)及非分泌性碳水化合物活性酶(B)的表达情况Fig.7 Secreted carbohydrate active enzyme(A) and non-secreted carbohydrate active enzyme(B) expression

2.6 抗菌关键次级代谢产物分析

差异表达基因结合次级代谢基因簇分析表明，在Hhs.015的28个次级代谢物基因簇中，七个基因簇在抗菌过程中上调表达的基因较多，分别是基因簇9、10、12、13、17、18、19，上调表达的基因个数分别为27、10、14、8、11、8、25(图8)。

图8 Hhs.015次级代谢物中的差异表达基因Fig.8 Different expression genes in Hhs.015 secondary metabolites

AntiSMASH[26]分析表明(表 4)，基因簇18与喷他霉素的相似性达到80%，该基因簇在拮抗Valsamali过程中有8个基因上调表达。喷他霉素属于多烯大环内酯类抗生素，具有较强的抗真菌特性，前期已经从Hhs.015中分离鉴定到pentaene macrolides，对Valsamali具有较好的抗性[14]。基因簇19与已知的macrotermycins具有96%的相似性，在拮抗过程中达到25个基因上调表达。曾有研究表明Macrotermycins A和C对人致病性金黄色葡萄球菌具有抗菌活性，它们还对白蚁真菌园的真菌寄生虫具有选择性的抗真菌活性[27]。然而剩下的Cluster 9、10、12、13、17与已知的物质相似性都比较低，表明Hhs.015在抗真菌过程中还有许多新物质在发挥作用，可见Hhs.015具有一定的天然抗生素开发潜力。

表4 上调表达基因集中的次级代谢物基因簇Table 4 Up-regulated gene in secondary metabolite gene cluster

2.7 qRT-PCR分析

从Hhs.015的差异表达基因中选择13个基因，利用qRT-PCR验证测序数据的可靠性。图9表明qRT-PCR检测结果与测序数据所分析的差异表达基因趋势一致。这13个基因包括4-羟苯基丙酮酸双加氧酶、丝氨酸3-脱氢酶、转录延伸因子GreA、包含CHAP域的蛋白、ATP结合蛋白、多糖生物合成酪氨酸自身激酶等等。

图9 qRT-PCR分析部分差异表达基因Fig.9 qRT-PCR analysis of some differentially expressed genes

3 讨论与结论

转录组分析在抗菌物质对微生物的抗性机制研究方面取得了不少成果，比如高梦莎研究发现杀菌剂臭氧通过破坏副溶血弧菌细胞膜的通透性、抑制活性氧清除酶的活性等过程杀死副溶血弧菌[28]。高亮等[29]研究壳聚糖处理辣椒疫霉后，通过抑制辣椒疫霉细胞膜、运动及糖代谢相关功能，从而起到了抑制辣椒疫霉生长的作用。本研究发现杨凌糖丝菌Hhs.015在拮抗Valsamali的过程中，Hhs.015的能量供应受到一定的阻碍，但是通过增强部分碳水化合物、氨基酸、脂类等物质代谢补充一部分能量的供应，为Hhs.015的生长及次级代谢物的合成提供必要的能量。Hhs.015的信号传递、细胞壁和细胞膜的合成等功能也受到Valsamali的抑制。与此同时，Hhs.015自身抗菌物质的合成和运输加强，比如一些次级代谢产物。ABC转运蛋白家族和MFS超家族几乎占据了生物所有转运蛋白的一半[30]。MFS转运蛋白具有营养吸收，代谢产物、有害物质排出等功能[31-33]，对于Hhs.015产生的抗生素运输和毒性物质排出具有重要作用。在拮抗过程中，ABC和MFS转运蛋白显著性上调表达是为了更好的适应环境而产生的一种自我保护机制，这对于提高产Hhs.015自身抗毒性物质和抗生素的运输具有重要意义。

Hhs.015在受到Valsamali的胁迫时，分泌型和非分泌型的碳水化合物结合结构域和糖苷水解酶类绝大部分上调表达。碳水化合物结合结构域可以提高酶的催化效率，为Hhs.015的生长以及抗生素的合成提供必要的能量。真菌细胞壁中多糖成分占到干质量的80%[34]，碳水化合物结合结构域可以特异性结合多糖底物，包括纤维素、几丁质、淀粉、糖原类等物质，糖苷水解酶的大量上调表达破坏了Valsamali细胞壁的合成，降低了Valsamali的毒性，同时能够高效的抑制Valsamali的生长[35]。还预测到在拮抗过程中有5种次级代谢物基因簇的基因大量上调表达，并且这5种次级代谢物与已知的物质相似性较低。糖丝菌属为稀有放线菌，近几年也不断的从糖丝菌属中发现了新的抗菌物质，比如Tianchimycin[36]、cyanogriside I 和cyanogriside J[37]，说明在Hhs.015抗真菌过程中可能还有新的物质在发挥着重要作用。

本研究通过转录组分析，揭示了杨凌糖丝菌Hhs.015拮抗Valsamali过程中发挥主要作用的通路表达情况，还预测到5种新的次级代谢物基因簇的基因在抗菌过程中大量上调表达。今后可以分离获得相关基因簇的次级代谢产物，有望获得高效的抗菌物质，为防治病原真菌提供新的材料。