解淀粉芽孢杆菌群体遗传特征分析

王依婷，李素珍，毛志涛，廖小平，许慧敏，邱艺华，闫雪崧，杨凤妍，李贞景，郭庆彬*，刘欢欢*

（1.天津科技大学 食品科学与工程学院，天津 300457；2.思念食品有限公司，河南 郑州 450000；3.中国科学院 天津工业生物技术研究所，天津 300308；4.滨州市沾化区海洋和渔业发展服务中心，山东 滨州 256600）

解淀粉芽孢杆菌（Bacillus amyloliquefaciens）是一类革兰氏阳性菌，最初被定义为可产α-淀粉酶的枯草芽孢杆菌，广泛用于生物防治、传统发酵食品等领域[1-3]。B.amyloliquefaciens 具有较广的抑菌谱，能通过自身代谢产生抑菌蛋白、聚酮化合物等来达到抑制病原菌生长繁殖的作用。刘冬等[4]从大豆根际土壤分离得到的菌株JDF3，可抑制大豆疫霉孢子囊的产生，对大豆疫病的防治效果达到70.70%。Rasiya 等[5]从印度红树内生细菌中分离到的菌株RKEA3，其产生的抗菌脂肽能够抑制Staphylococcus aureus、Streptococcus pyogenes和Cryptococcus neoformis。同时，B.amyloliquefaciens 也可用于酱油[6]、白酒[7]、豆豉[8]等传统食品的发酵。此外，B.amyloliquefaciens 也被发现对黄曲霉毒素B1 诱导的小鼠盲肠炎症具有保护作用[9]，表现出一定的益生活性。

由于B.amyloliquefaciens 在农业与食品领域中的广泛应用，越来越多的菌株已完成了基因组测序。截止目前，NCBI 数据库已收录超过150 个基因组序列，为更好地理解该菌株的生理生化特征以及工业化提供了支撑，同时也为解析该物种群体遗传的多样性和系统进化以及功能基因的多态性、基因组变异与表型关联、正向基因筛选等提供了数据基础。近些年发展起来的泛基因组学技术（pan-genome）[10]，通过群体内大量的基因组来完整描述一个物种全部的遗传信息[11]，侧重于分析特定菌株基因的存在/缺失对种群功能的获得/缺失变异（presence/absence variation，PAV）之间的关联。其中，该种群所有样本共有的基因称为核心基因，一般与物种生物学功能和主要表型特征相关；仅在部分样本中存在的基因称为附属基因，一般与物种对特定环境的适应性或特有的生物学特征相关，反映了部分物种的特性。泛基因组是全面解析物种多样性和研究表型差异背后分子机制的强有力工具[12]。

本文从泛基因组学角度出发，通过Prokka 基因组注释[13]、Roary 泛基因组构建[14]、ProCAST 蛋白功能注释、antiSMASH 次级代谢产物基因簇注释[15]以及ABRicate 抗生素耐受性基因注释等分析流程，构建B.amyloliquefaciens 的泛基因组并进行深入分析，挖掘B.amyloliquefaciens 物种的普遍性遗传特征及碳水化合物的利用、植物促生、生物防治等方面的遗传特异性，以期为深度利用B.amyloliquefaciens 这种微生物资源作参考。

1 材料与方法

1.1 材料

从NCBI 数据库（https：//www.ncbi.nlm.nih.gov/datasets/genomes）下载组装水平标识为“Complete”的66株B.amyloliquefaciens 基因组序列（fasta 格式，截止2022 年8 月17 日）。

1.2 方法

1.2.1 平均核苷酸一致性分析

平均核苷酸相似度（average nucleotide identity，ANI）是在核苷酸水平比较两个基因组亲缘关系的指标。FastANI 可将基因组序列分割为短序列片段[16]，使用Mashmap 算法计算同源映射并估计ANI，在本研究中，参数fragLen 设定为500 bp。之后对所有基因组的ANI 数值矩阵在R 语言中进行层级聚类，用iTOL v6[17]（https：//itol.embl.de/）绘制聚类树图。

1.2.2 泛基因组构建及分析

采用Prokka 注释流程[13]及默认参数，对所有基因组进行基因结构预测和功能注释，其中获得的gff 格式文件作为泛基因组的输入文件。Roary 通过blastp 对gff文件中的序列进行直系同源基因分析[14]，获得B.amyloliquefaciens 泛基因组（pan-genome）和核心基因组（core genome）。采用Origin 软件对泛基因组/核心基因组编码基因数目与基因组数目之间的关系进行拟合[18]，拟合公式如下。

y=A×xB+C

式中：y 为泛基因组编码基因数目；x 为新增基因组数目；A、B、C 为系数。对于核心基因组的变化，设新增核心基因组数目为x，泛基因组大小为y，A、B、C 为系数，拟合公式如下。

y=A×eBx+C

使用mafft 对核心基因组进行序列比对，并采用FastTree 及GTR 算法对比对文件构建基于核心基因组序列的系统发育树[19]。同时，对Roary 生成的泛基因组PAV 矩阵，进行层级聚类并构建树图。

1.2.3 核心基因组功能注释及富集分析

利用ProCAST 平台（https：//procast.biodesign.ac.cn）对Roary 分析获得的核心基因组蛋白序列进行功能注释。该管路是由中科院天津工业生物技术研究所开发的综合蛋白质注释的无服务器网络工具，涵盖蛋白家族、域、分类、结构、基因本体和途径信息等26 个数据库，并提供多个工具用于结果的可视化分析。

1.2.4 质粒、抗生素耐药性等基因的鉴定

采用ABRicate 流程（Seemann T，Abricate，https：//github.com/tseemann/abricate），整合CARD[20]、BacMet[21]、PHAGE、PlasmidFinder[22]、NCBI AMRFinderPlus[23]数据库，对66 个菌株基因组的所有蛋白编码序列（coding sequence，CDS）进行功能注释。其中，CARD 与NCBI AMRFinderPlus 为抗生素耐药性信息数据库；BacMet为抗菌杀菌剂和金属抗性基因数据库；PHAGE 为噬菌体基因数据库；PlasmidFinder 包含质粒复制子序列，可用于原始和组装测序数据的复制子序列分析。

1.2.5 次级代谢产物基因簇的鉴定

采用antiSMASH 本地化软件流程[15]，对66 个菌株基因组注释产生的gbk 文件进行次级代谢产物基因簇挖掘，参数采用--cb-general、--cb-subclusters、--cbknownclusters、--asf 及--pfam2go。

1.2.6 碳水化合物水解酶基因的鉴定

采用本地化的CAZy 数据库对66 个菌株基因组的所有蛋白编码序列（CDS）进行功能注释[24]。

1.2.7 解磷途径及吲哚乙酸（indoleacetic acid，IAA）合成途径相关基因的鉴定

根据京都基因与基因组百科全书（Kyoto Encyclopedia of Genes and Genomes，KEGG）中IAA 合成途径中的蛋白质酶编号[25]，以及解磷途径涉及的酶编号[26]，在UniProt 数据库中检索已校验的fasta 蛋白序列，并通过blast 建库，对66 个基因组蛋白序列进行比对，identity阈值设定为40%，鉴定潜在的IAA 合成和解磷途径。

2 结果与分析

2.1 平均核苷酸一致性分析结果

ANI 被定义为两个微生物基因组同源片段之间平均的碱基相似度，是从全基因组水平评估物种间亲缘关系，尤其是近缘物种的重要工具，通常以ANI>95%作为划分物种边界的阈值[16]。本研究通过使用fastANI对66 个B.amyloliquefaciens 基因组进行同源性评估，结果如图1 所示。

图1 基于ANI 分析的B.amyloliquefaciens 层级聚类树Fig.1 Hierarchical clustering based on ANI values from B.amyloliquefaciens genomes

如图1 所示，所有基因组可聚类为2 个主要分支Clade I 与II，其中，在分支内两两基因组之间的ANI数值均高于95%，而分支间的两个基因组ANI 数值均高于90%而低于95%[16]。因此，Clade I 和Clade II 之间已经出现较为明显的物种边界，未来可将这2 个分支进行独立命名。

2.2 泛基因组功能分析

利用Prokka 对所有基因组进行基因结构预测，结果见表1。

表1 B.amyloliquefaciens 基本信息Table 1 Basic information of B.amyloliquefaciens

由表1 可知，B.amyloliquefaciens 的基因组介于3.70～4.44 Mb，最小为RD7-7，最大为XJ5，GC 含量为45.50%～46.50%，编码基因数目为3 716～4 445 个，平均每个基因组编码3 997 个基因。结合Roary 泛基因组分析流程，对上述基因组的核心基因组和泛基因组进行鉴定，结果见图2。

图2 核心基因数量和总基因数量随基因组数目增长的变化趋势Fig.2 Variation trend of core gene number and total gene number with the increase of genome number

由图2 可知，B.amyloliquefaciens 泛基因组基因总数为13 051 个，核心基因组中包含基因2 165 个，占基因组平均编码基因数目的54.16%。

由于物种受到生存环境与生态位差异、有效群体规模、多态性水平差异等因素的影响，泛基因组呈现开放型和闭合型两种模式。开放型泛基因组特点是物种内随着测序菌株的增加，基因个数也随之增加，较难到达平台期；闭合型泛基因组特点是物种内随着有限的菌株个体增加，基因个数迅速达到平台期[27-28]。

在对泛基因组编码基因数目与基因组数目之间的拟合公式中，y=3 925x0.29-343.8（R2=0.998）。其中，0.29为非零正数，泛基因组中基因总数会随着新的基因组的加入而呈指数性增加，因此推断B.amyloliquefaciens的泛基因组呈现一定的开放性。

核心基因组编码基因数目与基因组数目之间的拟合公式为y=1 388e-0.1571x+2 235（R2=0.980）。根据拟合方程，随着B.amyloliquefaciens 基因组测序数量的增加，核心基因组编码基因数目可能稳定在2 235 个。上述分析表明，B.amyloliquefaciens 易从外界获得新的基因，但同时核心基因组具有较强的稳定性。

泛基因组的缺失与存在变异（PAV）进行聚类分析，泛基因组中核心基因序列的系统发育分析结果如图3 和图4 所示。

图3 B.amyloliquefaciens 基因存在/缺失变异（PAV）矩阵的层次聚类Fig.3 Hierarchical clustering based on gene presence/absence variation（PAV）matrix of B.amyloliquefaciens

图4 基于B.amyloliquefaciens 核心基因组序列的系统发育树Fig.4 Phylogenetic tree based on core genome sequence alignment of B.amyloliquefaciens

由图3 与图4 可知，此结果与ANI 结果一致，表明66 株菌在全基因组ANI 相似性、PAV 聚类以及核心基因组序列相似性中均可被明确划为两个分支，Clade II中的15 株菌与Clade I 中的51 株菌不属于一个种，将来可作为独立物种进行研究。

为进一步解析该物种所具有的核心基因功能特征，即菌株的普遍性遗传特征，本研究采用ProCAST中集成的21 个数据库对核心基因组蛋白序列进行了注释。其中NR、SwissPort 和PFAM 数据库的注释基因占核心基因数的比例均超过90%。

基于COG 数据库的功能注释结果如图5 所示。

图5 核心基因组的COG 注释结果Fig.5 Main annotations for COG annotation

如图5 可知，R（11.14%）、E（8.56%）、J（7.30%）、K（7.03%）、P（5.71%）、C（5.65%）、M（5.10%）、G（4.83%）等占据了主要的功能模块。核心基因组的PATHWAY注释结果见图6。

图6 核心基因组的PATHWAY 注释结果Fig.6 Main annotations for PATHWAY annotation

在PATHWAY 注释中（包含KEGG pathway、Reactome、MetaCyc 等数据库），Aminoacyl-tRNA biosynthesis、Mitochondrial translation elongation 以及Purine metablism 等模块聚集了大量的功能基因，显示出B.amyloliquefaciens 在转录、翻译、氨基酸代谢、离子运输、能量代谢等基础代谢方面功能保守。

CAZy 是一个整合了碳水化合物酶活性的数据库。B.amyloliquefaciens 的核心基因组注释结果见图7。

图7 核心基因组的CAZy 注释结果Fig.7 Main annotations for CAZy annotation

如图7 所示，糖苷转移酶类（glycoside transferases，GT）数量为50 个，糖苷水解酶（glycoside hydrolases，GH）数量为49，占CAZy 酶总数的37.04%和36.30%。GH 和GT 是参与糖的合成和代谢的主要酶类，其中GH13 家族具有α-淀粉酶活性，占CAZy 酶总数的5.19%；GH3 和GH51 两个家族具有β-葡聚糖酶活性，占CAZy 酶总数的3.70%；GH3 和GH43 两个家族具有木聚糖酶活性，共占CAZy 酶总数3.70%。表明B.amyloliquefaciens 具有较好的水解淀粉、β-葡聚糖和木聚糖潜力。这为B.amyloliquefaciens 的糖苷水解酶和糖苷转移酶类开发利用提供了遗传基础。已有研究表明，B.amyloliquefaciens 可以产生β-1，3-葡聚糖酶等病程相关蛋白，通过作用于真菌细胞壁的不同成分来破坏真菌细胞壁，以达到抵抗病原真菌侵染的目的[29-30]。核心基因组的基因本体论注释结果见图8。

图8 核心基因组的GO 注释结果Fig.8 Main annotations for GO annotation

如图8 所示，在基因本体论（Gene Ontology，GO）注释中，B.amyloliquefaciens 大量核心基因聚集在oxidoreductase activity（7.14%）、structural constituent of ribosome（5.42%）、integral component of membrane（2.97%）等涉及氧化还原、遗传信息传递和生物膜合成等生物学过程模块；在细胞组成模块中，主要涉及catalytic activity（9.25%）、protein binding（1.98%）、cytoplasm（0.66%）等；而在分子功能模块中，大量功能基因聚集在DNA binding（10.58%）、catalytic activity（9.25%）、transmembrane transporter activity（7.87%）等方面，表明B.amyloliquefaciens 在催化分解、蛋白质的合成以及转运方面较为保守。上述核心基因组的功能注释证明B.amyloliquefaciens 在基础代谢相关的基因中显示出一定的遗传和功能保守性，为研究该菌的普遍代谢特征提供了基础。

此外，利用TMHMM 数据库对蛋白质跨膜螺旋结构域进行预测，同时采用SIGNALP 数据库对氨基酸序列中的信号肽进行注释。将含有信号肽且无跨膜结构域的蛋白质识别为分泌蛋白。结果显示，在B.amyloliquefaciens 的核心基因组中，共识别出99 个分泌蛋白，占核心基因总数的4.57%，并且在分泌蛋白中发现了溶杆菌素的蛋白质序列，表明B.amyloliquefaciens普遍能够分泌溶杆菌素这一抗菌肽，赋予了该物种广谱抗菌的特性。

2.3 所有B.amyloliquefaciens 菌株基因组功能模块的PAV 分析

因B.amyloliquefaciens 在生物防治以及食品等领域应用广泛，因此对所有66 株菌的CAZy、解磷途径、IAA 合成途径、次级代谢产物基因簇以及抗生素耐药性基因进行全面比较分析，为比较B.amyloliquefaciens种群内生物功能提供理论基础。

将66 株B.amyloliquefaciens 与其他芽孢杆菌（例如：B.subtilis 168、B.pumilus Ha06YP001 等）对比，结果见图9。

图9 CAZyme 编码基因的分布Fig.9 CAZyme coding gene distribution

如图9 所示，B.amyloliquefaciens GH3、GH43 以及GH51 基因数较多。其中，GH3 和GH51 家族为β-葡聚糖酶，能催化水解谷物细胞壁中的β-葡聚糖[29]；GH43 是一类降解木聚糖的酶系，可降解自然界中大量存在的木聚糖类半纤维素。此外，B.amyloliquefaciens ARP23 等30 株菌的GH13 家族编码基因数普遍较多，而GH13 家族为支链淀粉酶，能够催化淀粉等多糖化合物中的α-1，6-糖苷键水解，形成直链淀粉。以上结果表明66 株B.amyloliquefaciens 具有丰富的β-葡聚糖和木聚糖降解酶的基因资源。依赖色氨酸的IAA 生物合成途径见图10。解磷途径和IAA 合成途径相关酶的基因簇分布见图11。

图10 依赖色氨酸的IAA 生物合成途径Fig.10 Tryptophan dependent IAA biosynthesis pathway

图11 解磷途径和IAA 合成途径相关酶的基因簇分布Fig.11 Gene distribution of enzymes related to phosphorus hydrolysis and IAA synthesis pathway

如图11 所示，吲哚丙酮酸（indole pyruvic acid，IPyA）途径为66 株菌所共有，而其它途径均不完整或缺失。如图11 所示，15 株菌（Clade II）在醛脱氢酶（NAD+）的基因存在与缺失上与其余51 株菌（Clade I）呈“互补”状态。在IPyA 介导的IAA 合成途径中，色氨酸通过氧化转氨作用形成吲哚丙酮酸，随后在苯丙酮酸脱羧酶作用下脱羧形成吲哚乙醛，再经过醛脱氢酶（NAD+）脱氢变成IAA，而IAA 对于植物生长具有重要影响。

此外，通过blast 分析比对，所有菌株的基因组中共含有5 种编码碱性磷酸酯酶（4 个Pho D 和1 个Pho A）基因和1 种植酸酶（phy）编码基因。由图11 可知，这些基因在66 株菌中均被发现。微生物的解磷过程十分复杂，目前被广泛接受的是有机磷（organophosphorus，OPB）酶解和无机磷（inorganic phosphorus，IPB）酸解两个方式[31]。碱性磷酸酶是解磷菌参与OPB 矿化的重要酶类，分为phoA、phoD、phoX 3 种类型[32]，其中在海洋和土壤微生物中phoD 的含量更高[33]。磷酸酶通过去磷酸化或断裂有机磷物质中的磷酯键释放有效磷，植酸酶则通过释放植酸磷来增加有效磷的含量[34]。

2.4 抗性基因的挖掘与次级代谢产物基因簇

为揭示不同菌株在次级代谢产物合成以及抗生素耐药性等方面的特征，本研究利用本地化的anti SMASH 及ABRicate 注释流程对每个B.amyloliquefaciens 的基因组进行了功能挖掘，结果如图12 所示。

图12 B.amyloliquefaciens 次级代谢产物和抗性、噬菌体以及质粒基因簇分布Fig.12 Distribution of B.amyloliquefaciens secondary metabolites resistance，phage and plasmid gene clusters

如图12 所示，在所有B.amyloliquefaciens 基因组中，共鉴定出5 种抗药性基因和3 种铜离子抗性基因，菌株覆盖率高达98.48%，只有B.amyloliquefaciens B408未鉴定出抗性基因。其中YFMP、YFMO、FABL、CUER数量最多，出现频次各为66 次；其次为rphC、satAbsub、clbA、rphB、tet（L）、CSOR、COPZ，共254 个。YFMP、YFMO、CUER、CSOR、COPZ 蛋白表现出对铜的抗性；FABL 蛋白表现出对酚类化合物的抗性；rphB 和rphC蛋白对利福霉素具有耐受性；clbA 和satA-bsub 对链霉素具有耐受性，且clbA 对林可酰胺以及大环内酯类抗生素也具有耐受性，tet（L）对四环素以及多西环素具有耐药性。

噬菌体和质粒的存在可以增强细菌对环境的适应性。通过PlasmidFinder 数据库注释发现，B.amyloliquefaciens HM618、TL106、LL3、IT -45、S499、LS1 -002-014s 和MBE1283 中含有质粒，覆盖率达10.61%；PHAGE 数据库对66 株B.amyloliquefaciens 进行噬菌体基因注释，共检测出16 种噬菌体基因，覆盖率达100%，涉及参与细菌素生产或免疫、小酸溶性孢子蛋白和硫氧还蛋白的合成等功能，为B.amyloliquefaciens 的抗菌、遗传稳定、环境耐受等特性提供遗传物质基础。

antiSMASH 注释结果表明，66 株B.amyloliquefaciens 共预测到815 个基因簇，每个基因组内约包含10～15 个次级代谢基因簇。其中，与已知基因簇同源性大于60%的有441 个，同源性为100%的有354 个，主要涉及杆菌烯（bacillaene）、杆菌肽（bacillibactin）和溶杆菌素（bacilysin）等的生产（图12），这可能是B.amyloliquefaciens 作为生防菌株的主要遗传基础。

上述分析表明，作为生防菌株，B.amyloliquefaciens有望通过次级代谢产物遏制病原菌来实现抑菌效果，并通过抗药性/噬菌体/质粒增强自我保护与环境适应性；B.amyloliquefaciens 对铜离子具有较好的抗性，表明在土壤受到铜污染时，具有自我保护的能力。江春玉等[35]从重金属污染土壤中分离到Cu2+抗性B.amyloliquefaciens，可显著增加培养液中游离Cu2+的含量，从而提高植物吸收Cu2+速率，增强Cu2+污染土壤治理的效率，但相关机制尚不明确。

3 结论

本研究通过对B.amyloliquefaciens 进行泛基因组构建以及对核心基因组的功能注释，从基因组层面评估了其作为生物防治菌株的遗传潜力。ANI 以及聚类分析表明66 株B.amyloliquefaciens 分为相对独立的两个分支；COG、GO、PATHWAY 等对核心基因组基础注释表明B.amyloliquefaciens 具有较为保守的能量生产、蛋白质合成、基因表达、化合物转运等功能，这为B.amyloliquefaciens 的代谢稳定性提供了遗传基础。

在功能基因注释上，CAZy 注释结果表明B.amyloliquefaciens 具有较强的水解多糖类物质的能力。PATHWAY 和SwissProt 注释结果表明B.amyloliquefaciens 可通过IPyA 途径合成IAA 以及通过分泌植酸和碱性磷酸酶降解环境中的磷，促进植物的生长。

在66 株B.amyloliquefaciens 的次级代谢产物合成基因簇分析中预测到大量的杆菌烯和杆菌肽合成基因簇。同时抗生素耐药性分析表明，多数B.amyloliquefaciens 基因组中含有大量编码利福霉素、链霉素、林可酰胺以及大环内酯类抗生素等耐药性基因以及铜抗性基因，上述基因元件进一步增强了该菌在遏制病原菌、联合用药、自我保护和土壤改良方面的优势和潜力。总之，本文通过构建B.amyloliquefaciens 泛基因组，系统揭示了该物种内普遍性及特异性遗传特征，为将来深度利用B.amyloliquefaciens 这种微生物资源奠定了基础。