Regulatory rewiring drives intraspecies competition in Bacillus subtilis

TL;DR

This study shows how changes in communication systems in bacteria can lead to competitive advantages, affecting microbial community dynamics.

Contribution

The study reveals that natural or selected disruptions in quorum-sensing pathways rapidly enhance competitive fitness in B. subtilis.

Findings

Disruption of the RapP-PhrP module in B. subtilis leads to faster growth and competitive exclusion.

Recurrent mutations in the sensor kinase comP confer a growth-linked competitive advantage in wild isolates.

A type VII secretion system toxin contributes to competitive fitness differences among isogenic strains.

Abstract

Intraspecies interactions shape microbial community structure and evolution, yet the mechanisms determining competitive outcomes among closely related strains remain unclear. The soil bacterium Bacillus subtilis is a model for microbial social interactions, where quorum-sensing systems regulate cooperation and antagonism. Here, we take a multifaceted approach to dissect the role of quorum-sensing regulation in competitive fitness. Isolate NCIB 3610 carries a signal unresponsive RapP-PhrP module that alters quorum-sensing control and promotes faster growth. Modelling and mutant analysis demonstrate that the small differences in growth rate conferred by RapP-PhrP3610 are sufficient to drive competitive exclusion. The importance of quorum sensing control is further exemplified by experimental evolution of distinct wild isolates, which revealed recurrent mutations in the sensor kinase comP, which phenocopy complete comP or comA deletions and confer a growth-linked competitive advantage. Key quorum sensing mechanisms are abandoned even in structured microbial communities, where it might be expected that communal traits are favoured. Furthermore, a phylogenomic survey of 370 B. subtilis genomes identified disruptive comP mutations in ~16% of isolates. However, growth rate alone does not explain all interaction outcomes as even isogenic strains with equivalent doubling times differ in competitiveness. Transcriptomic profiling and validation experiments implicated a type VII secretion system toxin as an additional effector. These findings reveal that disruption of quorum-sensing pathways, whether naturally or through selection, provides a rapid route to competitive advantage, highlighting a fundamental trade-off between communal signalling and individual fitness in microbial populations.

Microbial competition and cooperation are key in shaping the structure and evolution of microbial communities. Our study on Bacillus subtilis, a model for microbial social interactions, reveals how alterations in cell-cell communication can enhance competitive fitness. We show, through a combination of modelling, mutational analysis, and experimental evolution, that certain strains of B. subtilis gain competitive advantage by disrupting their quorum sensing pathways, which leads to faster growth and enhanced competitiveness. Such mutations are prevalent in ~16% of analysed genomes, underscoring their widespread evolutionary benefit under competitive conditions. Our research highlights the delicate balance between individual success and community cooperation in microbes, providing insights into microbial behaviour and potential applications in modifying microbial ecosystems.

Figures

Deletion of rapP in NCIB 3610 decreases competitive fitness in mixed biofilms.(A) Schematic depicting the interconnections between the ComQXPA quorum sensing system and the RapP-PhrP quorum sensing module in NCIB 3610 (left) and other wild type isolates that contain a signal responsive RapP (right). ComX and PhrP are produced as pre-peptides that are exported and processed into the extracellular environment. At high cell density ComX binds ComP and triggers phosphorylation of the cognate response regulator ComA, while PhrP is imported to the cytoplasm. In NCIB 3610 RapP carries a single amino acid substitution (N236T) which renders it insensitive to PhrP levels and thus the level of phosphorylated ComA remains low due to the dephosphorylation activity of RapP. In either a rapP deletion strain (not depicted) or in isolates where RapP is responsive to PhrP the level of phosphorylated ComA can rise in line with increases in cell density. (B) Relationship between competitive outcome and relative growth rates of in silico strains B1 and B2. (C-D) Model simulations with equal growth rates (C) or with B2 having a growth advantage (relative growth rate = 1.3) (D). Black areas show the computational domain (Ω). Scale bars represent 28 nondimensional units. Strains are shown individually (B1, green; B2, magenta) and merged (grey = overlap). (E) Colony biofilms of NCIB 3610 and NCIB 3610 ΔrapP grown for 48 h at 30°C. Scale bars, 0.5 cm. The ΔrapP strain expressed mTagBFP (NRS7280). (F) Doubling time (min) during exponential growth at 30°C with shaking. p ≤ 0.05, unpaired t test with Welch’s correction. Each point = one sample; error bars = SD. (G) Representative dual-strain colony biofilm after 48 h at 30°C. Fluorescence signals from mTagBFP (magenta) and GFP (green). Scale bar, 10 mm.1 (H) Quantification of competitive outcomes as % challenger NCIB 3610 (mTagBFP+) remaining providing a measure of relative strain density. Each point = one colony biofilm; error bars = SD. p < 0.0001, unpaired t test with Welch’s correction.

The impact of RapP on NCIB 3610 competitive fitness is isolate dependent.(A) Quantification of pairwise biofilm competitions shown as % challenger remaining (relative strain density). Challengers were NCIB 3610 (mTagBFP +, magenta; data from [14]) or NCIB 3610 ΔrapP (mTagBFP +, green). Partner strains are indicated on the y-axis. Each point = one biofilm; error bars = SD. p < 0.0001, unpaired t test. (B) RapP effect on competition outcomes, calculated as (% remaining after ΔrapP co-culture – % remaining after NCIB 3610 co-culture). The definition of a RapP effect is designated in the legend and was > 5% change in the interaction outcome and a p value of <0.05. (C) Maximum-likelihood phylogeny based on gyrA, rpoB, dnaJ, and recA, with bootstrap values shown as branch labels (adapted from [14]). Strains showing a RapP effect are marked as in (B). (D) Doubling time (min) of each strain during exponential growth at 30 °C with shaking. a p-value of ≤0.05 (*); ≤ 0.01 (); ≤ 0.001 (*) or ≤0.0001 () vs. NCIB 3610 (one-way ANOVA, Dunnett’s test). Each point = one sample; error bars = SD. (E) Presence/absence of rap genes in NCIB 3610, NRS6105, NRS6145, and NRS6202. (F) Representative colony biofilms grown for 48 h on biofilm-inducing medium at 30°C. Scale bars, 10 mm. (G) Quantification of % NCIB 3610 (gfp+) remaining after co-culture with the indicated strain for 48 h. “–ve” = NCIB 3610 alone; “+ve” = NCIB 3610 (gfp+) (NRS6942) alone. RapP-PhrP presence (+) or absence (–) in the chromosome is shown in the legend. Each point = one biofilm; error bars = SD. p thresholds as in (D).

Experimental evolution in dual-isolate colony biofilms.(A) Prevalence of NCIB 3610 (gfp+) after 24 h in dual biofilms, replotted from [14]. Data shown as % challenger remaining (relative strain density). (B) Representative confocal image of NCIB 3610 (magenta) with NRS6105 (green) after 24 h. Scale bar, 500 μm. (C) Colony-forming units (CFU) per biofilm after 24 h co-culture at 30 °C. p-value of ≤0.05 (*); ≤ 0.01 (); ≤ 0.001 (*) or ≤0.0001 (), one-way ANOVA with Tukey’s test. (D) Experimental design schematic: soil isolates (NRS6105, NRS6145, NRS6202) co-cultured with NCIB 3610; evolved clones re-isolated on spectinomycin plates. (E) Colony biofilms of ancestral and evolved isolates grown alone (“single”) or mixed with NCIB 3610 (gfp+) for 24 h. NCIB 3610 = green; soil isolates = magenta. P1–P6 denote passages; strain background and lineage indicated below. Scale bars, 3 mm. (F–H) Quantification of the prevalence of NCIB 3610 after 48 h in co-culture with ancestral and evolved lineages as measured by relative strain density: (F) NRS6105, (G) NRS6145, (H) NRS6202. Statistics as in (C). (I) Schematic of mutations identified in evolved isolates. Dashed lines and numbers indicate genomic positions; strain backgrounds are indicated.

Distribution of mutations in the B. subtilis comP coding region.(A) 370 genomes of B. subtilis with comP status indicated as intact or disrupted. (B) The position of the comP disruption with respect to the ComP amino acid sequence is indicated with the calculated length of the variant protein produced (using NCIB 3610 ComP sequence as the reference) indicated by a black circle. (C) The string of eight thymidines found in wild type comP that form a mutation hotspot is indicated. This sequence is not found in all pherotypes of comQXPA.

ComQXPA influences the outcome of intraspecies biofilm interactions.(A) Quantification of NCIB 3610 (gfp+) prevalence (relative density) after 48 h in dual colony biofilms with partner strains (mTagBFP+). % challenger remaining was measured by image analysis. p < 0.0001, one-way ANOVA with Sidak’s test. (B) Representative colony biofilms corresponding to (A). Scale bars, 10 mm (left) or 3 mm (zoom). (C) Confocal images of NCIB 3610 (gfp+) vs. NRS6105 or NRS6105 ΔcomP (mTagBFP+) at 24 h and 48 h. Scale bar, 100 μm. (D) Prevalence of NRS6105 (gfp+) in dual biofilms with indicated partner strains after 48 h. (E) Representative images for (D). Scale bars, 10 mm (left) or 3 mm (zoom). (F–H) Outcomes of dual biofilms with additional isolates. % challenger remaining was measured by image analysis. Significance determined by Welch’s t test (F–G) or ANOVA with Tukey’s test (H). (I) Doubling times of NCIB 3610 and indicated isolates during exponential growth at 30 °C. p ≤ 0.05, Welch’s t test. (J) Interaction outcomes with NCIB 3610 (mTagBFP+) as challenger and the indicated strains (gfp+). p < 0.0001, ANOVA with Tukey’s test. (K) Doubling times of NRS6942, NRS7279, and NRS7771 at 30 °C. p ≤ 0.01, ANOVA with Tukey’s test. Note - Data for NCIB 3610 in Fig 5A, 5J and 5H & NCIB 3610 ΔcomA in Fig 5J and 5H are the same dataset.

The Type VII section system toxin YwqJ is a key determinant.(A) Schematic of the ywq operon based on genomic information extracted from Subtiwiki [30]. (B) RNA sequencing log2FC values for the members of the ywq operon. (C) Heatmap depicting the % protein sequence identity of the top matched LXG containing protein from the suite of B. subtilis isolates (x-axis) relative to the indicated LXG protein encoded by NCIB 3610 (y-axis). The competition outcome of the different wild isolates when cocultured with NCIB 3610 ΔrapP which has an equivalent doubling time is shown (data from Fig 2A). (D) Interaction outcomes with ΔrapP ΔcomA (gfp+) (NRS7771) as challenger and the indicated strains (mTagBFP+). p < 0.0001, significance determined by Welch’s t test.

Full text

Introduction

Microbial life is inherently social. In nutrient-limited environments, microbes compete for survival but also cooperate through behaviours such as public goods sharing, division of labour, and collective defence actions. These interactions profoundly shape community structure and consequently ecosystem function, influencing not only bacterial populations but also their interactions with fungi, plants, and animals [1–3]. Competition is particularly intense among closely related strains that share ecological niches, where outcomes depend on both cooperative investments and antagonistic strategies.

The Gram-positive bacterium Bacillus subtilis is an established model for dissecting microbial sociality [4]. Within clonal biofilms, genetically identical cells differentiate into specialised subpopulations, coordinating cooperative matrix production [5] while also deploying antagonistic toxins to regulate growth and sporulation [6–9]. In natural environments, however, B. subtilis strains frequently encounter conspecific rivals. Such encounters can be shaped by kin discrimination, where strains cooperate with close relatives but compete with more distant lineages [10–12]. Although molecules such as polymorphic toxins and secondary metabolites have been implicated in competitive outcomes [10–13], strain dominance cannot always be predicted from genomic content alone [14].

Central to B. subtilis social regulation are quorum-sensing systems that modulate cooperative and competitive behaviours. The ComQXPA system (Fig 1A) exhibits extensive polymorphism, partitioning isolates into distinct “pherotypes” that can communicate privately [15,16]. Parallel Rap-Phr modules further regulate ComA activity antagonistically [17,18], and diversify social traits displayed by different isolates [19]. Of intense study is strain NCIB 3610, which harbours an atypical RapP-PhrP system that is encoded on plasmid pBS32 [20,21] (Fig 1A). NCIB 3610 RapP carries threonine at position 236 which renders it insensitive to PhrP levels and thus the level of phosphorylated ComA remains low at high cell density due to the dephosphorylation activity of RapP. In either a rapP deletion strain or in isolates where RapP is responsive to PhrP the level of phosphorylated ComA can rise in line with increases in cell density [10,20,21]. The overall consequence is that the RapP status alters the balance of cooperative regulation and tips the scales towards enhanced growth [22].

*Deletion of rapP in NCIB 3610 decreases competitive fitness in mixed biofilms.(A) Schematic depicting the interconnections between the ComQXPA quorum sensing system and the RapP-PhrP quorum sensing module in NCIB 3610 (left) and other wild type isolates that contain a signal responsive RapP (right). ComX and PhrP are produced as pre-peptides that are exported and processed into the extracellular environment. At high cell density ComX binds ComP and triggers phosphorylation of the cognate response regulator ComA, while PhrP is imported to the cytoplasm. In NCIB 3610 RapP carries a single amino acid substitution (N236T) which renders it insensitive to PhrP levels and thus the level of phosphorylated ComA remains low due to the dephosphorylation activity of RapP. In either a rapP deletion strain (not depicted) or in isolates where RapP is responsive to PhrP the level of phosphorylated ComA can rise in line with increases in cell density. (B) Relationship between competitive outcome and relative growth rates of in silico strains B1 and B2. (C-D) Model simulations with equal growth rates (C) or with B2 having a growth advantage (relative growth rate = 1.3) (D). Black areas show the computational domain (Ω). Scale bars represent 28 nondimensional units. Strains are shown individually (B1, green; B2, magenta) and merged (grey = overlap). (E) Colony biofilms of NCIB 3610 and NCIB 3610 ΔrapP grown for 48 h at 30°C. Scale bars, 0.5 cm. The ΔrapP strain expressed mTagBFP (NRS7280). (F) Doubling time (min) during exponential growth at 30°C with shaking. p ≤ 0.05, unpaired t test with Welch’s correction. Each point = one sample; error bars = SD. (G) Representative dual-strain colony biofilm after 48 h at 30°C. Fluorescence signals from mTagBFP (magenta) and GFP (green). Scale bar, 10 mm.1 (H) Quantification of competitive outcomes as % challenger NCIB 3610 (mTagBFP+) remaining providing a measure of relative strain density. Each point = one colony biofilm; error bars = SD. ***p < 0.0001, unpaired t test with Welch’s correction.

The mechanisms that confer competitive superiority to NCIB 3610, and whether susceptible strains can evolve counter-adaptations, remain poorly understood. Here, we combine experimental evolution, genetic reconstruction, transcriptomics, and comparative genomics to dissect the role of quorum-sensing regulation in competitive fitness. Our results reveal that the presence of the signal unresponsive RapP-PhrP module, or disruption of the ComQXPA system, underpin the strain’s ability to dominate diverse environmental isolates in colony biofilm competitions. Furthermore, disruption of the ComQXPA systems provides a recurrent evolutionary strategy both in laboratory and natural populations that uncouples cells from costly regulation, accelerates growth, and enhances competitiveness in environments dominated by conspecific rivals. However, intraspecies interactions cannot be reduced to growth rate alone: here we demonstrate that isogenic mutants with equivalent growth diverge in competitive ability. Transcriptomic profiling implicates toxins, including those secreted by type VII secretion systems, as additional drivers of these outcomes.

Results

Signal unresponsive RapP of B. subtilis NCIB 3610 impacts pairwise intraspecies competition outcomes

We previously demonstrated that B. subtilis strain NCIB 3610 dominates pairwise interactions when co-cultured with an array of B. subtilis isolates in a dual isolate colony biofilm setting [14]. B. subtilis quorum sensing is regulated by ComQXPA and a suite of Rap-Phr accessory systems (Fig 1A). The allele of RapP encoded by NCIB 3610 is associated with increased fitness in an otherwise isogenic context [23] and linked with an accelerated growth rate [22,23]. We therefore questioned whether it is NCIB 3610 RapP-PhrP (RapP-PhrP^3610^) that drives the successful outcome for NCIB 3610 during intraspecies interactions. We first tested this hypothesis by adapting an existing mathematical model [24] which describes the spatio-temporal dynamics of two isolates in a single growing biofilm. For simplicity, the partial differential equations model reduces the dynamics to local growth through a logistic growth function and spatial spread through negative density-dependent diffusion (see Methods). Our in-silico modelling revealed that a seemingly small difference in doubling times of isolates co-cultured in a single biofilm manifested in a large impact on the competition outcome when all else is assumed to be equal. The in-silico isolate with the shorter doubling time rapidly dominated the pairwise interaction (Fig 1B-1D; [24]).

To explore the role of RapP^3610^ in mediating the competitiveness of NCIB 3610, we constructed an NCIB 3610 ΔrapP mutant and noted a mucoid colony biofilm morphology (Fig 1E) [21], reminiscent of that formed by environmental isolates of B. subtilis (S1A Fig). Measuring the doubling time under planktonic growth conditions during exponential phase revealed that NCIB 3610 was on average ~1.4-fold faster compared to the otherwise isogenic NCIB 3610 ΔrapP mutant (Fig 1F). Consistent with growth contributing to the competition outcome, NCIB 3610 readily outcompeted NCIB 3610 ΔrapP in a dual isolate colony biofilm competition setting (Fig 1G-1H). In these experiments, we employed NCIB 3610 ΔrapP expressing mTagBFP (NRS7280) and co-cultured it with NCIB 3610 constitutively expressing GFP (NRS6942). As control experiments, we co-cultured isogenic variants of NCIB 3610 and NCIB 3610 ΔrapP which differed in whether they produced GFP or mTagBFP [25] as the “reporter”. As anticipated, each of the control strain combinations showed coexistence in the mature population (Figs 1G-1H, S1B, and S1C).

The impact of RapP on the competition outcome for an array of B. subtilis isolates

To delve deeper into the role of RapP in mediating the competitiveness of NCIB 3610 we examined pairwise interactions between the NCIB 3610 ΔrapP mutant (constitutively producing mTagBFP (NRS7280)) and a suite of different wild-type B. subtilis isolates, modified to constitutively produce GFP, that we had previously used [14]. We overlayed the data obtained with our previously published results that used NCIB 3610 (mTagBFP+) as the challenger strain [14] (Fig 2A) and calculated the impact of RapP on the outcome of each pairwise interaction (Fig 2B). Our findings demonstrate a strain-dependent contribution of RapP to NCIB 3610’s highly competitive phenotype. Some strains that could co-exist with NCIB 3610, albeit in minimal proportions, became more prevalent in the mature colony biofilm when cocultured with NCIB 3610 ΔrapP. These strains included NRS6202, NRS6105, NRS6145, NRS6107, NRS6116, NRS6118, NRS6121, and NRS6190. Other strains that were entirely outcompeted by NCIB 3610 remained virtually undetectable within the population (e.g., NRS6103, NRS6132, NRS6181, and NRS6187). By aligning the interaction outcome with strain phylogeny, we established that the RapP-mediated contribution to NCIB 3610 competitiveness was broadly associated with the relatedness of the soil isolate to NCIB 3610 (Fig 2C). The more closely related the co-cultured strain was to NCIB 3610 in terms of phylogeny, the greater the extent of RapP’s impact compared with more distantly related isolates.

The impact of RapP on NCIB 3610 competitive fitness is isolate dependent.(A) Quantification of pairwise biofilm competitions shown as % challenger remaining (relative strain density). Challengers were NCIB 3610 (mTagBFP +, magenta; data from [14]) or NCIB 3610 ΔrapP (mTagBFP +, green). Partner strains are indicated on the y-axis. Each point = one biofilm; error bars = SD. p < 0.0001, unpaired t test. (B) RapP effect on competition outcomes, calculated as (% remaining after ΔrapP co-culture – % remaining after NCIB 3610 co-culture). The definition of a RapP effect is designated in the legend and was > 5% change in the interaction outcome and a p value of <0.05. (C) Maximum-likelihood phylogeny based on gyrA, rpoB, dnaJ, and recA, with bootstrap values shown as branch labels (adapted from [14]). Strains showing a RapP effect are marked as in (B). (D) Doubling time (min) of each strain during exponential growth at 30 °C with shaking. a p-value of ≤0.05 (); ≤ 0.01 (); ≤ 0.001 () or ≤0.0001 (**) vs. NCIB 3610 (one-way ANOVA, Dunnett’s test). Each point = one sample; error bars = SD. (E) Presence/absence of rap genes in NCIB 3610, NRS6105, NRS6145, and NRS6202. (F) Representative colony biofilms grown for 48 h on biofilm-inducing medium at 30°C. Scale bars, 10 mm. (G) Quantification of % NCIB 3610 (gfp+) remaining after co-culture with the indicated strain for 48 h. “–ve” = NCIB 3610 alone; “+ve” = NCIB 3610 (gfp+) (NRS6942) alone. RapP-PhrP presence (+) or absence (–) in the chromosome is shown in the legend. Each point = one biofilm; error bars = SD. p thresholds as in (D).

To explore if the variable impact of rapP on the competition outcome was underpinned by differences in growth, we determined the doubling time during the exponential phase of planktonic growth for a selection of the strains (Fig 2D). We chose three strains that showed a competitive advantage against NCIB 3610 ΔrapP (NRS6105, NRS6145, NRS6202), and three strains (NRS6103, NRS6099, NRS6186) that showed no change in the interaction outcome. We determined that each strain had a significantly slower doubling time compared to NCIB 3610 but comparable to that of NCIB 3610 ΔrapP. Therefore, growth rate is not the sole determinant of the interaction outcome. This conclusion is supported by the outcome of the pairwise interaction screen when the closely phylogenetically related isolate NRS6202 replaced NCIB 3610/ NCIB 3610 ΔrapP as the challenger strain. The same pattern of interaction outcomes against the other isolates was measured as when NCIB 3610ΔrapP was the challenger (compare Figs 2A and S2A).

Rap-Phr protein pairs are found as discrete modules that are found in variable number in different wild isolates of B. subtilis [19]. This allows for plasticity in quorum sensing responses. Analysis of 370 complete Bacillus subtilis genomes (see methods) detected a complete rapP coding region and downstream phrP in 38 strains (41 genome sequences, accounting for duplicate sequencing of the same strain, e.g., NCIB 3610), and an additional 3 strains (3 genome sequences) carried a copy of rapP with a 5’ truncation. 23 of the strains harbouring an intact rapP-phrP, including NCIB 3610, carry rapP on pBS32 (S2B Fig), the remainder of the modules are integrated into the chromosome. Seven of the 38 rapP-encoding strains harbour the NCIB 3610 variant of the rapP coding region that is PhrP unresponsive. Metabolic activity can be rewired when signal unresponsive RapP^3610^ is transferred to a closely related domesticated strain of B. subtilis [22]. We sought to understand if introduction of the RapP-PhrP^3610^ coding region into more distantly related wild-type isolates of B. subtilis would also accelerate growth and correspondingly enhance competitiveness against NCIB 3610. In essence, this represents the reciprocal of deleting rapP from NCIB 3610 and observing a drop in competitiveness. We used NRS6105, NRS6145, and NRS6202, of which, NRS6145 encodes RapP-PhrP (the signal-responsive variant) (see S1 File and Fig 2E). We introduced the rapP-phrP^3610^ coding region under the control of the native NCIB 3610 rapP promoter at the ectopic amyE location on the chromosome. Introduction of RapP-PhrP^3610^ conferred an “NCIB 3610-like” colony biofilm morphology in each isolate (Fig 2F) and a decrease in the doubling time measured in planktonic cultures during exponential phase (Fig 2D). In addition, there was an increase in competitive fitness when co-cultured with NCIB 3610 (gfp+) compared to the respective parental strains. For these dual isolate co-cultures, the proportion of NCIB 3610 (gfp+) remaining in the mature population was measured using flow cytometry (Figs 2G and S2C). Collectively, these data demonstrate that RapP-PhrP^3610^ is dominant over the native RapP-PhrP encoded by NRS6145, and that introduction of RapP-PhrP^3610^ increases an isolate’s ability to endure an encounter with NCIB 3610.

Experimental evolution of B. subtilis soil isolates

Given the relative rarity of rapP (S2B Fig and S1 File) and the knowledge that the genetic landscape of strains can adapt under stress, we aimed to identify other spontaneous genetic adaptations that could enhance the ability of B. subtilis strains (normally outcompeted by NCIB 3610) to withstand such encounters. We selected NRS6105, NRS6145, and NRS6202 as each of these strains is largely outcompeted by NCIB 3610 within a pairwise interaction in a dual isolate colony biofilm competition setting (Fig 3A), but the isolates can persist at a low level (for example, Fig 3B), and each can be re-isolated as single colonies from the mixed population (Fig 3C). We used derivatives of the NRS6105, NRS6145, and NRS6202 that were engineered to be spectinomycin resistant and constitutively produce mTagBFP (the derived isolates are called NRS6936, NRS6951, and NRS7201 respectively and will be referred to for simplicity by the parent strain name (mTagBFP+)). We used a variant of NCIB 3610 that constitutively expressed GFP and was resistant to kanamycin as the challenger strain (NRS1473, NCIB 3610 (gfp+)) (S1 Table). We conducted 6 passages in which NRS6105 (mTagBFP+), NRS6145 (mTagBFP+), and NRS6202 (mTagBFP+) were co-incubated with NCIB 3610 (gfp+) in mixed isolate colony biofilms for 24h (Fig 3D). NRS6105 (mTagBFP+), NRS6145 (mTagBFP+), and NRS6202 (mTagBFP+) were selectively recovered each time using spectinomycin antibiotic resistance selection. We conducted this experiment concurrently for 10 lineages of each soil isolate (three of the lineages of NRS6202 (mTagBFP+) were stopped prematurely) (Figs 3E and S3A-S3C).

Experimental evolution in dual-isolate colony biofilms.(A) Prevalence of NCIB 3610 (gfp+) after 24 h in dual biofilms, replotted from [14]. Data shown as % challenger remaining (relative strain density). (B) Representative confocal image of NCIB 3610 (magenta) with NRS6105 (green) after 24 h. Scale bar, 500 μm. (C) Colony-forming units (CFU) per biofilm after 24 h co-culture at 30 °C. p-value of ≤0.05 (); ≤ 0.01 (); ≤ 0.001 () or ≤0.0001 (**), one-way ANOVA with Tukey’s test. (D) Experimental design schematic: soil isolates (NRS6105, NRS6145, NRS6202) co-cultured with NCIB 3610; evolved clones re-isolated on spectinomycin plates. (E) Colony biofilms of ancestral and evolved isolates grown alone (“single”) or mixed with NCIB 3610 (gfp+) for 24 h. NCIB 3610 = green; soil isolates = magenta. P1–P6 denote passages; strain background and lineage indicated below. Scale bars, 3 mm. (F–H) Quantification of the prevalence of NCIB 3610 after 48 h in co-culture with ancestral and evolved lineages as measured by relative strain density: (F) NRS6105, (G) NRS6145, (H) NRS6202. Statistics as in (C). (I) Schematic of mutations identified in evolved isolates. Dashed lines and numbers indicate genomic positions; strain backgrounds are indicated.

Of the 27 lineages generated, 6 lineages became better able to co-exist with NCIB 3610 (gfp+). After validation, the number of lineages was reduced to four, which comprised two lineages derived from the NRS6145 (mTagBFP+) parental strain, one lineage from NRS6105 (mTagBFP+), and one from NRS6202 (mTagBFP+) (Figs 3E and S3A-S3C). In each case, the change in the ability to co-exist with NCIB 3610 (gfp+) became apparent from the 4^th^ passage (Fig 3E-3H). For each of the evolved strains that displayed an enhanced ability to co-exist, the mixed isolate colony biofilm that formed showed a distinct spatial distribution of the two strains within the overall architecture when imaged (Fig 3E-3H). The evolved soil isolates dominated the middle of the biofilm while NCIB 3610 (gfp+) occupied the periphery, essentially encircling the evolved isolates. In addition to the change in the phenotype of the mixed biofilms, the single isolate colony biofilms of the evolved variants showed a different phenotype to their parental isolates, such that they appeared less mucoid and occupied a reduced footprint (S3D Fig).

We acquired short-read genome sequencing data for each of the evolved isolates at passage 3 and 4. Sequence analysis revealed that in all four lineages, there was a single genetic change present between the passages (S2 File). For three out of four evolved isolates, the mutations appeared within the comP coding region, a core component of the B. subtilis quorum sensing system (recall Fig 1A). There was one comP mutant from each of the strain backgrounds subjected to evolution. The comP coding region of the evolved isolate of NRS6202 contained a single adenine insertion at chromosome position 3,298,556 bp which resulted in a truncation of the comP coding region. In strain NRS6105 the nucleotides that encode amino acids 633–649 (inclusive) of ComP were excised. These 17 amino acids reside in a helical structure in proximity to the catalytic and ATP binding domain, specifically in helix α8 (S3E-S3F Fig). Although this deletion does not directly interfere with ATP binding, based on structural modelling (S3G-S3H Fig), and the side chain of His691 remains in proximity to the ATP binding site, it likely alters the overall conformation of the catalytic domain and its spatial relationship within the ComP dimer, potentially affecting the interface with the cognate response regulator and substrate ComA. Therefore, ComP_Δ633–649_ variant is expected to be defective of phosphorylating its cognate response regulator ComA. In the evolved isolate of NRS6145 (L1P4) the chromosome contained a single nucleotide mutation at position 3,101,476 bp which yielded an amino acid substitution in ComP - E510K. We again generated structural models for both the ComP^NRS6145^ (WT) and ComP^NRS6145_E510K^ mutant sequences using AlphaFold3 [26] (see methods). In the WT model, a potential hydrogen bond is observed between the side chain of E510 and R331, with both residues being highly conserved among ComP homologues (S3 File). Additionally, E394, located on an adjacent α-helix, is positioned within hydrogen-bonding distance to R331, suggesting a stabilising three-residue interaction network that may contribute to the relative packing of these helices. In contrast, inspection of the E510K mutant models indicates a reorganisation of a hydrogen-bonding network. The lysine side chain at position 510 does not maintain the interaction geometry with R331 as seen for wild type ComP^NRS6145^. Instead, K510 can form a hydrogen bond with nearby residues (e.g., E394), while R331 appears displaced, and its distance to residue 510 exceeds typical hydrogen bond cutoffs. This rearrangement suggests that the local tertiary structure of the three helices may be altered, effectively “locking” the region into an alternative conformation compared to the WT protein. Taken together, the modelling results point to a rewiring of hydrogen-bonding interactions between residues E/K510, R331, and E394, consistent with the evolved isolate being more competitive due to presumed inactivation of ComP. The last mutant analysed was NRS6145-L6P4 which contained a 4,542 bp deletion in the chromosome. The region deleted covers the majority of the rhgR1 and opuBC genes along with all the genes between them (rghR2, yvaP, yvaQ, and opuBD) (Fig 3I and S2 File).

Mutations within comP are detected in the genomes of B. subtilis isolates

Given the propensity of comP to acquire mutations during our evolution experiments, we asked whether comP mutations were present in natural populations of B. subtilis. We analysed 370 high-quality genomes (S2 Table), placed them within a phylogenetic framework (Fig 4A), and extracted comQXPA sequences with a focus on comP. As the comQXP region is highly polymorphic [27], we screened specifically for (i) internal nonsense mutations causing premature truncation, (ii) frameshift-inducing indels, (iii) in-frame deletions within the conserved 3′ region of comP, and (iv) transposon insertions [28]. Using these criteria, we identified 60 genomes carrying disrupted comP alleles, representing at least 16% of the dataset. This figure is likely an underestimate, as point mutations disrupting ComP function were not included. To validate this classification, we sourced seven strains from stock centres: four predicted to carry an intact comP allele (GCA_029537135, DSM 13109; GCA_037914925, DSM 10; GCA_024507895, BGSC 10A5; GCA_024508175, BGSC B98af) and three predicted to carry disrupted alleles (GCA_029536975, DSM 5611; GCA_029537155, DSM 1090; GCA_004101945, ATCC 11774/NCTC 8236). Short-read sequencing confirmed the predicted comP status in all cases, supporting the accuracy of our bioinformatic screen (S2 File).

Distribution of mutations in the B. subtilis comP coding region.(A) 370 genomes of B. subtilis with comP status indicated as intact or disrupted. (B) The position of the comP disruption with respect to the ComP amino acid sequence is indicated with the calculated length of the variant protein produced (using NCIB 3610 ComP sequence as the reference) indicated by a black circle. (C) The string of eight thymidines found in wild type comP that form a mutation hotspot is indicated. This sequence is not found in all pherotypes of comQXPA.

Computational analysis of the 60 comP mutations showed that disruptive changes occurred throughout the coding sequence (Fig 4B) and were distributed broadly across the phylogeny. One notable cluster of closely related strains carried mutations within a poly-thymidine tract (eight consecutive T residues; Fig 4C), likely reflecting a mutational hotspot caused by slipped-strand mispairing during replication. These 18 strains were isolated from a wide range of environments between 2004 and 2023, including rhizosphere soil (wheat, cucumber), plant tissues (e.g., Passiflora edulis, tobacco leaves), aquatic and anthropogenic sources (river water, aquaculture ponds, sewage sludge, human faeces), and marine invertebrates (gills of Teredo navalis). Most isolates originated from multiple provinces in China (Guangxi, Henan, Jiangsu, Guangzhou, Shanxi, Changzhou, Lanzhou, Baoding, Anshun), with additional isolates from agricultural soil in Canada and from South Korea. Thus, this cluster of isolates reflects both ecological and geographic diversity while retaining phylogenetic similarity.

Finally, we assigned isolates to ComQXPA pherotypes based on ComQ sequence (S4 Fig; [15]). Strains within the comP-mutant cluster all belonged to the NCIB 3610 ComQ clade, indicating a shared pherotype. Sequence divergence in the 5′ coding region of other pherotypes prevents the formation of an equivalent poly-thymidine tract, explaining the apparent restriction of this hotspot to the NCIB 3610 clade (S3 File). Collectively, our analysis demonstrates that comP mutations are widespread across natural and anthropogenic environments, occur across phylogenetic groups, and can arise convergently within specific pherotypes.

Impacts associated with mutation of the ComQXPA quorum sensing system

To validate the role of ComP in shaping intra-species interactions we reconstructed the comPΔ633–649 deletion (Fig 5A-5B) and a complete in-frame comP deletion in isolate NRS6105. Each comP mutant strain (mTagBFP+) was challenged against NCIB 3610 (GFP+) in dual-isolate colony biofilms. The outcomes of these competitions resembled those observed during our experimental evolution experiments (compare Figs 3 with 5A–5C). Both comP mutants were more prevalent than the parental NRS6105 strain, although NCIB 3610 maintained a competitive advantage. These results show that deletion of residues 633–649 produces a phenotype equivalent to a full comP deletion (Fig 5A–5B), revealing that, consistent with the in silico structural analysis (S3G-S3H Fig, ComP is inactivated.

*ComQXPA influences the outcome of intraspecies biofilm interactions.(A) Quantification of NCIB 3610 (gfp+) prevalence (relative density) after 48 h in dual colony biofilms with partner strains (mTagBFP+). % challenger remaining was measured by image analysis. **p < 0.0001, one-way ANOVA with Sidak’s test. (B) Representative colony biofilms corresponding to (A). Scale bars, 10 mm (left) or 3 mm (zoom). (C) Confocal images of NCIB 3610 (gfp+) vs. NRS6105 or NRS6105 ΔcomP (mTagBFP+) at 24 h and 48 h. Scale bar, 100 μm. (D) Prevalence of NRS6105 (gfp+) in dual biofilms with indicated partner strains after 48 h. (E) Representative images for (D). Scale bars, 10 mm (left) or 3 mm (zoom). (F–H) Outcomes of dual biofilms with additional isolates. % challenger remaining was measured by image analysis. Significance determined by Welch’s t test (F–G) or ANOVA with Tukey’s test (H). (I) Doubling times of NCIB 3610 and indicated isolates during exponential growth at 30 °C. p ≤ 0.05, Welch’s t test. (J) Interaction outcomes with NCIB 3610 (mTagBFP+) as challenger and the indicated strains (gfp+). ***p < 0.0001, ANOVA with Tukey’s test. (K) Doubling times of NRS6942, NRS7279, and NRS7771 at 30 °C. p ≤ 0.01, ANOVA with Tukey’s test. Note - Data for NCIB 3610 in Fig 5A, 5J and 5H & NCIB 3610 ΔcomA in Fig 5J and 5H are the same dataset.

We next assessed whether the impact of deleting comP in NRS6105 was mediated through its cognate response regulator ComA. To test this, we generated an NRS6105 ΔcomA mutant (GFP+) and challenged it against NCIB 3610 (mTagBFP+) in dual-isolate biofilms. Like the comP mutant, the ΔcomA strain displayed increased survival compared to the parental NRS6105 isolate (Fig 5A). These findings suggest that disruption of the ComP–ComA regulatory axis enhances competitiveness of NRS6105 in mixed biofilms. This phenotype mirrors the case in NCIB 3610, where deletion of rapP elevates ComA~P levels and reduces its competitiveness [20,22].

To assess whether the competitive advantage conferred by NRS6105 comP and comA mutations was specific to interactions with NCIB 3610, we also competed these mutants against their parental NRS6105 background. In both cases, the mutants strongly outcompeted the parental strain, occupying nearly the entire colony biofilm when starting from a 1:1 inoculum (Fig 5D–5E). Similar results were obtained for evolved isolates: comP mutants of NRS6145 and NRS6202 dominated their respective parental strains (Figs 5F–5G and S5A-S5B). In contrast, comP or comA mutations introduced into NCIB 3610 itself produced only modest changes when competed against the NCIB 3610 parent consistent with the presence of RapP^3610^ (Fig 5H).

Like RapP-PhrP^3610^ [22], ComQXPA is linked to a trade-off in growth and antibiotic production [29]. We next examined whether disruption of the ComQXPA system influences growth rate. Deletion of comP or comA reduced doubling times in all backgrounds tested except NCIB 3610 (Fig 5I). For NRS6105 ΔcomP, as well as evolved NRS6145 L1P4 and NRS6202 L6P4 strains, growth rates during exponential phase approached those of NCIB 3610 and differed significantly from their parental isolates. These findings are consistent with unequal growth being a major driver of exclusion of parental strains by their comP-mutant derivatives, and with accelerated growth underpinning the ability of these mutants to coexist with NCIB 3610.

Finally, to complete our interrogation, we investigated the relationship between ComA and RapP in NCIB 3610. Since NCIB 3610 already exhibits low ComA ~ P activity due to RapP-PhrP^3610^ [20] we hypothesised that comA would act epistatically to rapP. To test this, we constructed an NCIB 3610 ΔrapP ΔcomA double mutant and compared its biofilm phenotype with the single mutants. The double mutant more closely resembled the ΔcomA strain in biofilm morphology and the surface area (footprint) occupied by the biofilm (S5C-S5D Fig). In competition assays, the double mutant displayed enhanced fitness relative to ΔrapP alone but did not reach the coexistence level observed for ΔcomA (Figs 5J and S5E). Notably, this occurred despite recovery of the doubling time shown by the ΔcomA mutant to wild-type levels (Fig 5K). These results suggest that RapP’s impact on NCIB 3610 competitiveness is only partially mediated through ComA, and that equalising growth rates alone does not guarantee coexistence in colony biofilms.

Towards identification of specific competition determinants

To explore the molecular basis of NCIB 3610’s competitive advantage, we performed RNA sequencing to obtain the global transcriptional profiles of NCIB 3610 (gfp+) and its ΔrapP ΔcomA (gfp+) derivative. These strains were chosen because they show equivalent growth rates in planktonic culture, yet NCIB 3610 consistently outcompetes the ΔrapP ΔcomA mutant in dual-isolate colony biofilms (Fig 5). RNA was isolated from each strain after 20 hours of incubation under colony biofilm conditions and transcriptomic analysis revealed 469 differentially expressed transcripts between B. subtilis NCIB 3610 and its ΔrapP ΔcomA mutant, as defined by FDR < 0.01 and |logFC| > 1 (S6A Fig).

Despite the comparable growth rates in planktonic culture, KEGG pathway and GO-term analysis revealed down-regulation of growth and metabolism associated pathways in ΔrapP ΔcomA (S3 Table). This is consistent with the rewiring of metabolism that occurs when rapP^3610^ is deleted [22]. While metabolic rewiring may contribute to the interaction outcome, proteinaceous toxins and differential antibiotic expression have also been shown to influence intraspecies interactions. For example, toxins secreted by the Type VII secretion system (T7SS) influence competition among B. subtilis isolates grown in biofilms [13]. Therefore, we examined expression of 89 known toxin and toxin-associated transcripts (S4 Table).

We noted that transcripts from the yfj, pks, ywq, alb, and wap toxin-encoding operons were significantly down-regulated in the ΔrapP ΔcomA strain (S4 Table). Of these, the yfj operon has a minimal impact on intraspecies competition in biofilms (Kobayashi, 2023) and therefore was excluded from further analysis. The pks operon, controlling bacillaene synthesis, is linked to interspecies competition (Vargas-Bautista et al., 2014). The operon is present in all isolates in our analysis [14] and therefore its simple presence or absence does not correlate with the intraspecies outcomes we observe and was excluded. For the same reason [14], the alb operon, which contributes to the export of SboA, the bacteriocin subtilosin (Zheng et al., 1999) was excluded. The wap operon is polymorphic and contains the wapA toxin and wapI immunity genes. WapA mediates contact-dependent competition (Koskiniemi et al., 2013) but while the ΔwapAI mutant was outcompeted by the parental strain in liquid culture there was a lesser impact in biofilms (Kobayashi, 2023). This left the ywq operon (Fig 6A and S4 File) which is linked with intraspecies competition mediated by the Type VII secretion system [13]. The RNA sequencing analysis showed that expression of all members of the ywq operon were downregulated over 2-fold (Fig 6B). In additional support of YwqJ having a role as an intraspecies competition determinant, our bioinformatic analysis of the T7SS toxin genes (which contain LXG-domains) carried by NCIB 3610 revealed considerable variation in their prevalence across the environmental isolates used in this study (Figs 6C and S6B). More specifically, we noted that isolates that were better able to coexist with NCIB 3610 derived strains displayed greater similarity in their T7SS repertoires compared to isolates that were excluded. Therefore, to investigate the impact of ywqJ in mediating ΔrapP ΔcomA competitive fitness, we constructed a deletion of ywqJ in NCIB 3610 (mTagbfp+) (NRS8185). A dual isolate biofilm competition assay between ΔrapP ΔcomA (gfp+) (NRS7771) and the ywqJ mutant displayed higher survival rates for ΔrapP ΔcomA (gfp+) compared to a competition assay between ΔrapP ΔcomA (gfp+) and the parental NCIB 3610 (gfp+) strain (Figs 6D and S6C-S6D). Collectively, these results indicate that ywqJ contributes to intraspecies competition, but it is likely not to be the sole determinant of the interaction outcome as the two competing strains did not remain in equal prevalence after co-culture suggesting the involvement of other factors.

*The Type VII section system toxin YwqJ is a key determinant.(A) Schematic of the ywq operon based on genomic information extracted from Subtiwiki [30]. (B) RNA sequencing log2FC values for the members of the ywq operon. (C) Heatmap depicting the % protein sequence identity of the top matched LXG containing protein from the suite of B. subtilis isolates (x-axis) relative to the indicated LXG protein encoded by NCIB 3610 (y-axis). The competition outcome of the different wild isolates when cocultured with NCIB 3610 ΔrapP which has an equivalent doubling time is shown (data from Fig 2A). (D) Interaction outcomes with ΔrapP ΔcomA (gfp+) (NRS7771) as challenger and the indicated strains (mTagBFP+). ***p < 0.0001, significance determined by Welch’s t test.

Discussion

Quorum-sensing disruption as a strategy for competitive advantage

Intraspecies competition is a fundamental force shaping microbial community structure and evolution. Here we show that in B. subtilis, disruption of key quorum sensing pathways, either through natural variation or experimental evolution, confers a significant fitness advantage, even in structured biofilm communities where evolutionary pressure might be expected to favour communal traits. Specifically, we demonstrate that the RapP-PhrP module in strain NCIB 3610 and spontaneous mutations in the comP sensor kinase promote increased growth and competitive exclusion of otherwise closely related strains. However, growth is not the only determinant of the competition outcome as isolates with comparable growth rates still show antagonism. Antagonistic interactions were influenced by the profile of Type VII secretion system substrates.

In NCIB 3610, the RapP-PhrP^3610^ system is rendered unresponsive to its cognate signal due to a single amino acid substitution in RapP, leading to sustained inhibition of ComA [20]. The N236T mutation in RapP effectively disables signal responsiveness, decoupling quorum sensing from cell density [20]. This rewiring suppresses ComQXPA-dependent signalling, accelerates growth [20,22] and enhances competitive ability in biofilm environments. It also renders NCIB 3610’s behaviour atypical relative to true wild strains, particularly in social traits like secondary metabolite production, swarming, and biofilm formation [10,20,21]

Experimental evolution and parallels across bacterial systems

Experimental evolution of strains typically outcompeted by NCIB 3610 repeatedly yielded mutations in comP that mimic the effects of RapP^3610^; that is, reduced ComA activity, increased growth, and improved fitness during co-culture with NCIB 3610. These findings underscore a central theme: the suppression of communal signalling pathways can favour individual competitiveness. This trade-off between quorum-sensing–regulated cooperation and rapid growth is increasingly recognised across bacterial systems. In Pseudomonas aeruginosa, for example, lasR mutants that disrupt quorum sensing frequently emerge during chronic infections, displaying enhanced growth while losing virulence factor production and public goods sharing [31,32]. Additionally, in Streptomyces the inactivation of costly antibiotic biosynthesis pathways improves growth under nutrient-limiting conditions [33]. Spontaneous quorum sensing-defective mutants have been observed in diverse bacteria during laboratory culture, including Burkholderia glumae [34] and P. aeruginosa [35]. Consistent with these examples, our findings show that disruption of the comQXPA or the presence of the RapP-PhrP variant that constitutively antagonises the ComQXPA system can offer a competitive edge in structured microbial communities. Importantly, such quorum sensing disruptions are not restricted to experimental evolution. Our analysis of 370 B. subtilis genomes shows that disruptive comP mutations are not rare but present in at least 16% of wild isolates, particularly in those sharing the NCIB 3610 ComQXPA pherotype. It is possible, however, that in the NCIB 3610 ComQXPA pherotype, the mutation which frequently occurs within the string of eight homopolymeric bases is akin to phase variation and therefore may be reversible. Overall, the reoccurring detection of multiple different mutations in comP suggests a recurrent and naturally selected strategy whereby quorum sensing is attenuated to reduce the cost of social engagement under competitive pressure.

Our analysis further explores the complex regulatory interplay between the Rap–Phr and ComQXPA systems (Fig 1A). The competitive advantage conferred by RapP^3610^ is only partially dependent on ComA, as evidenced by the intermediate phenotype of the rapP comA double mutant. RapP also influences Spo0A activity by dephosphorylating Spo0F, thereby altering matrix production and sporulation [4]. Moreover, during the evolution experiments we isolated a derivative of NRS6145 with a deletion of rhgR1 (Fig 3), which encodes a repressor of several Rap phosphatases including RapD, RapG, and RapH [36]. Loss of rhgR1 would be expected to derepress these Rap proteins, enabling them to dephosphorylate their targets and perhaps further modulate ComA and Spo0A activity. Thus, in the NRS6145-L6P4 rhgR1 containing mutant, multiple Rap-mediated inputs may be simultaneously deregulated, expanding the regulatory pathways that shape competitive outcomes. As Rap-Phr systems integrate environmental cues such as nutrient availability and extracellular peptide signals [20,37], perturbations at this level are likely to propagate broadly.

Beyond growth: the role of antagonistic effectors

While growth rate is a strong predictor of competitive fitness, our findings make clear that it is not the sole determinant. Strains with equivalent doubling times exhibited contrasting outcomes depending on their regulatory configuration. For instance, NCIB 3610 and its rapP comA double mutant grew at similar rates in planktonic culture, but only the parental strain dominated in biofilms. Using RNA sequencing analysis, we were able to narrow one of the contributors of the competition outcome to differential expression of the Type VII secretion substrate YwqJ. Toxin-mediated interference, such as through LXG-domain proteins, can confer dominance even at slower growth rates [10,13], and differential expression on deletion of rapP and comA is consistent with the knowledge that the production of antibiotics or matrix components, often ComA-regulated [38,39], can tip the balance depending on ecological context.

Context-dependent fitness of quorum-sensing mutants

Crucially, the benefit of disabling ComA-based quorum sensing appears to be context-dependent. In our dual-B. subtilis strain model, where all competitors belong to the same species and are relatively resistant to each other’s secondary metabolites, the energetic cost of producing ComA-regulated traits likely outweighs their utility. In contrast, it has recently been shown that B. subtilis invests more heavily in secondary metabolites when confronted with phylogenetically distant species [29], precisely when such compounds are more likely to be effective. Conversely, B. subtilis evolved under fungal competition (Setophoma terrestris) acquired mutations in comQXPA that led to increased production of antifungal volatiles, enhancing inter-kingdom fitness [40]. These findings underscore that while quorum sensing repression is beneficial in intraspecies competition, it may be detrimental in broader ecological settings where toxin production is essential for survival.

Taken together, our study reveals that modulation of quorum sensing via genetic disruption or natural variation represents a powerful mechanism to mediate competition in microbial populations. It highlights the evolutionary plasticity of regulatory networks and suggests that the ecological value of cooperation is highly contingent on community composition and environmental structure. A systems-level understanding that integrates growth dynamics, regulatory complexity, spatial ecology, and inter-strain antagonism will be essential to predict and manipulate microbial community outcomes.

Materials and methods

Growth Media and additives

For routine growth of Bacillus subtilis, lysogeny broth (LB) liquid media was made using the following recipe: 1% (w/v) Bacto-peptone, 1% (w/v) NaCl, 0.5% (w/v) yeast extract. For solid plates, LB broth was supplemented with 1.5% (w/v) agar. When necessary, LB media cultures and plates were supplemented with antibiotics which were used at the following concentrations for B. subtilis: 10 μg/ml kanamycin, 100 μg/ml spectinomycin, 5 μg/ml chloramphenicol and 0.5 μg/ml erythromycin combined with 12.5 μg/ml lincomycin for “MLS”. Ampicillin was used at 100 μg/ml for work with E. coli.

Plasmid and strain construction

The strain used for storing of plasmids for cloning was Escherichia coli strain MC1061 [F’ lacIQ lacZM15 Tn10. For making mutations in the NCIB 3610 background, as this strain is not genetically competent, plasmids were first transformed into the laboratory strain 168 using standard protocols (S1 Methods).

Markerless in-frame deletion strain construction

For construction of strains with in-frame deletions without the insertion of antibiotic resistance cassettes into their genomes, pMiniMAD-based plasmids [41] were used as detailed for the parental pMAD plasmid [42] (S5-S7 Tables and S1 Methods).

Setting up single strain colony biofilms

B. subtilis isolates were streaked out on LB agar plates and incubated at 37°C overnight. The following day single colonies were grown in 3 ml of LB broth at 37°C with agitation. The cultures were grown to an OD_600_ of 1 and 10 µl of the cultures were spotted onto MSgg media plates. The plates were incubated at 30°C for 48 h before imaging. Biofilm imaging was performed using a Leica MZ16 FA stereoscope, LAS version 2.7.1. The imaging data were imported to OMERO [43].

Biofilm co-culture assays and surface area calculation

For dual isolate biofilm assays, cultures of the individual strains to be used were set up in 5 ml of LB and incubated at 37^o^C with agitation overnight. The following morning, day cultures were set up by inoculating 3 ml of LB with 200 µl of the overnight cultures. The day cultures were incubated at 37^o^C with agitation. The growth of the cultures was monitored, and all cultures were normalised to an OD_600_ of 1. After normalisation, cultures were mixed at a 1:1 ratio as required. 5 µl drops of the culture mixtures were spotted onto MSgg agar plates and 5 µl of the individual normalised cultures were included as controls. The plates were incubated at 30^o^C and imaged after 24, 48 and 72 hours as required. Fluorescence imaging was performed using a Leica fluorescence stereoscope (M205FCA) with a 0.5 × 0.2 NA objective. Imaging files were imported to OMERO [43]. To quantify the biofilm footprint, we used OMERO.Insight v.5.8.1 (Allan et al., 2012). Each biofilm was highlighted as a region of interest, and the area occupied was calculated. Relative strain densities of GFP and mTagBFP-expressing cells in mixed biofilm assays were determined by analysing fluorescence image data [24]. ImageJ/Fiji [44] was used to run the macro as previously described [24].

Confocal Analysis of Colony Biofilms

Confocal imaging was done using methodology similar to that described previously [45]. See S1 Methods.

Recovery of isolates from dual isolate colony biofilms

A colony biofilm was removed from the agar plate with a sterile plastic loop and placed in 500 µL of phosphate-buffered saline (PBS): 137 mM NaCl, 2.7 mM KCl, 8 mM Na_2_HPO_4_, and 2 mM KH_2_PO_4_. The biofilm was disrupted by passage through a 1 mL syringe with a 23G needle 6 times, a process repeated with a 27G needle. Before plating, the disrupted biofilm was serially diluted in PBS and 100 µL was spread onto LB agar plates containing the antibiotic to which the desired strain was resistant. Colonies recovered were counted after overnight incubation at 37^o^C and enumerated where required.

Isolation of cells from dual isolate colony biofilms for flow cytometry

Mixed isolate colony biofilms were set up by growing the cultures to an OD_600_ of 1 and mixing the strains of interest at a 1:1 ratio as described above. After 48h, a colony biofilm was removed from the agar plate with a sterile plastic loop and placed in 500 µL of 4% formaldehyde before processing (S1 Methods).

Experimental evolution

Strains NRS6936, NRS6951, and NRS7201 (NRS6105, NRS6145 and NRS6202 respectively adapted to express mTagBFP and containing spectinomycin resistant) and NRS1473 (NCIB 3610 expressing GFP, kanamycin resistant), were used for the experimental evolution detailed in the S1 Methods.

Genomic DNA extraction and sequencing

Genomic DNA was extracted using the QIAprep Spin Miniprep Kit (cat. Nos. 27104 and 27106) and following the “Quick-Start Protocol”. DNA concentration was measured using a Qubit 4 Fluorometer or a NanoDro One/One^C^ Microvolume UV-Vis Spectrophotometer. Whole Genome Sequencing was performed by GENEWIZ Germany GmbH (Part of Azenta Life Sciences), MicrobesNG (UK) or SeqCentre (USA) as indicated in the strain table (S1 Table). See S1 Methods for more information.

Bacillus subtilis phylogenetic tree

Sequence selection and preparation for phylogenetic analysis was conducted using custom Python code co-ordinated through a Snakemake workflow [46]. A database of ‘complete’ genome B. subtilis sequences was constructed from the European Nucleotide Archive (with the genome represented in <5 contigs, to allow for plasmid sequences) which had been assigned an NCBI taxonomy ID of 1423, including child taxa, resulting in a set of 478 genomes (Data downloaded on 6^th^ June 2025) (S2 Table). Genome sequences were reannotated with Bakta 1.11.0 [47] (database version 5.1) to ensure consistent annotation. Genome completeness was assessed using Busco 5.8.3 [48] using the ‘bacillales_odb10’ lineage, and the database filtered to remove genomes with a Busco completeness of <98%. Taxonomic classification was carried out using GTDB-TK 2.4.1 (GTDB reference data version r226) [49] and misclassified sequences removed from the database, retaining only those classified as B. subtilis. The final database contained 370 genome sequences.

Analysis of comP sequences in Bacillus subtilis genomes sequences

Operon sequences from the start of c**omQ to the end of c**omA were extracted from the annotated genome sequences and analysed (S1 Methods and S8 Table).

RapP and LXG protein sequence analysis

Proteins in the pfam family PF04740 LXG domain of WXG superfamily [50] were identified in each strain using HMMER [51].

Measurement of Growth Rate

To calculate the doubling time of cultures, strains were streaked out and incubated for 16 hours at 30^o^C. Liquid overnight cultures were established in 5 ml of LB and incubated at 37^o^C with agitation at 250 rpm for 16 hours. Day cultures were established by dilution (200 μl into 5 ml) and grown at 37^o^C with shaking for 3 hours and 10 mins. The tubes were centrifuged, and the cell pellets were resuspended in 1 ml of MSgg. The OD_600_ was measured, and the culture normalised to an OD_600_ of 1 using MSgg. Cultures for growth measurement were set up to an initial OD_600_ of 0.01 using 25 ml of pre-warmed MSgg and incubated at 30^o^C with 250 rpm agitation. OD_600_ values were recorded every 30–60 minutes the doubling time calculated using the exponential growth phase measurements.

Structural modelling of ComP

The protein sequences of ComP from strain NRS6145 (wild type, WT) and the mutant variant carrying the E510K substitution were submitted to the AlphaFold2 server [26] (S8 Table and S1 Methods). The same process was used for NRS6105.

RNA isolation and RNAseq analysis

RNAseq analysis was conducted on strains NRS6942 and NRS7771, grown for 20 h at 30^o^C on MSgg agar plates. Five biological replicates were used, using 4 technical replicates pooled from different plates. The RNA isolation protocol was adapted [52] using a Qiagen RNeasy RNA Isolation kit (S1 Methods). Library construction and RNA sequencing was performed by Edinburgh Genetics using Stranded Total RNA Prep with Ribo-Zero Plus Microbiome kit. RNA was sequenced using Illumina Paired Ends 100 bp, with>=12M PE reads per sample. The data processing is detailed in the S1 Methods.

Mathematical model of dual-isolate biofilms

We used a mathematical model to test hypotheses on the impact of an isolate’s doubling time on its competitiveness during biofilm growth. Adapted from an existing mathematical model [24], the model describes the spatio-temporal dynamics of two isolates in a single growing biofilm (S1 Methods).

Data analysis

Graphs were constructed using GraphPad Prism 9 or R/ggplot2. Statistical analysis methods are detailed in the legends.

Supporting information

S1 FigRemoval of rapP reverts NCIB 3610 colony biofilm architecture to that of soil isolates.(A) Colony biofilm morphology of NCIB 3610 ΔrapP and other soil isolates of B. subtilis. Strains were grown on biofilm-inducing media for 48h at 30^o^C before imaging. The strain and genotype are indicated. The scale bars represent 0.5 cm. All soil isolates shown here are variants constitutively expressing GFP and the image of NCIB 3610 rapP is of an mTagBFP-expressing variant. (B) Representative outcome after 48 hours incubation in a dual isolate colony biofilm incubated at 30^o^C for 48 hours. (C) Quantification of the interaction outcome represented as %challenger remaining. The cocultured strain is indicated on the x-axis. Each data point represents a single colony biofilm and is derived from a combination of technical and biological repeats. The error bars represent the standard deviation of the mean.(TIF)

S2 FigOccurrence of rapP and interaction outcome after pairwise screen.(A) Quantification of pairwise biofilm competitions shown as % challenger remaining. Challenger was NRS6202 (mTagBFP+). Partner strains are indicated on the y-axis. Each point = one biofilm; error bars = SD. (B) Prevalence of the rapP coding region in the genomes of 370 Bacillus subtilis isolates. B. subtilis isolates indicating isolates carrying rapP, and the state of the N236T residue associated with reduced signal responsiveness (C) Representative images of single isolate and dual isolate biofilms linked to Fig 2G. The scale bars represent 10 mm. The strains are as indicated. Images were taken after 48 hours incubation at 30^o^C.(TIF)

S3 FigLineages of Bacillus subtilis subjected to laboratory evolution.(A-C) Colony biofilms of the evolving soil isolates clones mixed with NCIB 3610 grown on MSgg agar plates for 24h at 30^o^C before imaging. “mixed” shows the NCIB 3610 (green) and soil isolate (magenta) co-culture biofilms. Panel A shows lineages for NRS6105; panel B shows lineages for NRS6145; panel C shows lineages for NRS6202. P1-P6 on the left-hand side indicates the passage in the experimental evolution experiment and the soil isolate background and lineage shown in each panel are indicated at the bottom, where “L” indicates the lineage and “P” the passage number. The scale bars represent 5 mm (A and C) and 2 mm (B). (D) The area of the colony biofilms of the evolving soil isolates was quantified with respect to passage number (x-axis) (the images in the “single” panel Fig 3). n = 1. (E-F) Structural models of ComP wild type (WT) and E510R mutant. (E) AlphaFold-predicted models of the ComP dimer in the WT (pink) and E510R (purple) forms were superimposed and displayed as cartoon representations. The transmembrane (TM) helices are shown at the top and boxed to indicate the approximate position of the lipid bilayer. The periplasmic sensor domain lies between the TM helices, while the cytoplasmic portion consists of the HAMP, DHp, and the catalytic ATP-binding domain typical of histidine kinases. The yellow circle marks the location of residue 510 within at the membrane interface. (F) Close-up view of the membrane interface domain circled in (E). Side chains of key residues are shown as sticks and coloured as in panel A. Residues R331, E394, and position 510 are labelled. In the WT model (pink), E510 forms potential interactions with R331 and possibly E394 on the adjacent helix. In the E510R mutant (purple), the substitution disrupts this local hydrogen bond network, resulting in a rearranged side-chain interaction pattern (dashed lines). (G-H) Structural prediction of ComP dimer from NRS6105 (Δ633–649) vs NRS6105 WT. (G) Predicted structural models were superimposed on the Calpha of one of the two ATP binding domains (circle). The WT protein structure is shown in pink and the mutant Δ633–649 in purple. ADP is shown in sticks, and the deleted 633–649 domain is coloured cyan in the WT protein. The second ATP-binding domain does not superpose between the WT and Δ633–649 structure. (H) Close-up of the ATP-binding domains of the WT and Δ633–649 mutant. ADP and the side chain of His691 are shown in sticks. The overall ATP-binding site remains largely unchanged between the two proteins.(TIF)

S4 FigDistribution of disrupted comP sequences across phylogeny of 370 B. subtilis comQ sequences.The comQ gene was selected as a conserved member of the comQXPA operon. The phylogeny is rooted on an outgroup consisting of a Bacillus amyloliquefaciens sequence. Bootstrap values are represented by differently sized circles on the tree branches.(TIF)

S5 FigComQXPA influences the outcome of interactions in otherwise isogenic contexts.(A-B) Representative examples of colony biofilm interactions. The scale bars represent 10 mm or 3 mm. The strains in the coculture are indicated. Images were taken after 48 hours incubation at 30^o^C. (C) Colony biofilm morphology after growth on biofilm-inducing media for 48h at 30^o^C. The strain and genotype are indicated. All strains are GFP+ derivatives of NCIB 3610 (NRS6942, NRS7770, NRS7279 & NRS7771). The scale bars represent 1.0 cm. (D) Footprint measurements for colony biofilms of the genotypes detailed in the legend of the graph. Ordinary One-way ANOVA Tukey’s multiple comparisons test, with a single pooled variance, was used for analysis. The asterisks represent statistical significance with a p-value of ≤0.05 (); ≤ 0.01 (); ≤ 0.001 () or ≤0.0001 (****). (E) Colony biofilm morphology after growth on biofilm-inducing media for 48h at 30^o^C. The scale bars represent 10 mm. The strains in the coculture are indicated.(TIF)

S6 FigIdentification of a specific competition determinant.(A) Volcano plot of the 469 differentially expressed genes identified by RNA sequencing. The toxins and other antimicrobials that are expressed to a higher level in NCIB 3610 compared to the ΔrapP ΔcomA strain are indicated. (B) Phylogenetic tree based upon Muscle alignment of the top matched YwqJ LXG containing protein sequence from the suite of B. subtilis isolates relative YwqJ encoded by NCIB 3610. The LXG domain containing protein YeeF from NCIB 3610 was used as an outgroup. The heatmap represents the competition outcome of the different wild isolates when cocultured with NCIB 3610 ΔrapP which has an equivalent doubling time is shown (data from Fig 5K). (C) and (D) A representative outcome of a dual isolate colony biofilm morphology after growth on biofilm-inducing media for 48h at 30^o^C. The strains and genotypes are indicated. The upper panel is bright field, and the lower is the FITC signal false coloured green. The scale bars represent 1.0 cm.(TIF)

S1 TableStrains used in this study.(XLSX)

S2 TableGenomic sequences used in this study.(XLSX)

S3 TableFunctional terms/pathways derived from GO-terms and KEGG pathways.(XLSX)

S4 TableList of 89 toxins and antibiotics.Gene symbols provided were matched to the SubtiWiki gene symbols, mapped to NCBI annotations.(XLSX)

S5 TablePrimers used in this study.(XLSX)

S6 TablePlasmids used in this study.(XLSX)

S7 TableSequences of constructs synthesized.(XLSX)

S8 TableProtein sequences used in this study.(XLSX)

S1 FileJalview RapP alignments.File that can be opened using open-source software containing all the sequence information linked to the analysis of RapP. https://www.jalview.org/(JVP)

S2 FileBreseq Analysis.Results from the analysis of short read sequence data for a given strain against a complete genome sequence to identify genomic different. Each tab in the file indicates a different analysis where the test strain identification number is given first before the reference genome strain number used for analysis.(XLSX)

S3 FileJalview com alignments.File that can be opened using open-source software containing all the sequence information linked to the analysis of comQXPA. https://www.jalview.org/(JVP)

S4 FileJalview YwqJ alignments.File that can be opened using open-source software containing all the sequence information linked to the analysis of YwqJ. https://www.jalview.org/(JVP)

S1 MethodsRegulatory Rewiring Drives Intraspecies Competition in Bacillus subtilis.(DOCX)