Biosynthetic gene cluster profiling predicts the positive association between antagonism and phylogeny in Bacillus

Abstract

Understanding the driving forces and intrinsic mechanisms of microbial competition is a fundamental question in microbial ecology. Despite the well-established negative correlation between exploitation competition and phylogenetic distance, the process of interference competition that is exemplified by antagonism remains controversial. Here, we studied the genus Bacillus, a commonly recognized producer of multifarious antibiotics, to explore the role of phylogenetic patterns of biosynthetic gene clusters (BGCs) in mediating the relationship between antagonism and phylogeny. Comparative genomic analysis revealed a positive association between BGC distance and phylogenetic distance. Antagonistic tests demonstrated that the inhibition phenotype positively correlated with both phylogenetic and predicted BGC distance, especially for antagonistic strains possessing abundant BGCs. Mutant-based verification showed that the antagonism was dependent on the BGCs that specifically harbored by the antagonistic strain. These findings highlight that BGC-phylogeny coherence regulates the positive correlation between congeneric antagonism and phylogenetic distance, which deepens our understanding of the driving force and intrinsic mechanism of microbial interactions.

Fig. 1. Phylogram of the tested Bacillus genomes.

The maximum likelihood (ML) phylogram of 4268 Bacillus genomes was based on the sequences of 120 ubiquitous single-copy proteins. The phylogenetic tree shows that Bacillus species can be generally clustered into the subtilis (light green circle; further includes subtilis (dark green) and pumilus (blue) subclades as shown in the branches), cereus (red), megaterium (yellow), and circulans (gray) clades. For detailed information of the species, please refer to the phylogenetic tree in Supplementary Fig. 1.

Fig. 2. Biosynthetic gene cluster (BGC) distribution is correlated with phylogeny in the genus Bacillus.

a The numbers of BGCs in the 4268 Bacillus genomes from different clades as defined by antiSMASH. In the violin plot, the centre line represents the median, violin edges show the 25th and 75th percentiles, and whiskers extend to 1.5× the interquartile range. b Hierarchal clustering among the 545 representative Bacillus genomes based on the abundance of the different biosynthesis gene cluster families (GCFs). Each column represents a GCF, which was classified through BiG-SCAPE by calculating the Jaccard index (JI), adjacency index (AI), and domain sequence similarity (DSS) of each BGC; the color bar on the top of the heatmap represents the BGC class of each GCF, where PKS includes classes of PKSother and PKSI, PKS-NRPS means PKS-NRPS Hybrids, Others includes classes of saccharides, terpene, and others. Each row represents a Bacillus genome, and the abundance of each GCF in different genomes is shown in the heatmap. The left tree was constructed based on the distribution pattern of GCFs, which showes a similar pattern to the phylogram in Fig. 1. c The correlation between the BGC and phylogenetic distance of the 545 representative Bacillus genomes (P < 2.2 × 10⁻¹⁶, R² = 0.2847). Linear model (LM) was used for the correlation analysis and adjustments were made for R² calculation; one-sided F test was applied for multiple comparisons (F-statistic: 5.669 × 10⁴ on 1 and 142,436 DF). Source data are provided as a Source Data file.

Fig. 3. Colony antagonism phenotype is positively correlated with the phylogenetic and BGC distance within Bacillus species.

a–d Inhibition of colonies of B. amyloliquefaciens ACCC19745, B. pumilus ACCC04450, B. thuringiensis YX7, and B. mobilis XL40 against Bacillus from different clades. The number below the abbreviations indicates the average 16 S rRNA gene phylogenetic distance of the target strains with the corresponding antagonistic strain. Each inhibition assay includes three biological replicates and the average is shown as points in the boxplots (n = 50, 49, 51, and 51 for a, b, c, and d, respectively) and used for the correlation analysis. In all boxplots, the centre line represents the median, box edges show the 25th and 75th percentiles, and whiskers extend to 1.5× the interquartile range; boxplots with different letters are statistically different according to the two-sided Duncan’s multiple range tests (P < 0.05), where the significance for a, b, c, and d are 9.30 × 10⁻⁵, 9.68 × 10⁻⁵, 0.007, and 0.008, respectively. e, f Correlation between the antagonism phenotype (diameter of the inhibition zone) and 16 S rRNA phylogenetic distance (e, F-statistic: 65.58 on 1 and 387 DF, P = 7.338 × 10⁻¹⁵, R² = 0.1427) or predicted BGC distance (f, F-statistic: 47.79 on 1 and 387 DF, P = 1.971 × 10⁻¹¹, R² = 0.1076) among all the tested paired strains. For strains whose genomes have not been completely sequenced, we referred to the Bacillus genomes in the NCBI database that shared the highest 16 S rRNA similarity. The error bands indicate the 95% confidence intervals. Linear model (LM) was used for the correlation analysis and adjustments were made for R² calculation; one-sided F test was applied for multiple comparisons. Source data are provided as a Source Data file.

Fig. 4. Congeneric inhibition by fermentation supernatants is positively correlated with the phylogenetic and BGC distance in Bacillus.

a Heatmap showing the antagonistic profiles of the fermentation supernatant of 17 antagonistic strains (in the left column) on the 40 target strains (in the top line). The maximum likelihood (ML) phylogenetic tree was constructed based on the 16 S rRNA sequence of the 17 antagonistic strains: subt, B. subtilis RZ30; vele, B. velezensis SQR9; halo, B. halotolerans CF7; pumi, B. pumilus ACCC04450; safe, B. safensis LY9; lich, B. licheniformis CC11; sono, B. sonorensis YX13; aqui, B. aquimaris XL39; mari, B. marisflavi XL37; gins, B. ginsengihumi ACCC05679; cere, B. cereus ACCC10263; mobi, B. mobilis XL40; wied, B. wiedmannii XL36; bing, B. bingmayongensis KF27; hori, B. horikoshii ACCC02299; flex, B. flexus DY11; mega, B. megaterium ACCC01509. Each inhibition assay includes three biological replicates and the average is shown in the heatmap and used for the correlation analysis. b, c Correlation between the antagonism phenotype (diameter of the inhibition zone) and 16 S rDNA phylogenetic distance (b, F-statistic: 115.6 on 1 and 593 DF, P < 2.2 × 10⁻¹⁶, R² = 0.1618) or the predicted BGC distance (c, F-statistic: 56.34 on 1 and 593 DF, P = 2.241 × 10⁻¹³, R² = 0.08523) among all the tested paired strains. For strains whose genomes have not been completely sequenced, we referred to the Bacillus genomes in the NCBI database that shared the highest 16 S rRNA similarity. d, e Correlation between the antagonism-phylogeny association (d, F-statistic: 10.23 on 1 and 15 DF, P = 0.005992, R² = 0.3657) or antagonism-BGC distance association (e, F-statistic: 4.57 on 1 and 15 DF, P = 0.0494, R² = 0.1824) and the (predicted) quantity of BGCs in antagonistic strains. For strains whose genomes have not been completely sequenced, we referred to the average quantity of BGCs in this species. The color of the dots represents the clade/subclade which the antagonistic strains belong to. The error bands indicate the 95% confidence intervals. Linear model (LM) was used for the correlation analysis and adjustments were made for R² calculation; one-sided F test was applied for multiple comparisons. Source data are provided as a Source Data file.

Fig. 5. Contribution of BGCs to antagonizing Bacillus species from different clades by B. velezensis SQR9.

The heatmap shows the contribution of each BGC product (on the top) to the inhibition of each target strain (on the left), which was calculated as the percentage of the decreased inhibition zone of the corresponding BGC-deficient mutant compared with wild-type. The maximum likelihood (ML) tree on the left was constructed based on the 16 S rRNA sequences of the 24 target strains: halo B. halotolerans CF7, lich B. licheniformis LY2, sono B. sonorensis YX13, alti B. altitudinis LY37, safe B. safensis LY9, pseu-A B. pseudomycoides ACCC10238, pseu-B B. pseudomycoides TZ8, albu B. albus XL388, mobi B. mobilis XL40, cere B. cereus ACCC10263, wied B. wiedmannii CF23, funi B. funiculus ACCC05674, luci B. luciferensis XL165, aqui B. aquimaris XL39, mari B. marisflavi XL37, gins B. ginsengihumi ACCC05679, coag B. coagulans ACCC10229, badi B. badius ACCC60106, fort B. fortis ACCC10219, terr B. terrae TL19, arya B. aryabhattai XL26, mega B. megaterium ACCC01509, flex B. flexus DY11, kore B. koreensis ACCC05681. The abbreviations on the top (except Sfp) represent different BGC products: Srf surfactin, Bmy bacillomycin D, Fen fengycin, Dhb bacillibactin, Bac bacilysin, Mln macrolactin, Bae bacillaene, Dfn difficidin, FA an antimicrobial fatty acid, bacillunoic acid. Specifically, Sfp (phospho-pantheinyltransferase) is not an antibiotic but is necessary for modification of the above antibiotics and ensuring their activity, except for bacilysin^,; here, the contribution of Sfp to antagonism means the relative reduction of inhibition by SQR9Δsfp compared to that of wild-type SQR9 towards different targets. “Sum” represents the overall contribution of the nine BGCs (excludes Sfp) to the antagonism against different targets. The cross represents the BGCs presence in the testing strains based on their genomic information, or for the non-sequenced strains if more than 80% of the corresponding Bacillus species genomes possessed the cluster. The curves in the right box show the antagonistic phenotype and phylogenetic distance of B. velezensis SQR9 with each of the target strains. Each inhibition assay included six biological replicates and the average contribution was shown in the heatmap. Source data are provided as a Source Data file.