Interaction variability shapes succession of synthetic microbial ecosystems

Liu, Feng; Mao, Junwen; Kong, Wentao; Hua, Qiang; Feng, Youjun; Bashir, Rashid; Lu, Ting

doi:10.1038/s41467-019-13986-6

Download PDF

Article
Open access
Published: 16 January 2020

Interaction variability shapes succession of synthetic microbial ecosystems

Feng Liu^1,2,3,4^na1,
Junwen Mao^1,2,5^na1,
Wentao Kong^1,2,
Qiang Hua^3,4,
Youjun Feng⁶,
Rashid Bashir^1,7,8 &
…
Ting Lu^1,2,9,10,11

Nature Communications volume 11, Article number: 309 (2020) Cite this article

7848 Accesses
29 Citations
47 Altmetric
Metrics details

Subjects

Abstract

Cellular interactions are a major driver for the assembly and functioning of microbial communities. Their strengths are shown to be highly variable in nature; however, it is unclear how such variations regulate community behaviors. Here we construct synthetic Lactococcus lactis consortia and mathematical models to elucidate the role of interaction variability in ecosystem succession and to further determine if casting variability into modeling empowers bottom-up predictions. For a consortium of bacteriocin-mediated cooperation and competition, we find increasing the variations of cooperation, from either altered labor partition or random sampling, drives the community into distinct structures. When the cooperation and competition are additionally modulated by pH, ecosystem succession becomes jointly controlled by the variations of both interactions and yields more diversified dynamics. Mathematical models incorporating variability successfully capture all of these experimental observations. Our study demonstrates interaction variability as a key regulator of community dynamics, providing insights into bottom-up predictions of microbial ecosystems.

Short-range interactions govern the dynamics and functions of microbial communities

Article 10 February 2020

Non-additive microbial community responses to environmental complexity

Article Open access 22 April 2021

Rare and localized events stabilize microbial community composition and patterns of spatial self-organization in a fluctuating environment

Article Open access 25 January 2022

Introduction

Microbial communities are assemblages of multi-species microorganisms that live and interact with each other. They regulate the biogeochemical cycling of the planet¹, fertilization of crops², and metabolism of our human body³, thereby affecting profoundly the environment, agriculture and human health. As the dynamics of a community often underlies its emergent properties, being able to predict ecosystem succession is central to the elucidation of community organization and function⁴, rational design of artificial ecosystems^5,6,7 and introduction of intervention⁸. To that end, rapidly developed is bottom-up analysis^9,10,11,12, which aims to determine community behaviors from the characterization of microbial interactions—such as competition and cooperation—that are ubiquitous among microbes¹³. This approach is conceptually compelling, because it offers a systematic solution to capture emergent ecosystem properties and is potentially generalizable for both native and synthetic communities. In practice, although there are exciting successes in selected cases such as the utilization of the generalized Lotka-Volterra model for specific microbiomes^14,15,16, predicting community dynamics from the bottom up remains as a grand challenge in general⁴.

One possible major cause of this challenge is pointed to the discrepancy between current modeling scheme and experimental observations. To date, bottom-up models commonly consider microbial interactions invariant and presume that determining interaction strength autonomously specifies community behaviors. Associated with these model developments are experimental efforts, which largely focus on the identification and measurement of interaction strength^4,13. However, increasingly overwhelming evidences show that microbial interactions are highly variable, rather than static, in nature. Indeed, microbial interactions often change with environmental cues such as pH, nutrient and stress^{17,18,19,20,21,22,23,24,25,26,27}. For example, the mutualism between Escherichia coli and Rhodopseudomonas palustris is moderated by the toxicity of organic acids in the culture²¹; the antagonism from Pseudomonas aeruginosa to Staphylococcus aureus increases with ion depletion²⁴. Microbial interactions are also subjective to the population of microorganisms generating the interactions as well as the presence of other species^{20,28,29,30,31,32}. For instance, Lactococcus lactis produces nisin to suppress pathogens such as Staphylococcus aureus through quorum sensing of its own population³¹ while Enterococcus faecalis secrets cytolysin when sensing the presence of target cells³². Additionally, as cellular interactions are typically fulfilled through the production of metabolites and proteins—biochemical processes that are fundamentally stochastic^33,34,35, there are intrinsic fluctuations for all microbial interactions including those ‘constant’. Recognizing this characteristic of microbial interactions, a handful of mathematical frameworks have been proposed to consider interaction variations^{36,37,38,39,40}. In contrast, there is a lack of systematic experimental investigations that quantify the degree of variability for given microbial interactions. It also remains unclear to what extent such variations drive ecosystem succession and alter community structures and characteristics. Accordingly, it is unknown how the incorporation of variability into modeling shapes the predictive power of bottom-up mathematical modeling.

Here we hypothesize that variations of microbial interactions are a key modulator of community behaviors and characterizing and incorporating the variability empowers predictive understanding of ecosystem succession from the bottom up. To test the hypotheses, we design and build a set of synthetic three-strain microbial consortia, which involve both cooperation and competition, and use them as our experimental model systems. Compared to native ecologies, such synthetic communities possess a significantly reduced degree of complexity while offering the feasibility for mechanistic dissection and quantitative measurement^{7,41,42,43,44,45,46,47,48,49}. In parallel, we develop mathematical models with an explicit incorporation of interaction variability to analyze ecosystem succession. For the consortium containing a variable cooperation, we quantify the variability of cooperation, elucidate the alteration of ecosystem dynamics arising from the variations, and demonstrate the power of variability-incorporated modeling in capturing community development. For the ecosystem whose cooperation and competition both fluctuate, more complex ecosystem dynamics arises but characterizing the variabilities again lead to successful succession predictions. Together, our results elucidate the role of interaction variability in regulating community dynamics, providing fundamental insights into bottom-up understanding of microbial ecosystem succession.

Results

Creation of a cooperation between synthetic populations

We started by engineering a cooperation in synthetic populations because it is ubiquitous among microorganisms and critical to their organization^50,51,52. Specifically, the interaction involves two engineered Lactococcus lactis strains, Cα and Cβ, both of which harness the biosynthetic pathway of lcnG, a Class II two-subunit lactococcus bacteriocin (Fig. 1a and Methods)^53,54. Here, Cα constitutively expresses lagA, lagD, and lagE which encode the peptide α precursor, ABC transporter and accessory protein of the pathway respectively, allowing the strain to synthesize and secret the α subunit of lcnG. Similarly, Cβ constitutively expresses the pathway’s β precursor gene lagB, transporter gene lagD and accessory protein gene lagE, enabling the synthesis and secretion of the β peptide, the other subunit of lcnG. In the extracellular milieu, the two subunits α and β self-assemble into a bioactive antimicrobial which inhibits the growth of L. lactis strains. Through this fashion of division of labor, Cα and Cβ achieve a cooperation for successful lcnG production. As lcnG inhibits all L. lactis strains, the immunity gene lagC was introduced into both Cα and Cβ to confer them an immunity. Additionally, two reporter genes, yemGFP and mCherry, were loaded into Cα and Cβ respectively to enable the quantification of ecosystem dynamics.

**Fig. 1: An engineered cooperation between two bacterial strains.**

To validate the cooperation, we conducted inhibition zone assays using the supernatants of Cα and Cβ monocultures (Methods). Four supernatant combinations, including blank culture (GM17 media supplemented with chloramphenicol) (−/−), supernatant of Cα monoculture ( + /−), supernatant of Cβ monoculture (−/ + ), and mix of the two supernatants ( + / + ), were loaded into single wells in the solid agar plated with a lawn of lcnG-sensitive cells (L. lactis NZ9000 loaded with a chloramphenicol resistant gene)⁵⁵. Upon 8 hours (h) of incubation at 30 °C (Methods), only the well loaded with the supernatant mix ( + / + ) produced an inhibition zone (Fig. 1b). In addition, the four supernatant combinations were loaded into separate wells in agar plates covered with lcnG sensitive cells. We found that a clear inhibition zone formed between adjacent wells only when they were loaded with the Cα and Cβ supernatants respectively (Fig. 1c). Together, the results confirmed two pieces of information. First, combination of Cα and Cβ produced the active lcnG but individual strains alone did not. Second, peptides α and β were able to autonomously assemble into an active bacteriocin without the need for any assistance.

Characterizing the variations of cooperation

Driven by diverse biotic and abiotic factors, microbial interactions are highly variable. Here we aimed to experimentally determine the variability of the cooperation originating from its division of labor nature. Specifically, to quantitate how labor partition, reflected by the cooperator ratio, affects the strength of cooperation, we grew the monocultures of Cα₀ and Cβ₀ (the reporter-free version of Cα and Cβ), mixed their supernatants with varied ratios while fixing the total volume (30 μL), and further used the mixes to perform inhibition zone experiments (Methods). Our results (Fig. 2a, top row) showed that, across the ratios from 30:1 to 1:30, the size of the inhibition zones varies from small to large and back to small with the largest occurring at the 1:1 ratio. As zone size correlates positively to a mix’s lcnG level and hence the strength of cooperation, the results suggested that initial labor partition can induce significant variations of the cooperation. Using the supernatants of Cα₀-Cβ₀ co-cultures with varied initial relative population abundance, we also observed the same dependence of the lcnG level on the initial population partition (Fig. 2a, bottom row). To quantitatively determine the variations, we further measured the relative lcnG level of each sample by normalizing its inhibition zones with that of Cαβ, a lcnG-producing strain (L. lactis MG1363 loaded with the complete lcnG pathway) (Fig. 2b, Methods, Supplementary Fig. 1).

**Fig. 2: Characterization of cooperation variations.**

In the above experiments, a low-high-low pattern of the cooperation strength (i.e., lcnG productivity) was observed with the maximum at the 1:1 labor partition. From a molecular perspective, we reasoned that the bell-shape variation of the cooperation is rooted in the 1:1 stoichiometric ratio of the α and β subunits. We further speculated that the 1:1 initial ratio happened to be optimal because Cα₀ and Cβ₀ have a comparable growth rate and a comparable subunit productivity and, thus, the final ratio of the two peptides in a co-culture is solely determined by their initial densities. Supporting the speculation, the inhibition zone assays and lcnG quantifications for the Cα-Cβ co-culture (Supplementary Fig. 2b, c) showed that the maximum was shifted towards a higher initial Cβ abundance due to a slower growth rate of Cβ than Cα (Supplementary Fig. 3a). Such a growth reduction is a common consequence of heterologous protein production in microorganisms that is well documented in literature^56,57. Meanwhile, the maximum for the monoculture mixes remained at 1:1 ratio (Supplementary Fig. 2a, c) because, in monocultures, Cα and Cβ had the same amount of nutrient and, thereby, produced a comparable level of subunits. To further confirm that the maximal cooperation is characterized by the stoichiometric ratio of the lcnG subunits, we derived new strains with altered peptide productivities, including Cα_0v2 with a doubled α productivity compared to Cα, and Cβ_0v1/3 and Cβ_0v1/5 whose β productivities are reduced to 1/3 and 1/5 of that of Cβ₀ (Supplementary Fig. 4). In theory, the combinations of Cα₀ and Cβ_0v1/3, Cα₀ and Cβ_0v1/5, Cα_0v2 and Cβ_0v1/3, and Cα_0v2 and Cβ_0v1/5 would shift the optimal ratios from 1:1 to 1:3, 1:5, 1:6, and 1:10 respectively, which were subsequently confirmed by the experiments (Fig. 2c).

For any microbial ecosystems, there are intrinsic random fluctuations of cellular populations arising from various stochastic processes^{33,34,35,58,59}, which motivated us to quantify the variability of cooperation from intrinsic stochasticity by specifically examining the effects of sampling of initial populations (i.e., genetic drift)^60,61. We first inoculated the consortium from an initial culture (1:1 Cα-to-Cβ ratio, 1.0 optical density at 600 nm (OD)) into fresh media through serial 1:10 dilution to generate samples with different initial ODs (10⁻², 10⁻⁴, 10⁻⁶ and 10⁻⁸) (Fig. 2d and Methods). Here, 10⁻⁸ was selected as the minimal initial OD since it is the minimal density for both Cα and Cβ to stably grow in monocultures (Supplementary Fig. 5). Next, we grew the cultures for the defined incubation time, collected their supernatants, and measured the lcnG levels (Methods). Figure 2e shows the lcnG levels of the co-cultures normalized by the productivity of Cαβ, a single-strain lcnG producer we created. We found that, with the reduction of initial OD, the mean lcnG productivity of the consortium decreased monotonically but, in the meanwhile, the sample-to-sample variation increased. Such a trend was also clearly observed in the inhibition zones formed by the co-cultures (Fig. 2f and Supplementary Fig. 6). By conducting the same experiments for the reporter-free version of the consortium, the Cα₀-Cβ₀ ecosystem (Supplementary Fig. 7), we confirmed that random sampling continued to serve as a key inducing factor of the variation.

Notably, partition alteration and random sampling are two independent sources of cooperation variations; however, they are intrinsically connected. The both alter cooperation strength by varying the division of labor among the cooperators, but the former is a controlled, deterministic alteration of the partition while the latter is unintended, stochastic alteration.

Ecosystem successions driven by cooperation variations

To elucidate the consequences of interaction variations on community behaviors, we designed a three-strain consortium composed of the two cooperators (Cα and Cβ) and a competitive third strain (Ks). Here, Ks is an engineered strain capable of constitutively secreting lactococcin A (lcnA)⁶², a bacteriocin that effectively kills all L. lactis strains including Cα and Cβ unless immunized. Experimentally, Ks was built by introducing the lcnA pathway into the L. lactis MG1363 (Fig. 3a). To efficiently count Ks in the three-strain ecosystem, it was inserted with constitutively expressed gusA3, a beta-glucuronidase gene that enables colorimetric quantification upon the supplementation of X-Gluc⁶³. As Ks does not contain the lcnG immunity gene lagC, it is sensitive to lcnG cooperatively produced by Cα and Cβ. Therefore, the community involves a cooperation between Cα and Cβ and a competition of Ks with Cα and Cβ (Fig. 3b).

**Fig. 3: Succession of a three-strain ecosystem driven by cooperation variations from labor partition.**

Meanwhile, we hypothesized that predicting ecosystem behaviors requires mathematical models that take in account the variability of interaction. To test the hypothesis, we devised a dynamic ecosystem model using a systematic, bottom-up fashion (Methods, Supplementary Information 1). Briefly, we first constructed models of Cα and Cβ monocultures (Supplementary Equations 1-2, Supplementary Fig. 3b), and used them as modules to derive a model of Cα-Cβ co-culture incorporating variability (Supplementary Equation 4) with their parameters specified with experiments (Supplementary Fig. 8). We then constructed a model of lcnA-producing Ks monoculture (Supplementary Equation 5) and determined its parameters experimentally (Supplementary Fig. 9a, b). Finally, by merging the models of the Cα-Cβ co-culture model with the Ks monoculture and characterizing the inhibitions between the modules (Supplementary Equations 6-7, Supplementary Fig. 9c-f), we obtained a model for the Cα-Cβ-Ks consortium (Supplementary Equation 9). For batch fermentations starting from a fixed (Cα + Cβ):Ks ratio (2:1), the model predicted that the ecosystem evolves temporally into Ks dominance at imbalanced cooperator partitions and Ks subordinate at close partitions (Fig. 3c), suggesting that the variations of cooperation can modulate the succession of the consortium dramatically.

To test the predictions, we experimentally assembled the consortium by mixing Cα, Cβ and Ks with altered initial ratios but a fixed total OD (10⁻²) as our model construction (Methods). Consistent with the predictions, our fermentations showed that the Ks percentage indeed increased over time and eventually became dominant when the Cα:Cβ ratio was extremely imbalanced (30:1 and 1:30); in contrast, when their abundances were close (3:1, 1:1, and 1:3), the Ks percentage declined monotonically and diminished eventually (Fig. 3d). Opposite with the Ks percentage, the relative abundances of Cα and Cβ both decreased at imbalanced Cα:Cβ ratios (30:1 and 1:30) but increased when the ratio is close (3:1, 1:1, and 1:3) (Supplementary Fig. 10a, b). Meanwhile, although the total OD of the ecosystem remained largely consistent regardless of the Cα:Cβ ratios (Supplementary Fig. 8, circles), the lcnG level changed significantly (Supplementary Fig. 8, bars): it remained limited throughout the fermentations in imbalanced scenarios but accumulated rapidly at the balanced cases, consistent with our previous characterization (Fig. 2a–c). The correspondence between higher lcnG production (Supplementary Fig. 8) and lower Ks abundance (Fig. 3c, d) suggested a strong correlation between cooperation variations and diversified ecosystem succession.

To further confirm that it is a causal relationship between interaction variations and diversified succession, we designed a Ks variant, named Kr, which is resistant to lcnG co-produced by Cα and Cβ and deficient in producing lcnA. Experimentally, Kr was established by removing the lcnA biosynthetic pathway from Ks while introducing the lcnG immunity gene lagC (Fig. 3e, Supplementary Fig. 11). The mixture of Cα, Cβ and Kr formed a control consortium where the cooperation and competition are both abolished (Fig. 3f). Our mathematical model (Methods, Supplementary Equation 11) predicted that, upon the same alterations of the initial Cα:Cβ ratio, the structure of the consortium remains largely invariant (Fig. 3g), which was subsequently verified by co-culture experiments (Fig. 3h, Supplementary Fig. 10c, d). In this case, varying initial cooperator partition continued to generate the variation of lcnG production as in the Cα-Cβ-Ks consortium (Fig. 2a–c); however, due to Kr’s resistance to lcnG, the variation was ‘insulated’ and not propagated to ecosystem succession, leading to the invariant community structure. These results confirmed that it is indeed the variations of interaction that cause the diversification of ecosystem development.

To examine if the modulation of ecosystem succession by interaction variability is specific to the origin of variation, we investigated the dynamics of the Cα-Cβ-Ks consortium upon fluctuations in sampling, another source of variation we characterized. Mathematically, we used a bottom-up strategy to create a corresponding dynamic model by introducing sampling-induced variations into the previous Cα-Cβ co-culture module (Supplementary Equation 4) and Ks monoculture module (Supplementary Equation 5) with experimentally derived parameters (Supplementary Tables 7 and 8) and combining the resulting modules into a single model (Methods). Using the model that encapsulates variations from sampling, we conducted multiple repeats of computational batch fermentations for the consortium starting with the 1:1:1 ratio but different initial ODs. In parallel, we experimentally mixed the strains Cα, Cβ and Ks in 1:1:1 ratio to form the consortium, inoculated them into ten replicates at specific initial total ODs (10⁻², 10⁻⁴, 10⁻⁶, and 10⁻⁸) and measured their population dynamics over time (Methods). Notably, due to the difference in initial conditions, the culturing time and sampling time were altered accordingly to enable a consistent and proper comparison (Supplementary Fig. 12).

Our model predictions and subsequent experiments showed that the Ks fraction consistently declined over time at high initial ODs (e.g., 10⁻² and 10⁻⁴) across all replicates (Fig. 4a–f), suggesting that Cα and Cβ robustly outperformed Ks. Conversely, when the initial OD was low (10⁻⁶ and 10⁻⁸), the consortium exhibited two divergent modes of succession (Fig. 4c–d, g–h): the Ks fraction declined monotonically as in the high initial OD cases; alternatively, it increased over time and dominated the population. Furthermore, comparison of all four cases suggested that the chance of Ks dominance increased with reducing initial OD. To directly visualize such succession outcomes, we further collected the co-cultures at the end of fermentations and performed colorimetric assays by adding X-Gluc to the co-cultures (Methods). Because Ks encodes beta-glucuronidase which can produce a clear blue green color, the colors of the treated supernatants (Fig. 4i) reflected the Ks dominance in individual experiments. Linking to the characterization of sampling-induced variations (Fig. 2e, Supplementary Tables 7, 8), these results confirmed that increasing sampling-induced variations equally drives the consortium into divergent outcomes. Supporting the statement, we computationally turned off the sampling-induced variations in our mathematical model and found that, without the variations, Ks persistently declined regardless of initial ODs (Fig. 4a–d, bold lines). The theory-experiment consistency suggested that incorporating variability into ecosystem modeling provides a predictive capacity over community behaviors.

**Fig. 4: Divergent ecosystem dynamics directed by cooperation variations from sampling.**

Construction of pH-dependent competition and cooperation

To examine if variability-modulated ecosystem succession is general to different cellular interactions, we designed a new cooperation-competition consortium composed of the cooperators—Cα and Cβ—and Kp, a strain that opposes Cα and Cβ and resists their killing in a pH-dependent manner. Using L. lactis MG1363 as the host, the pH-dependent Kp-to-Cα/Cβ inhibition was created by applying a pH-inducible promoter P₇₇₄⁶⁴ to control lcnA, the precursor gene of the lcnA pathway and the constitutive promoters P₁ and P₂ to drive the rest genes (lceA, lcmA, and lciA) in the pathway (Fig. 5a). Similarly, the pH-dependent resistance to Cα/Cβ-to-Kp inhibition was enabled by using the promoter P₇₇₄ to drive lagC, the immunity gene of lcnG co-produced by Cα and Cβ. Additionally, gusA3⁶³ was constitutively expressed to enable colorimetric quantification of the strain.

**Fig. 5: A synthetic strain whose interactions vary with environmental pH.**

Previous studies showed that the promoter P₇₇₄ is active when the environmental pH is below 6.5 but switched to be inactive when above 7⁶⁴. Thus, Kp’s lcnA production (i.e., inhibition over Cα and Cβ) and LagC production (i.e., resistance to killing by Cα and Cβ) are no longer constant but, instead, vary with the environment. Importantly, these strains are all derived from L. lactis which naturally produces a large amount of lactic acid and, thus, can lower the pH of culture in fermentation, which suggests that both the Kp-to-Cα/Cβ and Cα/Cβ-to-Kp inhibitions can be highly dynamic even in simple batch fermentation.

To validate the interactions, we grew Kp monoculture under three settings: pH ≥ 7, pH ≤ 6, and no pH control (Methods) (Fig. 5b and Supplementary Fig. 13). Our results (Fig. 5c and Supplementary Fig. 14) showed that the size of the lcnA inhibition zones (Methods) remains undetectable during the pH ≥ 7 fermentation, suggesting no lcnA production. By contrast, in the pH ≤ 6 fermentation, lcnA was detected as early as 3 h after fermentation and the culture yielded the highest lcnA level. For the case of no pH control, lcnA was detected after 4 h of fermentation and eventually accumulated to a medium level. These results confirmed that lcnA production (Kp-to Cα/Cβ inhibition) is highly correlated with the environmental pH. To confirm the pH-dependence of the Cα/Cβ-to-Kp inhibition, we cultured Kp and Ks in pH-defined media mixed with the supernatant of the corresponding Cα-Cβ or Cα-Cβ’ co-cultures (Cβ’ is a Cβ variant deficient in β subunit production), and compared their relative growth rates (Methods, Supplementary Fig. 9e, f and 15). The results (Fig. 5d) showed that Kp grew better than Ks in all three conditions despite a higher load due to lagC production, demonstrating that Kp gained a resistance to lcnG. Additionally, the results (Fig. 5d, blue bars) showed that Kp grew best at pH ≤ 6 and worst at pH ≥ 7, confirming pH modulates Kp’s resistance to lcnG. Equivalently, the results suggested that both the growth of Cα and Cβ monoculture (Supplementary Fig. 16) and the strength of Cα-Cβ cooperation are indeed modulated by pH (Supplementary Fig. 17).

These results further implied that, subject to environmental pH, the Cα-Cβ-Kp ecosystem can exhibit multiple modes of interaction. When pH ≥ 7, Cα and Cβ cooperatively inhibit Kp (Fig. 6a); when pH ≤ 6, Kp inhibits Cα and Cβ (Fig. 6b); by contrast, when there is no pH control, Cα and Cβ inhibit Kp initially but, later, are suppressed by Kp (Fig. 6c).

**Fig. 6: Ecosystem dynamics jointly regulated by cooperation variations and competition variations.**

Dynamics jointly regulated by multiple interaction variations

To illustrate how the Cα-Cβ-Kp consortium evolves upon both cooperation and competition variations, we assembled a dynamic community model (Supplementary Equation 17) from the bottom up (Methods). Then we used the model to explore the succession of the consortium when its interaction strengths vary due to simultaneous pH and cooperator partition alterations. In parallel, we performed Cα-Cβ-Kp co-culture experiments under the conditions identical to the computational test (Methods).

For the consortium starting from a fixed initial (Cα + Cβ):Kp ratio (2:1), the model predicted and subsequent experiments confirmed that, when pH was controlled above 7, Kp abundance declined gradually over fermentation for different initial Cα:Cβ ratios but reached to the lowest at 1:1 (Fig. 6d, g). In contrast, when pH was below 6, Kp became increasingly dominant over time for all Cα:Cβ ratios but augmented the most at unbalanced cases (30:1 and 1:30) (Fig. 6e, h). When there was no pH control, the consortium succession displayed two distinct patterns: At unbalanced ratios, Kp evolved to be dominant as the case of pH ≤ 6 but, at close ratios, Kp abundance declined over time as the pH ≥ 7 case (Fig. 6f, i).

Although seemingly diversified, these successions can be elucidated by considering the multiple interaction variations caused by pH modulation and labor partition. First, the systematic shift of the Kp abundance from consistent decrease (Fig. 6d, g) to consistent increase (Fig. 6e, h) and divergent development (Fig. 6f, i) originated from the pH-induced variation of the interactions: At pH ≥ 7, Kp-to-Cα/Cβ inhibition was abolished but Cα/Cβ-to-Kp inhibition remained potent; at pH ≤ 6, Kp-to-Cα/Cβ inhibition became effective (Supplementary Fig. 18) but Cα/Cβ-to-Kp inhibition was significantly reduced; when there was no pH control, Cα and Cβ inhibited Kp at beginning but were later suppressed by Kp (Fig. 6a–c). Second, within a single pH setting, the final Kp abundance was lower at close Cα:Cβ partitions (e.g. 1:1) than at imbalanced (e.g., 30:1 and 1:30) because Cα and Cβ had a stronger lcnG productivity when their partitions are close (Supplementary Fig. 17). These results showed that for ecosystems containing multiple variable interactions, at least for those we tested, their succession is determined jointly by all of the variations but not by any one of them. Of note, the divergent dynamics in the absence of pH control (Fig. 6f, i) exemplified the superposition of pH and labor partition effects: At imbalanced Cα:Cβ partitions, Cα and Cβ had the potential to kill Kp at beginning but their lcnG yield was too low; later, Kp gained the lcnG resistance and further secreted lcnA to kill Cα and Cβ, leading to the monotonic increase of Kp abundance. In contrast, at the close partitions, Cα and Cβ produced significant lcnG to efficiently inhibit Kp during the initial fermentation and the lcnG remained in the culture continued to suppress Kp even though the interaction topology was later altered.

To further demonstrate this finding, we conducted additional assays for the consortium by varying pH and initial total OD. Here, the model and the experimental setups were the same as previous except for the initial conditions (Methods). Accordingly, the sources of variation became pH modulation and random sampling. Our results showed that, at pH ≥ 7, the Kp abundance consistently reduced at high initial ODs (Fig. 7a–f), owing to the cooperative inhibition of Cα and Cβ to Kp; However, at low initial ODs, it could also remain largely invariant in some replicates (Fig. 7c–d, g–h) since increasing variations at random sampling abolished the cooperation (Fig. 2e–f). By contrast, in the absence of sampling-induced variations, the Kp abundance always declined regardless of initial ODs (Fig. 7a–d, bold lines). At pH ≤ 6, Kp became increasingly dominant regardless of initial ODs (Fig. 7i–l, m–p), because Kp constitutively suppressed Cα and Cβ and such a suppression was not affected by the fluctuations of Cα:Cβ ratios. These results are consistent with the ecosystem succession when sampling-induced variations were eliminated (Fig. 7i–l, bold lines). When there was no pH control, Kp declined minorly to a plateau at high initial ODs (Fig. 7q–r, u–v), attributed to the factors that Cα and Cβ collaborated to suppress Kp initially but were later suppressed by Kp due to pH reduction. At low initial ODs, it declined as the high initial OD case or diverged to be dominant (Fig. 7s–t, w–x) because increasing randomness diminished the Cα/Cβ-to-Kp inhibition at the beginning but did not affect Kp-to-Cα/Cβ suppression later. For comparison, in the absence of variation, there was no Kp dominance under all initial conditions (Fig. 7q–t, bold lines). These results demonstrated again that it is the joint regulation from multiple interaction variations that determines the dynamics of the ecosystems.

**Fig. 7: Model-predicted and experimentally measured succession of the Cα-Cβ-Kp ecosystem for varied pH conditions and initial densities.**

To quantitatively evaluate the capacity of the variability-incorporated modeling scheme, we calculated relative errors, defined as the differences between simulation and experimental results divided by experimental measures, for all of the simulations and experimental data above. The results (Supplementary Fig. 19) showed that, for most of the comparisons, the simulations agree quantitatively with experimental measures with the mean absolute relative errors (MARE) falling within the range of (0, 0.2) (panels a–d, g). For a subset of the cases, the models have larger relative errors but yet qualitatively agree well with the experimental findings (panels e, f, h). The encouraging agreements between the simulations and experiments demonstrated that incorporating interaction variability into ecosystem modeling is a promising strategy for quantitative and predictive understanding of complex community behaviors. Meanwhile, the discrepancies in certain cases suggest that the current models may need to consider additional processes involved in the experimental ecosystems, such as the nonlinearity observed during parameter fitting in Supplementary Fig. 8, in order to achieve a better modeling-experiment agreement.

Discussion

Microbial interactions are often modeled invariant²⁰; however, in nature, they constantly fluctuate over time and such fluctuations in strength are shown to be profound to ecosystem behaviors^{17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32}. Using synthetic microbial consortia as simple and reliable platforms, we showed that increasing variations of interaction diversifies an ecosystem’s succession into distinct outcomes. We also showed that, when there are multiple variable interactions^18,19,25, these variations collectively, but not a single one, regulate the behaviors of a community. Together, our results established interaction variability as a critical modulator of ecosystem behaviors.

Our synthetic ecosystems are relatively simple as they contain only a competition and a cooperation. Yet, their dynamics can be dramatically modulated by the variations of interaction. In native ecologies, there are significantly more cellular interactions and more microbial and environmental factors that modulate cellular interactions. We thus speculate that interaction variability is not specific to our synthetic systems and potentially a universal determinant for microbial ecosystem succession.

Searching for assembly rules has been invaluable to our understanding of community organization^65,66,67. In our study, we observed the diversification of ecosystem succession with increasing variations, which illustrates the intrinsic complexity and context-dependence nature of community behaviors. This finding implies that, in certain scenarios, qualitative rules may not be sufficient to specify the assembly of a community; instead, quantitative and systematic characterizations are needed. The result might also partially explain the current difficulty in predicting microbial ecosystem behaviors, an existing key challenge in microbial ecology⁴.

Lastly, our study shows encouraging consistence between experiment and mathematical modeling, which illustrates the promise of variability-cast ecosystem modeling for bottom-up predictions of the structure and dynamics of microbial ecosystems. Together, this work provides fundamental insights into the organization of microbial communities and also the de novo engineering of microbial consortia for various biotechnological applications.

Methods

Strain and growth conditions

All strains are derived from L. lactis MG1363 and grown at 30 °C in M17 broth supplemented with 0.5% (w/v) glucose and 5 μg mL⁻¹ of chloramphenicol (GM17/Cm). Tween 80 was added at a final concentration of 0.1% (v/v) when necessary. Cell cultures are adjusted with 2 M NaOH solution every hour to maintain pH above 7. To achieve pH ≤ 6, cell cultures are adjusted by 1 M HCl every hour. Strains used in this study are described in Supplementary Table 1.

Plasmid construction

All plasmids used in this study were developed from a L. lactis-E. coli shuttle vector, pleiss-Nuc⁶⁸, and described in Supplementary Table 1. Oligos for plasmid construction are listed in Supplementary Table 2. To generate the plasmid pleiss-lcnG for lactococcin G production, a 5-kb fragment of lcnG gene cluster including lagA, lagB, lagC, lagD, and lagE was amplified from the genome of L. lactis LMG 2081 using primers of lcnG-F and lcnG-R⁵³, and subsequently assembled with a fragment of pleiss-Nuc amplified with primers P_g-F and P_g-R using Gibson assembly. The plasmid pleiss-lcnG was then transformed into L. lactis MG1363 to obtain the lcnG producing strain Cαβ. To construct lcnG subunit expression plasmids pleiss-Cα₀ and pleiss-Cβ₀, the α or β coding gene was deleted from the plasmid pleiss-lcnG by reverse PCR and Gibson assembly using two pairs of primers: Cα₀-F/Cα₀-R and Cβ₀-F/Cβ₀-R. The resulting plasmids were transformed into L. lactis MG1363 to generate strains Cα₀ and Cβ₀. To enable screening and counting of cells with different subunits, a gfp or rfp reporter gene, yemGFP or mCherry was introduced using primers P_cα-F/P_cα-R and Cα-F/Cα-R, generating the plasmids pleiss-Cα and pleiss-Cβ. These plasmids were subsequently transformed into L. lactis MG1363 to construct the reporter version of α and β peptide producer Cα and Cβ. As a control, plasmid pleiss-Cβ’ was generated by deleting the β precursor gene lagB from the plasmid pleiss-Cβ with primers P_Cβ’-F, P_Cβ’-F, Cβ’-F and Cβ’-R. The resulting plasmid was then transformed into L. lactis MG1363 to get a β-free variant Cβ’. To increase the productivity of α peptide in Cα₀, an additional copy of expression cassette of α under the control of P₄ promoter was inserted to pleiss-Cα₀ using primers P_cα0-P4α-F/ P_cα0-P4α-R and Cα₀-P_4a-F/ Cα₀-P_4a-R⁶⁹, generating the variant Cα_0v2. Cβ_0v1/3 and Cβ_0v1/5, with approximately 1/3 and 1/5 of the productivity of Cβ₀ respectively, were created through inserting simple short repeat sequence (AT)_n into the spacer region of ribosome binding site of β peptide to weaken the translational initiation rate with primers P_Cβ0v1/3-F/ P_Cβ0v1/3-R, Cβ_0v1/3-F/Cβ_0v1/3-R, P_Cβ0v1/5-F/ P_Cβ0v1/5-R, Cβ_0v1/5-F/Cβ_0v1/5-R⁷⁰. The lactococcin A producing plasmid pleiss-lcnA was constructed by assembling the lcnA gene cluster from the plasmids pFI2396 and pFI2148 with pleiss-Nuc vector⁶². To simplify the detection of the lcnA-producing strain, a reporter gene gusA3 was amplified from the plasmid pTRK892⁶³ and then inserted into pleiss-lcnA using primers P_lcnA-gusA3-F/P_lcnA-gusA3-R and lcnA-gusA3-F/lcnA-gusA3-R. The resulting plasmid pleiss-lcnA-gusA3 was transformed into L. lactis MG1363 to generate the lcnA-producing strain Ks. Kr, as a control strain for Ks with immunity to lcnG and a gusA3 reporter, was created by assembling the immunity gene lagC of lcnG and GusA3 using the primers of P_IG-GusA3-F/P_IG-GusA3-R and GusA3-F/GusA3-R. To create the strain with pH-dependent lcnA production and lcnG resistance, plasmid pleiss-P₇₇₄-lcnA-gusA3 was firstly created by replacing the lcnA’s native promoter P_lcna of the plasmid pleiss-lcnA-gusA3 with a pH-inducible promoter P₇₇₄⁶⁴. Subsequently, pleiss-P₇₇₄-lcnA-gusA3 was assembled with the fragment of P₇₇₄-lagC using two pairs of primers: P_p774-lcnA-F and P_p774-lcnA-R, P₇₇₄-lcnA-F and P₇₇₄-lcnA-R; P_p774-lagC-F and P_p774-lagC-R, P₇₇₄-lagC-F and P₇₇₄-lagC-R, resulting in the final plasmid pleiss-P₇₇₄-lcnA-P₇₇₄-lagC-gusA3. The final plasmid was then transformed into L. lactis MG1363 to generate the strain Kp.

Measurement of lcnG productivity

The agar diffusion assay was performed using a protocol adapted from a previous study⁵⁵. Specifically, cultures of lcnG producing cells were grown in GM17/Cm/Tween broth at 30 °C overnight under corresponding culture conditions. The overnight cultures were inoculated in fresh media at 1:50 dilution and grown to the early stationary phase. Supernatants were obtained by centrifuging at 10,000× g for 10 min. Then, 30 μL samples were added into the wells in a double-layer agar (15 mL of GM17/Cm/Tween with 0.75% agar for each layer) in which the bottom layer was seeded with 50 μL of overnight culture of inducer strain L. lactis NZ9000/pleiss-Nuc. After incubation at 30 °C for 8 h, the inhibition zones were characterized by the blank circles around wells. To establish a standard curve of relative lcnG concentration for quantitatively evaluating samples, the concentration of lcnG in the supernatant of Cαβ culture at the early stationary phase is defined as 100%. And then the cell-free supernatant of Cαβ culture was diluted with fresh GM17/Cm/Tween media to the relative lcnG concentrations of 90%, 80%, 70%, 60%, 50%, 40%, 30%, 20%, 10%, 5%, 2%, and 1%. Next, 30 μL of samples were added into the wells in the top layer agar. After incubation at 30 °C for 8 h, the inhibition zones emerged and a standard curve was drawn by measuring the diameters of inhibition zones produced by different relative concentrations of lcnG. Using this curve, the relative concentrations of lcnG from tested samples were estimated.

Characterization of variations from labor partition

Monocultures of Cα and Cβ (Cα₀ and Cβ₀) were inoculated in GM17/Cm/Tween liquid medium at 30 °C overnight. Then the optical density at 600 nm (OD) of monocultures were measured and the co-cultures of Cα and Cβ (Cα₀ and Cβ₀) were mixed at a variety of ratios (30:1, 10:1, 3:1, 1:1, 1:3, 1:10, 1:30) with a start total OD₆₀₀ of 10⁻². After 8 h, the co-culture cells were centrifuged at 10,000× g for 10 min to obtain the co-culture supernatants. To test the cooperative strength in supernatant mixtures, overnight monocultures were firstly transformed to fresh GM17/Cm/Tween broth at the initial OD₆₀₀ of 10⁻² individually. Then supernatants were extracted after growth for 8 h by centrifuging and mixed at different ratios as mentioned above. Finally, 30 μL of samples were used to determine the activity of lcnG by agar diffusion assays to examine the variability of cooperation in co-culture supernatants and monoculture supernatant mixtures. Additionally, to evaluate the role of subunit stoichiometric ratio in cooperation variation, monoculture supernatants were also prepared from Cα₀, Cα_0v2, Cβ_0v1/3, and Cβ_0v1/5 as described above and mixed at ratios from 30:1, 10:1, 3:1, 1:1, 1:3, 1:5, 1:6, 1:10 to 1:30. Among them, ratios of 1:5 and 1:6 were used for determining the optimal cooperativity between Cα₀ and Cβ_0v1/5 (Cα_0v2 and Cβ_0v1/3). Then, activities of lcnG in these combinations were determined by agar diffusion assay. Similarly, to measure the cooperative inhibition from Cα and Cβ to Ks (Kp), three steps were involved. First, starting media were prepared by mixing 2-fold concentrated GM17/Cm/Tween medium with equal volume of supernatants from Cα-Cβ co-culture at the early stationary phase that were filtered with sterile 0.22 μm filter and adjusted to necessary pH. Second, overnight Ks (Kp) monoculture was inoculated with an initial OD of 10⁻² to the starting media. Third, the relative abundances of Ks were measured with fluorescence microscope from the starting to 8 h of incubation, during which pH was controlled properly if needed. For comparison, starting media were also prepared by mixing 2-fold concentrated GM17/Cm/Tween medium with equal volume of supernatants of Cα-Cβ’ co-culture for culturing and measuring Ks (Kp). Each experimental condition has three replicates (n = 3).

Determination of the maximum dilution rate

To determine the minimum OD₆₀₀ that is required for cell growth in a fresh medium, monocultures of Cα₀, Cβ₀, Cα, Cβ and Ks were serially diluted to an OD₆₀₀ of 10⁻² to 10⁻¹³ at 1:10 dilution and the survival rate of each OD₆₀₀ was calculated by counting the growing cultures in twenty replicates of each OD₆₀₀. After incubation at 30 °C for 24 h, the tubes with cell growth were counted to calculate the survival rate of different initial OD₆₀₀. Our experiment (Supplementary Fig. 5) showed that 10⁻⁸ is the minimum OD₆₀₀ with a 100% survival rate for almost all strains. An inoculation with an initial OD₆₀₀ lower than 10⁻⁸ could result in failure in growth and would disturb the studies of small number fluctuations in the community. Therefore, the initial OD₆₀₀ after dilution in the serial dilution experiments should be higher than 10⁻⁸.

Characterization of cooperation variations from sampling

To examine the effect of sampling on the variability of Cα-Cβ cooperation, overnight cultures of Cα and Cβ were washed twice with sterile PBS buffer (pH = 7) and re-suspended in PBS buffer. Then, the Cα and Cβ suspensions were adjusted to an OD₆₀₀ of 1.0 with PBS buffer, and mixed together at 1:1 ratio. The resulting suspension was used as a start culture and diluted to the OD₆₀₀ of 10⁻², 10⁻⁴, 10⁻⁶ and 10⁻⁸ through serial 1:10 dilutions with PBS buffer. The total volume of a start culture was set at 5 mL. After being prepared, all samples were centrifuged at 10,000× g for 15 min to remove supernatants, and 5 mL of fresh GM17/Cm media were subsequently added for cell growth. At the end of incubation, the supernatants were obtained. Subsequently, inhibition zone assays were conducted to determine the strength of cooperation. As cultures with different initial ODs require different incubation times, for each initial condition we chose sampling time based on the corresponding growth profile so that the ODs at each time point are comparable across the samples (Supplementary Fig. 12). Such samplings enable a consistent and proper comparison. Ten trials were performed for each initial condition.

Three-strain cooperator-varying experiments

The initial total OD₆₀₀ of three strains Cα, Cβ and Ks (Kr or Kp) was set at 10⁻². The start abundance of Ks (Kr or Kp) was fixed at 33.3% in the population but the ratios of Cα and Cβ were set at 30:1, 10:1, 3:1, 1:1, 1:3, 1:10, and 1:30. During incubation, samples were taken every two hours for measuring their ODs and the relative numbers of green (Cα), red (Cβ), and non-fluorescent cells were counted under an AMG EVOS FL fluorescence microscope using green, red and bright field channels. Notably, for the pH-controlled three-strain experiments, the overnight culture of Kp was washed twice with sterile PBS buffer (pH = 7) and then inoculated with an OD of 10⁻² into fresh medium with proper different pH controls (pH ≥ 7, pH ≤ 6 and no pH control). During the course of fermentation, culture samples were collected every hour to determine environmental pH using pH meter and Kp’s lcnA productivity using the inhibition zone assay. Each culture condition has three replicates (n = 3). The pH and inhibition zones were calculated as mean ± s.d.

Three-strain sampling experiments

For the three-strain ecosystem of Cα, Cβ, and Ks (Kp), Cα and Cβ PBS suspensions were individually adjusted to OD₆₀₀ of 1 and then mixed at 1:1 ratio; Ks (Kp) PBS suspension was also adjusted to OD₆₀₀ of 1. Subsequently, both the Cα-Cβ mixture and Ks (Kp) were diluted to the OD₆₀₀ of 10⁻², 10⁻⁴, 10⁻⁶, and 10⁻⁸ with PBS buffer in 1:10 dilution fold. At each concentration, the Cα-Cβ mixture and Ks (Kp) monoculture suspensions were mixed at a 2:1 ratio to reach a final ratio of 1:1:1 for Cα, Cβ and Ks (Kp). The total volumes of both monoculture and co-culture suspensions were set at 5 mL. The prepared samples were centrifuged at 10,000× g for 15 min to remove supernatants, and 5 mL of fresh GM17/Cm media were subsequently added for cell growth. The co-cultures were then incubated at 30 °C and propagated for 8 (10⁻²), 9 (10⁻⁴), 14 (10⁻⁶) and 20 (10⁻⁸) hours which was needed for entering the early stationary phase. To consistently determine cell numbers in different initial ODs (e.g., 10⁻², 10⁻⁴, 10⁻⁶, 10⁻⁸), the initial relative abundances of Cα, Cβ and Ks (Kp) were measured by colony forming unit (CFU) counting instead of flow cytometry, because flow cytometry typically requires a minimum cell number (10⁴ cells per mL) in samples. For the samples with an initial OD₆₀₀ of 10⁻² or 10⁻⁴, the mixed suspensions were diluted to around OD₆₀₀ of 5×10⁻⁶ with sterile PBS buffer and 1 mL of the diluted suspension was plated on the GM17/Cm plate supplemented with X-Gluc for counting Ks (Kp). For the samples with an initial OD₆₀₀ of 10⁻⁶, 1 mL of the mixed suspensions were directly plated without dilution on the GM17/Cm plate supplemented with X-Gluc. For the samples with an initial OD₆₀₀ of 10⁻⁸, 5 mL of the mixed suspensions were plated without dilution on the GM17/Cm plate supplemented with X-Gluc. During the period of incubation, the samples were taken to measure OD₆₀₀ and the fraction of Ks (Kp) strain was calculated by counting the cells under microscope. For the ecosystem of Cα, Cβ and Kp, the above procedure remained the same except that the diluted co-culture samples were inoculated into pH-defined media whose pH was adjusted every hour. Ten trials were performed for each initial condition.

Quantification of lcnA inhibition

To measure the inhibition from Ks to Cα or Cβ, overnight Ks and Cα or Cβ monoculture were inoculated into fresh GM17/Cm/Tween media at 1:1 ratio based on OD₆₀₀. Then, the relative abundances of Ks were measured from the starting to 8 h of incubation with fluorescence microscope. The similar procedure was also used to determine the inhibition from Kp to Cα and Cβ in pH-defined media (pH ≥ 7, pH ≤ 6 and no pH control). Each culture condition has three replicates (n = 3).

GusA3 enzyme assay

To directly visualize the Ks in three-strain system, GusA3 protein was used as a reporter to produce a blue green color. At the end of growth in three-strain consortium, a final concentration of 2 mM of 5-bromo-4-chloro-3-indolyl- β-D-glucuronide (X-Gluc) was added into the cultures and the blue green color formation was monitored. After incubation at 37 °C for 1 h, the supernatants were obtained by centrifuging the cultures at 10,000× g for 10 min and images were taken.

Mathematical modeling

In concert with the experimental ecosystem assembly, a bottom-up strategy was utilized to construct the dynamic models of the synthetic ecosystems. Briefly, we first created and characterized growth models of monocultures (Cα, Cβ, Ks, Kr, Kp), then models for the cooperative Cα-Cβ species as well as their interactions, and finally assembled individual modules into integrated models that represent complete ecosystems (e.g., Cα-Cβ-Ks, Cα-Cβ-Kr and Cα-Cβ-Kp consortia). During model construction, ordinary differential equations were used to quantitatively describe the kinetics of three major classes of variables: nutrient availability, cell populations and bacteriocins mediating cellular interactions. Parameters were determined using data in the literature or by fitting the models to our experiments. MATLAB software was used to simulate the models, produce plots, and fit data for the models. A detailed description of the models is available in Supplementary Information.

Statistical analysis

All of the experiments were performed for multiple times. Replicate numbers of the experiments (n) are indicated in the figure legends. Sample sizes were chosen based on standard experimental requirement in molecular biology. Data are presented as mean ± s.d.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability

Strains and plasmids constructed in this study are available from the corresponding author upon request. Data supporting the results in this paper are available within the paper and its supplementary information files. The source data of Figs. 2b,c,e, 3c,d,g,h, 4a-h, 5b-d, 6d-6i, 7 and Supplementary Figs. 1a, 2c, 3, 5, 7b and 8-19 are provided as a Source Data file. All other relevant data are available from the author upon reasonable request.

Code availability

Custom MATLAB codes developed in this study are available from the corresponding author upon request.

References

Falkowski, P. G., Fenchel, T. & Delong, E. F. The microbial engines that drive Earth’s biogeochemical cycles. Science 320, 1034–1039 (2008).
Article ADS CAS PubMed Google Scholar
Berendsen, R. L., Pieterse, C. M. & Bakker, P. A. The rhizosphere microbiome and plant health. Trends Plant Sci. 17, 478–486 (2012).
Article CAS Google Scholar
Cho, I. & Blaser, M. J. The human microbiome: at the interface of health and disease. Nat. Rev. Genet. 13, 260–270 (2012).
Article CAS PubMed PubMed Central Google Scholar
Widder, S. et al. Challenges in microbial ecology: building predictive understanding of community function and dynamics. ISME J. 10, 2557–2568 (2016).
Article PubMed PubMed Central Google Scholar
Brenner, K., You, L. & Arnold, F. H. Engineering microbial consortia: a new frontier in synthetic biology. Trends Biotechnol. 26, 483–489 (2008).
Article CAS PubMed Google Scholar
Johns, N. I., Blazejewski, T., Gomes, A. L. & Wang, H. H. Principles for designing synthetic microbial communities. Curr. Opin. Microbiol. 31, 146–153 (2016).
Article PubMed PubMed Central Google Scholar
Großkopf, T. & Soyer, O. S. Synthetic microbial communities. Curr. Opin. Microbiol. 18, 72–77 (2014).
Article PubMed PubMed Central CAS Google Scholar
Mimee, M., Citorik, R. J. & Lu, T. K. Microbiome therapeutics - advances and challenges. Adv. Drug Deliv. Rev. 105, 44–54 (2016).
Article CAS PubMed PubMed Central Google Scholar
Song, H. S., Cannon, W. R., Beliaev, A. S. & Konopka, A. Mathematical modeling of microbial community dynamics: A methodological review. Processes 2, 711–752 (2014).
Article Google Scholar
Cordero, O. X. & Datta, M. S. Microbial interactions and community assembly at microscales. Curr. Opin. Microbiol. 31, 227–234 (2016).
Article PubMed PubMed Central Google Scholar
Gonze, D., Coyte, K. Z., Lahti, L. & Faust, K. Microbial communities as dynamical systems. Curr. Opin. Microbiol. 44, 41–49 (2018).
Article PubMed Google Scholar
Zomorrodi, A. R. & Segrè, D. Synthetic ecology of microbes: mathematical models and applications. J. Mol. Biol. 428, 837–861 (2016).
Article CAS PubMed Google Scholar
Faust, K. & Raes, J. Microbial interactions: from networks to models. Nat. Rev. Microbiol. 10, 538–550 (2012).
Article CAS PubMed Google Scholar
Mounier, J. et al. Microbial interactions within a cheese microbial community. Appl. Environ. Microbiol. 74, 172–181 (2008).
Article CAS PubMed Google Scholar
Marino, S., Baxter, N. T., Huffnagle, G. B., Petrosino, J. F. & Schloss, P. D. Mathematical modeling of primary succession of murine intestinal microbiota. Proc. Natl. Acad. Sci. USA 111, 439–444 (2014).
Article ADS CAS PubMed Google Scholar
Stein, R. R. et al. Ecological modeling from time-series inference: insight into dynamics and stability of intestinal microbiota. PLoS Comput. Biol. 9, e1003388 (2013).
Article PubMed PubMed Central CAS Google Scholar
Garcia, F. C., Bestion, E., Warfield, R. & Yvon-Durocher, G. Changes in temperature alter the relationship between biodiversity and ecosystem functioning. Proc. Natl Acad. Sci. USA 115, 10989–10994 (2018).
Article PubMed CAS Google Scholar
Ratzke, C. & Gore, J. Modifying and reacting to the environmental pH can drive bacterial interactions. PLoS Biol. 16, e2004248 (2018).
Article PubMed PubMed Central CAS Google Scholar
Rivett, D. W. et al. Resource-dependent attenuation of species interactions during bacterial succession. ISME J. 10, 2259–2268 (2016).
Article PubMed PubMed Central Google Scholar
Hart, S. F. M. et al. Uncovering and resolving challenges of quantitative modeling in a simplified community of interacting cells. PLoS Biol. 17, e3000135 (2019).
Article CAS PubMed PubMed Central Google Scholar
LaSarre, B., McCully, A. L., Lennon, J. T. & McKinlay, J. B. Microbial mutualism dynamics governed by dose-dependent toxicity of cross-fed nutrients. ISME J. 11, 337–348 (2017).
Article CAS PubMed Google Scholar
Bachmann, H., Molenaar, D., Kleerebezem, M. & van Hylckama Vlieg, J. E. High local substrate availability stabilizes a cooperative trait. ISME J. 5, 929–932 (2011).
Article CAS PubMed Google Scholar
Liu, J., Wu, C., Huang, I. H., Merritt, J. & Qi, F. Differential response of Streptococcus mutans towards friend and foe in mixed-species cultures. Microbiology 157, 2433–2444 (2011).
Article CAS PubMed PubMed Central Google Scholar
Nguyen, A. T., Jones, J. W., Ruge, M. A., Kane, M. A. & Oglesby-Sherrouse, A. G. Iron depletion enhances production of antimicrobials by Pseudomonas aeruginosa. J. Bacteriol. 197, 2265–2275 (2015).
Article CAS PubMed PubMed Central Google Scholar
Kreth, J., Merritt, J., Shi, W. & Qi, F. Competition and coexistence between Streptococcus mutans and Streptococcus sanguinis in the dental biofilm. J. Bacteriol. 187, 7193–7203 (2005).
Article CAS PubMed PubMed Central Google Scholar
Liu, F., Mao, J., Lu, T. & Hua, Q. Synthetic, context-dependent microbial consortium of predator and prey. ACS Synth. Biol. 8, 1713–1722 (2019).
Article CAS PubMed Google Scholar
Andrade-Domínguez, A., Salazar, E., Vargas-Lagunas, M., Kolter, R. & Encarnación, S. Eco-evolutionary feedbacks drive species interactions. ISME J. 8, 1041–1054 (2014).
Article PubMed Google Scholar
Buffie, C. G. et al. Precision microbiome reconstitution restores bile acid mediated resistance to Clostridium difficile. Nature 517, 205–208 (2015).
Article ADS CAS PubMed Google Scholar
Palmer, A. C., Angelino, E. & Kishony, R. Chemical decay of an antibiotic inverts selection for resistance. Nat. Chem. Biol. 6, 105–107 (2010).
Article CAS PubMed PubMed Central Google Scholar
van der Putten, W. H., Klironomos, J. N. & Wardle, D. A. Microbial ecology of biological invasions. ISME J. 1, 28–37 (2007).
Article PubMed Google Scholar
Kuipers, O. P., Beerthuyzen, M. M., de Ruyter, P. G., Luesink, E. J. & de Vos, W. M. Autoregulation of nisin biosynthesis in Lactococcus lactis by signal transduction. J. Biol. Chem. 270, 27299–27304 (1995).
Article CAS PubMed Google Scholar
Coburn, P. S., Pillar, C. M., Jett, B. D., Haas, W. & Gilmore, M. S. Enterococcus faecalis senses target cells and in response expresses cytolysin. Science 306, 2270–2272 (2004).
Article ADS CAS PubMed Google Scholar
Kaern, M., Elston, T. C., Blake, W. J. & Collins, J. J. Stochasticity in gene expression: from theories to phenotypes. Nat. Rev. Genet. 6, 451–464 (2005).
Article CAS PubMed Google Scholar
Elowitz, M. B., Levine, A. J., Siggia, E. D. & Swain, P. S. Stochastic gene expression in a single cell. Science 297, 1183–1186 (2002).
Article ADS CAS PubMed Google Scholar
Kiviet, D. J. et al. Stochasticity of metabolism and growth at the single-cell level. Nature 514, 376–379 (2014).
Article ADS CAS PubMed Google Scholar
Ramkrishna, D. & Mahoney, A. W. Population balance modeling. Promise for the future. Chem. Eng. Sci. 57, 595–606 (2002).
Article CAS Google Scholar
Shu, C. C., Chatterjee, A., Dunny, G., Hu, W. S. & Ramkrishna, D. Bistability versus bimodal distributions in gene regulatory processes from population balance. PLoS Comput. Biol. 7, e1002140 (2011).
Article MathSciNet CAS PubMed PubMed Central Google Scholar
Hellweger, F. L. & Bucci, V. A bunch of tiny individuals—Individual-based modeling for microbes. Ecol. Model 220, 8–22 (2009).
Article Google Scholar
Klitgord, N. & Segre, D. Environments that induce synthetic microbial ecosystems. PLoS Comput. Biol. 6, e1001002 (2010).
Article ADS PubMed PubMed Central CAS Google Scholar
Resat, H., Bailey, V., McCue, L. A. & Konopka, A. Modeling microbial dynamics in heterogeneous environments: Growth on soil carbon sources. Microb. Ecol. 63, 883–897 (2012).
Article CAS PubMed Google Scholar
Shou, W., Ram, S. & Vilar, J. M. Synthetic cooperation in engineered yeast populations. Proc. Natl. Acad. Sci. USA 104, 1877–1882 (2007).
Article ADS CAS PubMed Google Scholar
Gore, J., Youk, H. & van Oudenaarden, A. Snowdrift game dynamics and facultative cheating in yeast. Nature 459, 253–256 (2009).
Article ADS CAS PubMed PubMed Central Google Scholar
Balagaddé, F. K. et al. A synthetic Escherichia coli predator-prey ecosystem. Mol. Syst. Biol. 4, 187 (2008).
Article PubMed PubMed Central Google Scholar
Ozgen, V. C., Kong, W., Blanchard, A. E., Liu, F. & Lu, T. Spatial interference scale as a determinant of microbial range expansion. Sci. Adv. 4, eaau0695 (2018).
Article ADS Google Scholar
Mee, M. T., Collins, J. J., Church, G. M. & Wang, H. H. Syntrophic exchange in synthetic microbial communities. Proc. Natl. Acad. Sci. USA 111, E2149–E2156 (2014).
Article ADS CAS PubMed Google Scholar
Kerr, B., Riley, M. A., Feldman, M. W. & Bohannan, B. J. Local dispersal promotes biodiversity in a real-life game of rock-paper-scissors. Nature 418, 171–174 (2002).
Article ADS CAS PubMed Google Scholar
Chuang, J. S., Rivoire, O. & Leibler, S. Simpson’s paradox in a synthetic microbial system. Science 323, 272–275 (2009).
Article CAS PubMed Google Scholar
Wintermute, E. H. & Silver, P. A. Emergent cooperation in microbial metabolism. Mol. Syst. Biol. 6, 407 (2010).
Article PubMed PubMed Central Google Scholar
Youk, H. & Lim, W. A. Secreting and sensing the same molecule allows cells to achieve versatile social behaviors. Science 343, 1242782 (2014).
Article PubMed PubMed Central CAS Google Scholar
West, S. A., Griffin, A. S., Gardner, A. & Diggle, S. P. Social evolution theory for microorganisms. Nat. Rev. Microbiol. 4, 597–607 (2006).
Article CAS PubMed Google Scholar
Wingreen, N. S. & Levin, S. A. Cooperation among microorganisms. PLoS Biol. 4, 1486–1488 (2006).
Article CAS Google Scholar
Celiker, H. & Gore, J. Cellular cooperation: insights from microbes. Trends Cell Biol. 23, 9–15 (2013).
Article CAS PubMed Google Scholar
Mirkovic, N. et al. Lactococcus lactis LMG2081 produces two bacteriocins, a non lantibiotic and a novel lantibiotic. Appl. Environ. Micobiol 82, 2555–2562 (2016).
Article CAS Google Scholar
Moll, G. et al. Lactococcin G is a potassium ion-conducting, two-component bacteriocin. J. Bacteriol. 178, 600–605 (1996).
Article CAS PubMed PubMed Central Google Scholar
Kong, W., Kapuganti, V. S. & Lu, T. A gene network engineering platform for lactic acid bacteria. Nucleic Acids Res. 44, e37 (2016).
Article PubMed Google Scholar
Bentley, W. E., Mirjalili, N., Andersen, D. C., Davis, R. H. & Kompala, D. S. Plasmid-encoded protein: The principal factor in the “metabolic burden” associated with recombinant bacteria. Biotechnol. Bioeng. 35, 668–681 (1990).
Article CAS PubMed Google Scholar
Glick, B. R. Metabolic load and heterologous gene expression. Biotechnol. Adv. 13, 247–261 (1995).
Article CAS PubMed Google Scholar
Lidstrom, M. E. & Konopka, M. C. The role of physiological heterogeneity in microbial population behavior. Nat. Chem. Biol. 6, 705–712 (2010).
Article CAS PubMed Google Scholar
Stegen, J. C., Lin, X., Konopka, A. E. & Fredrickson, J. K. Stochastic and deterministic assembly processes in subsurface microbial communities. ISME J. 6, 1653–1664 (2012).
Article CAS PubMed PubMed Central Google Scholar
Masel, J. Genetic drift. Curr. Biol. 21, R837–R838 (2011).
Article CAS PubMed Google Scholar
Hallatschek, O., Hersen, P., Ramanathan, S. & Nelson, D. R. Genetic drift at expanding frontiers promotes gene segregation. Proc. Natl Acad. Sci. USA 104, 19926–19930 (2007).
Article ADS CAS PubMed Google Scholar
Fernández, A., Horn, N., Gasson, M. J., Dodd, H. M. & Rodríguez, J. M. High-level coproduction of the bacteriocins nisin A and lactococcin A by Lactococcus lactis. J. Dairy Res. 71, 216–221 (2004).
Article PubMed CAS Google Scholar
Callanan, M. J., Russell, W. M. & Klaenhammer, T. R. Modification of Lactobacillus β-glucuronidase activity by random mutagenesis. Gene 389, 122–127 (2007).
Article CAS PubMed Google Scholar
Madsen, S. M., Arnau, J., Vrang, A., Givskov, M. & Israelsen, H. Molecular characterization of the pH–inducible and growth phase–dependent promoter P170 of Lactococcus lactis. Mol. Microbiol. 32, 75–87 (1999).
Article CAS PubMed Google Scholar
Cody, M. L. & Diamond, J. M. Ecology and evolution of communities. Nature 260, 204 (1976).
Article Google Scholar
Gotelli, N. J. How do communities come together? Science 286, 1684–1685 (1999).
Article CAS Google Scholar
Weiher, E. & Keddy, P. Ecological assembly rules: perspectives, advances, retreats. Cambridge University Press (2001).
Le Loir, Y., Gruss, A., Ehrlich, S. D. & Langella, P. A nine-residue synthetic propeptide enhances secretion efficiency of heterologous proteins in Lactococcus lactis. J. Bacteriol. 180, 1895–1903 (1998).
Article PubMed PubMed Central Google Scholar
Zhu, D. et al. Isolation of strong constitutive promoters from Lactococcus lactis subsp. lactis N8. FEMS Microbiol. Lett. 362, fnv107 (2015).
Article PubMed CAS Google Scholar
Egbert, R. G. & Klavins, E. Fine-tuning gene networks using simple sequence repeats. Proc. Natl. Acad. Sci. USA 109, 16817–16822 (2012).
Article ADS CAS PubMed Google Scholar

Download references

Acknowledgements

We thank Prof. Jon Nissen-Meyer and Prof. Todd Klaenhammer for their generous gifts of DNA materials. This work was supported by the National Science Foundation (1553649), Office of Naval Research (N000141612525) and Department of Energy (DE-SC0019185). F.L. was supported by the China Scholarship Council. J.M. was supported by the China Scholarship Council and the National Natural Science Foundation of China (11575059 and 11847315). Q.H. was supported by the Research Program of the State Key Laboratory of Bioreactor Engineering and 111 Project (B18022).

Author information

These authors contributed equally: Feng Liu, Junwen Mao

Authors and Affiliations

Department of Bioengineering, University of Illinois at Urbana-Champaign, Urbana, IL, USA
Feng Liu, Junwen Mao, Wentao Kong, Rashid Bashir & Ting Lu
Carl R. Woese Institute for Genomic Biology, University of Illinois at Urbana-Champaign, Urbana, IL, USA
Feng Liu, Junwen Mao, Wentao Kong & Ting Lu
State Key Laboratory of Bioreactor Engineering, East China University of Science and Technology, Shanghai, China
Feng Liu & Qiang Hua
School of Bioengineering, East China University of Science and Technology, Shanghai, China
Feng Liu & Qiang Hua
Department of Physics, Huzhou University, Huzhou, China
Junwen Mao
School of Medicine, Zhejiang University, Hangzhou, China
Youjun Feng
Micro and Nanotechnology Laboratory, University of Illinois at Urbana-Champaign, Urbana, IL, USA
Rashid Bashir
Carle Illinois College of Medicine, Urbana, IL, USA
Rashid Bashir
Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL, USA
Ting Lu
Center for Biophysics and Quantitative Biology, University of Illinois at Urbana-Champaign, Urbana, IL, USA
Ting Lu
National Center for Supercomputing Applications, University of Illinois at Urbana-Champaign, Urbana, IL, USA
Ting Lu

Authors

Feng Liu
View author publications
You can also search for this author in PubMed Google Scholar
Junwen Mao
View author publications
You can also search for this author in PubMed Google Scholar
Wentao Kong
View author publications
You can also search for this author in PubMed Google Scholar
Qiang Hua
View author publications
You can also search for this author in PubMed Google Scholar
Youjun Feng
View author publications
You can also search for this author in PubMed Google Scholar
Rashid Bashir
View author publications
You can also search for this author in PubMed Google Scholar
Ting Lu
View author publications
You can also search for this author in PubMed Google Scholar

Contributions

T.L. conceived the project; T.L., Q.H., Y.F. and R.B. designed research; F.L. performed experiments; J.M. developed mathematical models; W.K. contributed new reagents/analytic tools; F.L., J.M. and T.L. analyzed data; F.L., J.M. and T.L. wrote the paper with input from all other authors.

Corresponding author

Correspondence to Ting Lu.

Ethics declarations

Competing interests

The authors declare no competing interests.

Additional information

Peer review information Nature Communications thanks the anonymous reviewers for their contribution to the peer review of this work.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Supplementary Information

Reporting Summary

Source data

Source Data File

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Cite this article

Liu, F., Mao, J., Kong, W. et al. Interaction variability shapes succession of synthetic microbial ecosystems. Nat Commun 11, 309 (2020). https://doi.org/10.1038/s41467-019-13986-6

Download citation

Received: 22 March 2019
Accepted: 09 December 2019
Published: 16 January 2020
DOI: https://doi.org/10.1038/s41467-019-13986-6

This article is cited by

Advances in Microbe-Driven Ecological Restoration on Bauxite Residue Disposal Areas
- Shengguo Xue
- Xinfeng Qin
- Feng Zhu
Reviews of Environmental Contamination and Toxicology (2024)
Nitrogen transfer and cross-feeding between Azotobacter chroococcum and Paracoccus aminovorans promotes pyrene degradation
- Xia Wang
- Ying Teng
- Yongming Luo
The ISME Journal (2023)
Microbial Biofilms for Environmental Bioremediation of Heavy Metals: a Review
- Zainab Syed
- Monika Sogani
- Kumar Sonu
Applied Biochemistry and Biotechnology (2023)
Integrative Circuit-Host Modeling of a Genetic Switch in Varying Environments
- Jordan J. Sickle
- Congjian Ni
- Ting Lu
Scientific Reports (2020)

Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.