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

通过平均核苷酸一致性分析、基因预测、泛基因组构建、功能注释等流程，以 66 株解淀粉芽孢杆菌（Bacillusamyloliquefaciens）的基因组为对象，进行该物种的泛基因组分析。平均核苷酸一致性分析表明，上述菌株可划为具有明显物种边界的两个独立分支；获得的泛基因组共包含 13 051 个基因，其中核心基因 2 165 个；碳水化合物活性酶（carbohydrate-active enzymes，CAZy）功能注释显示，B.amyloliquefaciens 相比 Bacillus 属内其他菌株在分解淀粉、茁-葡聚糖和木聚糖上具有较大优势；在促植物生长方面，B.amyloliquefaciens 通过吲哚丙酮酸途径合成吲哚乙酸，并通过碱性磷酸脂酶和植酸酶分解磷元素；抗生素耐药性分析发现，多数菌株含链霉素以及大环内酯类抗生素抗性基因。此外，在所有基因组中均发现杆菌烯和杆菌肽等次级代谢产物基因簇。上述研究进一步明确 B.amyloliquefaciens 基因组水平的分类学特点、核心基础代谢特征，以及菌株在碳水化合物代谢、植物促生、生物防治等方面的优势及潜力，可为开发利用该菌作参考。关键词：解淀粉芽孢杆菌；泛基因组；吲哚乙酸；解磷；耐药性；次级代谢产物Population Genetic Characteristics of Bacillus amyloliquefaciensWANG Yiting1，LI Suzhen2，MAO Zhitao3，LIAO Xiaoping3，XU Huimin1，QIU Yihua1，YAN Xuesong4，YANG Fengyan1，LI Zhenjing1，GUO Qingbin1*，LIU Huanhuan1*（1.College of Food Science and Engineering，Tianjin University of Science 驭 Technology，Tianjin 300457，China；2.Synear Foods Co.，Ltd.，Zhengzhou 450000，Henan，China；3.Tianjin Institute of IndustrialBiotechnology，Chinese Academy of Sciences，Tianjin 300308，China；4.Binzhou Municipal Marine andFisheries Development Service Center of Zhanhua District，Binzhou 256600，Shandong，China）粤遭泽贼则葬糟贼：Through processes such as average nucleotide identity analysis，gene prediction，pan-genomeconstruction，and functional annotation，the genome of 66 strains of Bacillus amyloliquefaciens was subjected topan-genome analysis.The average nucleotide identity analysis showed that the above strains could be divided intotwo distinctbranches withclear speciesboundaries.The obtainedpan-genome contained13051 genes，including2 165 core genes.Carbohydrate-active enzymes（CAZy）functional annotation showed that B.amyloliquefacienshad significant advantages over other Bacillus species in the degradation of starch，茁-glucan，and xylan.Regarding plant growth promotion，B.amyloliquefaciens synthesized indole-3-acetic acid through the indole-3-pyruvic acid pathway and utilized alkaline phosphatase and phytase to release phosphorus.Antibiotic resistanceanalysis revealed the presence of genes conferring resistance to streptomycin and macrolide antibiotics in moststrains.Additionally，secondary metabolite gene clusters such as baculene and bacitracin were identified in allgenomes.This study further elucidated the taxonomic characteristics and core metabolic features of the B.amyloliquefaciens genome，as well as the advantages and potential of strains in carbohydrate metabolism，plantgrowth promotion，and biocontrol.These findings can serve as a reference for the better development andutilizationof this bacterium.生物工程160食品研究与开发圆园23 年 9 月第 44 卷第 17 期解淀粉芽孢杆菌（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以及ABRi-cate 抗生素耐受性基因注释等分析流程，构建 B.amy-loliquefaciens 的泛基因组并进行深入分析，挖掘 B.amy原loliquefaciens 物种的普遍性遗传特征及碳水化合物的利用、植物促生、生物防治等方面的遗传特异性，以期为深度利用 B.amyloliquefaciens 这种微生物资源作参考。1材料与方法1.1材料从 NCBI 数据库（https：/www.ncbi.nlm.nih.gov/da-tasets/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 v617（https：/itol.embl.de/）绘制聚类树图。1.2.2泛基因组构建及分析采用 Prokka 注释流程13及默认参数，对所有基因组进行基因结构预测和功能注释，其中获得的 gff 格式文件作为泛基因组的输入文件。Roary 通过 blastp 对gff文件中的序列进行直系同源基因分析14，获得 B.amy原loliquefaciens 泛基因组（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 生成的泛基因运藻赠 憎燥则凿泽：Bacillus amyloliquefaciens；pan-genome；indole-3-acetic acid；phosphate-solubilizing；drug re原sistance；secondary metabolite引文格式：王依婷，李素珍，毛志涛，等.解淀粉芽孢杆菌群体遗传特征分析J.食品研究与开发，2023，44（17）：160-170.WANG Yiting，LI Suzhen，MAO Zhitao，et al.Population Genetic Characteristics of Bacillus amyloliquefaciensJ.FoodResearch and Development，2023，44（17）：160-170.生物工程161食品研究与开发圆园23 年 9 月第 44 卷第 17 期表 1B.amyloliquefaciens 基本信息Table 1Basic information of B.amyloliquefaciens组 PAV 矩阵，进行层级聚类并构建树图。1.2.3核心基因组功能注释及富集分析利用 ProCAST 平台（https：/）对 Roary 分析获得的核心基因组蛋白序列进行功能注释。该管路是由中科院天津工业生物技术研究所开发的综合蛋白质注释的无服务器网络工具，涵盖蛋白家族、域、分类、结构、基因本体和途径信息等 26 个数据库，并提供多个工具用于结果的可视化分析。1.2.4质粒、抗生素耐药性等基因的鉴定采用 ABRicate 流程（Seemann T，Abricate，https：/ CARD20、BacMet21、PHAGE、PlasmidFinder22、NCBI AMRFinderPlus23数据库，对 66 个菌株基因组的所有蛋白编码序列（codingsequence，CDS）进行功能注释。其中，CARD 与 NCBIAMRFinderPlus 为抗生素耐药性信息数据库；BacMet为抗菌杀菌剂和金属抗性基因数据库；PHAGE 为噬菌体基因数据库；PlasmidFinder 包含质粒复制子序列，可用于原始和组装测序数据的复制子序列分析。1.2.5次级代谢产物基因簇的鉴定采用 antiSMASH 本地化软件流程15，对 66 个菌株基因组注释产生的 gbk 文件进行次级代谢产物基因簇挖掘，参数采用-cb-general、-cb-subclusters、-cb-knownclusters、-asf 及-pfam2go。1.2.6碳水化合物水解酶基因的鉴定采用本地化的 CAZy 数据库对 66 个菌株基因组的所有蛋白编码序列（CDS）进行功能注释24。1.2.7解磷途径及吲哚乙酸（indoleacetic acid，IAA）合成途径相关基因的鉴定根据京都基因与基因组百科全书（KyotoEncyclope-dia of Genes and Genomes，KEGG）中 IAA 合成途径中的蛋白质酶编号25，以及解磷途径涉及的酶编号26，在UniProt 数据库中检索已校验的 fasta 蛋白序列，并通过blast 建库，对 66 个基因组蛋白序列进行比对，identity阈值设定为 40%，鉴定潜在的 IAA 合成和解磷途径。2结果与分析2.1平均核苷酸一致性分析结果ANI 被定义为两个微生物基因组同源片段之间平均的碱基相似度，是从全基因组水平评估物种间亲缘关系，尤其是近缘物种的重要工具，通常以 ANI95%作为划分物种边界的阈值16。本研究通过使用 fastANI对 66 个 B.amyloliquefaciens 基因组进行同源性评估，结果如图 1 所示。如图 1 所示，所有基因组可聚类为 2 个主要分支Clade I 与 II，其中，在分支内两两基因组之间的 ANI数值均高于 95%，而分支间的两个基因组 ANI 数值均高于 90%而低于 95%16。因此，Clade I 和 Clade II 之间已经出现较为明显的物种边界，未来可将这 2 个分支进行独立命名。2.2泛基因组功能分析利用 Prokka 对所有基因组进行基因结构预测，结果见表 1。由表 1 可知，B.amyloliquefaciens 的基因组介于Clade IClade II图 1基于 ANI 分析的 B.amyloliquefaciens 层级聚类树Fig.1Hierarchical clustering based on ANI values fromB.amyloliquefaciens genomes菌株编号登录号基因大小/Mb提交日期基因个数*菌株编号登录号基因大小/Mb提交日期基因个数*TA208bGCF_000195515.13.920114 069KC41aGCF_008824205.14.120194 128DSM7bGCF_000196735.14.020104 127WF02aGCF_013122255.14.020203 990LL3bGCF_000204275.14.020114 151205bGCF_013328835.14.020204 173XH7bGCF_000221645.13.920114 077R8-25aGCF_013341255.14.020203 946IT-45aGCF_000242855.23.920133 927T-5aGCF_014622745.14.020204 006CC178aGCF_000494835.13.920133 867INH2-4baGCF_014791945.14.020203 932LFB112aGCF_000508265.13.920133 895MOH1-5baGCF_014792065.14.020203 879KHG19aGCF_000835145.14.020153 897HM618aGCF_014854735.14.020204 010MBE1283aGCF_001483885.14.020163 958PP19bGCF_014930895.13.820203 833生物工程162食品研究与开发圆园23 年 9 月第 44 卷第 17 期注：*为 Prokka 注释流程识别出的基因个数；a 为 Clade I；b 为 Clade II。续表 1B.amyloliquefaciens 基本信息Continue table 1Basic information of B.amyloliquefaciens菌株编号登录号基因大小/Mb提交日期基因个数*菌株编号登录号基因大小/Mb提交日期基因个数*S499aGCF_001586105.13.920163 911DGL1aGCF_015999565.13.920203 867UMAF6639aGCF_001593765.14.020163 941GXU-1aGCF_017357905.14.020213 943UMAF6614aGCF_001593785.14.020163 932XJ5aGCF_017603845.14.420214 445B15aGCF_001596755.14.020163 961G10aGCF_017815555.13.920213 903RD7-7bGCF_001705195.13.720163 716Ba13bGCF_018223605.13.920213 929Y14aGCF_001874385.14.020163 947L-17aGCF_018363035.13.920213 912LM2303aGCF_001889285.14.020163 981SN16-1aGCF_018502525.14.020213 953WS-8aGCF_001922005.13.920163890GKT04aGCF_019396925.14.120214 139SRCM101267bGCF_002173635.14.120174 240JP3042aGCF_019856555.14.120213 967MT45bGCF_002209305.13.920173 94335MbGCF_019880285.13.920214 023ALB69aGCF_003149695.14.020183 987Bam1bGCF_019880325.14.020214 071ALB65aGCF_003149715.14.020183 988JP3144aGCF_019880425.14.220214 111ALB79aGCF_003149795.14.020183 874ELA1901024aGCF_020311635.14.120214 086SH-B74aGCF_003312895.14.120184 087ZF57aGCF_020971785.13.920213 884YP6bGCF_003868675.14.020184 072LS1-002-014saGCF_021233135.14.220214 273HK1bGCF_004006115.14.020194 153WF2020aGCF_022492915.14.020223 988FS1092aGCF_004421045.14.220194 240ZKY01aGCF_022559645.14.220224 180ZJU1aGCF_007362635.14.120194 008W0101aGCF_022833015.14.020223 934HbGCF_007827165.14.020194 072GL18aGCF_023115375.13.920223 868X030aGCF_007923205.14.020193 901A15.1aGCF_023238285.13.920223 890ARP23aGCF_008245105.14.020194 079CQN-2aGCF_023278425.14.120224 113DH8030aGCF_008534375.14.020193 954TL106aGCF_023920185.14.020224 005V417aGCF_008807995.13.920193 919TPS17aGCF_024138555.14.020223 928V167aGCF_008808015.13.920193 864B408aGCF_024452365.13.920223 8763.704.44 Mb，最小为 RD7-7，最大为 XJ5，GC 含量为45.50%46.50%，编码基因数目为 3 7164 445 个，平均每个基因组编码 3 997 个基因。结合 Roary 泛基因组分析流程，对上述基因组的核心基因组和泛基因组进行鉴定，结果见图 2。由图 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 易从外界获得新的基图 2核心基因数量和总基因数量随基因组数目增长的变化趋势Fig.2Variation trend of core gene number and total genenumber with the increase of genome number10 0005 000002060基因组数40保守基因数总基因数生物工程163食品研究与开发圆园23 年 9 月第 44 卷第 17 期因，但同时核心基因组具有较强的稳定性。泛基因组的缺失与存在变异（PAV）进行聚类分析，泛基因组中核心基因序列的系统发育分析结果如图 3 和图 4 所示。由图 3 与图 4 可知，此结果与 ANI 结果一致，表明66 株菌在全基因组 ANI 相似性、PAV 聚类以及核心基因组序列相似性中均可被明确划为两个分支，Clade II中的 15 株菌与 Clade I 中的 51 株菌不属于一个种，将来可作为独立物种进行研究。为进一步解析该物种所具有的核心基因功能特征，即菌株的普遍性遗传特征，本研究采用 ProCAST中集成的 21 个数据库对核心基因组蛋白序列进行了注释。其中 NR、SwissPort 和 PFAM 数据库的注释基因占核心基因数的比例均超过 90%。基于 COG 数据库的功能注释结果如图 5 所示。如图 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。在 PATHWAY 注释中（包含 KEGG pathway、Re原actome、MetaCyc 等数据库），Aminoacyl-tRNA biosyn原thesis、Mitochondrial translation elongation 以 及 Purinemetablism 等模块聚集了大量的功能基因，显示出 B.amyloliquefaciens 在转录、翻译、氨基酸代谢、离子运输、能量代谢等基础代谢方面功能保守。Clade IIClade I图 3B.amyloliquefaciens 基因存在/缺失变异（PAV）矩阵的层次聚类Fig.3Hierarchical clustering based on gene presence/absencevariation（PAV）matrix of B.amyloliquefaciens图 4基于 B.amyloliquefaciens 核心基因组序列的系统发育树Fig.4Phylogenetic tree based on core genome sequencealignment of B.amyloliquefaciens图 5核心基因组的 COG 注释结果Fig.5Main annotations for COG annotationClade IIClade I生物工程图 6核心基因组的 PATHWAY 注释结果Fig.6Main annotations for PATHWAY annotationI：lipid transport and metabolismT：signal transduction mechanismsF：nucleotide transport and metabolismO：posttranslational modification，protein turnover，chaperonesL：replication，recombination and repairH：coenzyme transport and metabolismG：carbohydrate transport and metabolismM：cell wall/membrane/envelope biogensisC：energy production and conversionP：inorganic ion transport and metabolismK：transcriptionJ：translation，ribosomal structure and biogenesisE：amino acid transport and metabolismS：function unknownR：general function prediction onlyOtherQ：secondary metabolites biosynthesis，transport and catabolismU：intracellular trafficking，secretion，and vesicular transportV：defense mechanismsD：cell cycle control，cell division，chromosome partitioning30252015105035major pathway of rRNA processing in the nucleolus and cytosol4-hydroxy-2（1H）-quinolone biosynthesisvoltage gated Potassium channelsresponse of EIF2AK4（GCN2）to amino acid deficiencypurine metabolismrespiratory electron transportregulation of pyruvate dehydrogenase（PDH）complexbranched-chain amino acid catabolismpyrimidine metabolismdetoxification of Reactive Oxygen Speciesneutrophil degranulationalanine，aspartate and glutamate metabolismaminoacyl-tRNA biosynthesismitochondrial translation elongationmitochondrial translation initiationmitochondrial translation terminationglycine，serine and threonine metabolismheme biosynthesissrp-dependent cotranslational protein targeting to membrane5-aminoimidazole ribonucleotide biosynthesis Iselenocompound metabolismtetrapyrrole biosynthesis I（from glutamate）defective MTR causes methylmalonic aciduria and homocystinuria type cblGpyrimidine deoxyribonucleotide phosphorylationl-cysteine biosynthesis II（tRNA-dependent）amino sugar and nucleotide sugar metabolismfolate transformations II（plants）terpenoid backbone biosynthesisreactomemetacyckeggglutamate and glutamine metabolismPentose phosphate pathway164食品研究与开发圆园23 年 9 月第 44 卷第 17 期如图 8 所示，在基因本体论（Gene Ontology，GO）注释中，B.amyloliquefaciens 大量核心基因聚集在 oxi原doreductase activity（7.14%）、structural constituent of ri原bosome（5.42%）、integral component ofmembrane（2.97%）等涉及氧化还原、遗传信息传递和生物膜合成等生物学过程模块；在细胞组成模块中，主要涉及 catalyticactivity（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.amy原loliquefaciens 的核心基因组中，共识别出 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。CAZy 是一个整合了碳水化合物酶活性的数据库。B.amyloliquefaciens 的核心基因组注释结果见图 7。如图 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。图 7核心基因组的 CAZy 注释结果Fig.7Main annotations for CAZy annotation图 8核心基因组的 GO 注释结果Fig.8Main annotations for GO annotation生物工程PL1S：polysaccharide lyasesAAs：auxiliary activitiesCEs：carbohydrate esterasesCBMs：carbohydrate-binding modulesGHs：glycoside hydrolasesGTs：glycosy transferases