Microbial diversity and ecological complexity emerging from environmental variation and horizontal gene transfer in a simple mathematical model

Abstract

Background

Microbiomes are generally characterized by high diversity of coexisting microbial species and strains, and microbiome composition typically remains stable across a broad range of conditions. However, under fixed conditions, microbial ecology conforms with the exclusion principle under which two populations competing for the same resource within the same niche cannot coexist because the less fit population inevitably goes extinct. Therefore, the long-term persistence of microbiome diversity calls for an explanation.

Results

To explore the conditions for stabilization of microbial diversity, we developed a simple mathematical model consisting of two competing populations that could exchange a single gene allele via horizontal gene transfer (HGT). We found that, although in a fixed environment, with unbiased HGT, the system obeyed the exclusion principle, in an oscillating environment, within large regions of the phase space bounded by the rates of reproduction and HGT, the two populations coexist. Moreover, depending on the parameter combination, all three major types of symbiosis were obtained, namely, pure competition, host-parasite relationship, and mutualism. In each of these regimes, certain parameter combinations provided for synergy, that is, a greater total abundance of both populations compared to the abundance of the winning population in the fixed environment.

Conclusions

The results of this modeling study show that basic phenomena that are universal in microbial communities, namely, environmental variation and HGT, provide for stabilization and persistence of microbial diversity, and emergence of ecological complexity.

Background

The extensive recent efforts in metagenomics have led to the realization that microbiomes, for example, the well-studied gut microbial communities, are surprisingly stable on the species level over long periods of time [1,2,3]. Changes in microbiome content caused by such factors as antibiotics intervention or diet change can lead to obesity, inflammatory bowel disease, and other pathological conditions, so that microbiome stability is essential for human health [4]. Bacterial and archaeal species within the microbiome have metabolic networks that complement each other, and are thought to cooperate [5,6,7]. However, there are also indications of strong interspecies competition in complex microbial communities such as soil [8].

With the development of single cell metagenomics, finer structure of microbiomes was discovered, demonstrating that most of the constituent bacterial species are represented by multiple, closely related strains [9,10,11,12]. The coexistence of multiple strains has been shown to be common, and moreover, strain composition tends to be stable over extended periods of time [10, 13, 14]. It is generally assumed that the observed stability of the strain composition can be accounted for by models with no competition between strains [15]. Whereas prevalence of cooperation has been demonstrated for different microbial species, metabolic networks of strains of the same species are closely similar, so that competition for resources is expected to occur. Indeed, attempts on administration of probiotics have shown that bacteria, occupying the same niche as the species that is already dominant in the microbiome, are eliminated [16, 17].

A key process in the evolution of prokaryotes is the extensive horizontal gene transfer (HGT) that occurs via transformation, transduction, or conjugation (reviewed in [18]). The size of the transferred DNA fragments varies from several hundred nucleotides to segments of many kilobases encompassing multiple genes. Horizontal gene transfer is thought to be essential for the survival of microbial populations [19] and appears to be largely responsible for rapid adaptation to new environments and even for the emergence of major clades of bacteria and archaea with distinct lifestyles, such as acetoclastic methanogens or extreme halophiles [20, 21]. The rate of HGT is highly non-uniform across the tree of life and across environments [22]. Indeed, there is growing evidence that HGT rate can be influenced by the environmental changes. For example, it has been shown that natural competence of the bacterium Staphylococcus aureus is induced by reactive oxygen species, antibiotics, and host defenses [23, 24]. Importantly, the HGT rate is much higher between closely related organisms compared to distantly related ones thanks, primarily, to the high rate of homologous recombination [25,26,27].

In recent years, microbiome stability and variation have become a major problem in microbiology, with important implications for human health [28, 29]. Horizontal gene transfer plays a crucial role in the preservation of genome diversity [30,31,32,33,34,35,36] but the effects of environmental variations on the microbiome composition remain poorly understood. Evolutionary processes in populations in time-varying environments can drastically differ from those in fixed environments [37,38,39,40,41,42,43,44,45,46,47,48]. The complexity of microbial communities including multiple strains that persist for extended periods of time and widespread HGT among them call for developing theoretical models describing the relationships between strains within complex microbial communities, which are expected to compete with each other, but also constantly exchange genes via HGT.

Here we describe the simplest possible model that includes competition and HGT between two microbial populations differing by a single gene allele in a fixed or oscillating environment. The time-varying environment drastically changes the evolutionary scenarios observed in a fixed environment. We show that in this setting environmental variation plays a crucial role in the preservation of genome diversity. The coarse-graining modeling presented here allows us to uncover different types of symbiotic relationships that arise due to environmental variation, spanning all possible scenarios, from pure competition to mutualism [49,50,51,52].

Results

Environmental oscillations enable coexistence of competing populations

In our model, oscillations in the environment induce a new evolutionary game between the competing populations (Fig. 1; for a detailed description of the model, see the “Methods” section). The new game emerges as a result of coarse-graining of the effects of environmental variations on the competing populations. Specifically, the fast-oscillating terms induce an "effective potential," a new game that alters the slow-time dynamics of the fractions of the two populations and total abundance of both populations.

Fig. 1

Evolutionary dynamics of two competing populations with HGT in a fixed and in a fast-oscillating environment. The left part depicts the competition between two populations (blue and red circles, respectively) in the fixed (averaged) environment. The fractions of each population and the total abundance of both populations are described by the replicator dynamics with density dependent payoff matrices (the matrix on the left). Here, $r_{A}$ and $r_{B}$ are reproduction rates of the populations A and B in the fixed (averaged) environment, respectively. The competition within and between the two populations are given by $a_{ij}$, $i,j=A,B$ rates, respectively. N is the total abundance of both populations, and $\gamma$ is the gene transfer balance in the fixed (averaged) environment. In the balanced gene transfer case $\gamma =0$, the exclusion principle applies so that only blue circles are present at equilibrium. The right part depicts the competition in the oscillating environment. Environmental variations induce a new evolutionary game between the competing populations (the matrix on the right), in addition to the already existing one representing the interaction in the fixed (averaged) environment (the matrix on the left). The fraction of each population and the total abundance of both populations is, again, given by the replicator dynamics. However, in this case, the fitness of each population is comprised of two terms obtained from two different evolutionary games. $\xi$ and $\kappa$ represent environmental variations in the model which are defined by the product of oscillatory parts of the reproduction and gene transfer balance rates (16). Environmental oscillations enable the coexistence of both populations and can induce synergistic interaction between them, resulting in a greater total abundance of both populations at equilibrium compared to the abundance at equilibrium in the fixed (averaged) environment

The coarse-grained variations of the total population structure and abundance are defined by additive fitness obtained from two distinct games, one of which represents the evolutionary process in the fixed/averaged environment, and the other, emerging game represents the coarse-grained impact of the oscillating environment (16). The fitness of each population in the new game is expressed through the overlap in the oscillations of gene transfer balance rate $\tilde{\gamma }$ and the reproduction rate of the competitor over the period of environmental variations (for details, see the “Methods” section).

The new payoffs $\xi ,\kappa$ nullify when $\omega \rightarrow \infty$ for both the bounded oscillations of reproduction and gene transfer rates, as it follows from (16). That is, if environmental oscillations are too fast, then only the averages of reproduction and gene transfer balance determine the evolutionary outcome. The impact of the environmental variations is trivial if the oscillating part of the reproduction rates and gene transfer rate is given by the same time dependency, that is if $r_{i}(\tau ), \gamma (\tau ) \propto g(\tau )$, then $\xi ,\kappa =0$. Moreover, if only the reproduction rates $\tilde{r}_{i}$, $i=A,B$ or only gene transfer balance $\tilde{\gamma }$ oscillate, the environmental variations do not produce new dynamical effects compared to the competition outcome in fixed (averaged) environment. To substantially influence population growth and composition compared to the fixed environment, environmental variations must differentially affect the reproduction and gene transfer balance rates.

Environmental oscillation can result in coexistence of the two populations in a regime for which it is impossible in a fixed environment, that is, when gene transfer is balanced ($\gamma =0$) between competing populations and the allele difference affects only the reproduction rates ($a_{ij}=a$) (for details, see the “Methods” section). In a fixed/averaged environment, the exclusion principle applies so that only the population with the higher reproduction rate survives, giving the total abundance of both populations as $N_{f}^*= \frac{\textrm{max}[r_{A},r_{B}]}{a}$.

Assuming that $a_{ij}=a$ and gene transfer is balanced on average ($\gamma =0$), we substitute $a_{AB} \rightarrow a-\xi$, $a_{BA}\rightarrow a + \kappa$ in (11, 12) and find the following conditions for the existence and stability of the coexistence of the two populations in the coarse-grained dynamics (14, 15).

$$\begin{aligned}{} & {} \kappa <a \left( \frac{r_{B}}{r_{A}}-1\right) , \quad \xi >a \left( 1-\frac{r_{A}}{r_{B}}\right) ,\end{aligned}$$

(1)

$$\begin{aligned}{} & {} \left( 1-\frac{\xi }{a}\right) \left( 1+\frac{\kappa }{a}\right) <1. \end{aligned}$$

(2)

Let us assume that the average reproduction rates of both populations are almost equal, $\frac{r_{A}}{r_{B}}\gtrapprox 1$. Then, from (1), it follows that, in a varying environment coexistence is possible for $\xi >0$ and $\kappa <0$, whereas, from (11), it follows that in a fixed environment coexistence is impossible if $a_{ij}=a$ and $\gamma =0$. From (1) it follows that environmental variations dampen the between-population competition such that the competition within a population becomes more severe than the competition between the populations. This is the case because environmental variation shifts the between-population competition rates $a_{AB}=a-\xi <a$ and $a_{BA}=a+\kappa <a$, as follows from emerged fitness terms (16), thus reducing between-population competition compared to the within-population competition.

Synergistic interaction between the competing population can arise due to the environmental oscillations. By synergy, we refer to the situations when the presence of the two populations in the equilibrium state increases the total abundance of both populations compared to the total abundance at equilibrium in the fixed environment. The synergy emerges because the coexistence of the two populations increases the utilization of the available resources in the environment. Note that environmental variations alter the stability of equilibria with trivial compositions, that is, when only A or B population is present in the environment, but the total abundance of both populations is not affected by the environmental oscillations in these states. That is, the total abundance of both populations is at best equal to $N_{f}^*$ in these equilibrium states, whereas synergy emerges only in the coexistence regime. The total abundance of both populations in the coexisting equilibrium state, obtained from the coarse-grained representation of the oscillating environment, is

$$\begin{aligned} N_{c}^*=\frac{r_{B}\xi -r_{A}\kappa }{a(\xi -\kappa )+\xi \kappa } \end{aligned}$$

(3)

The conditions (1) and (2) provide $N_{c}^*>0$. The coarse-grained approach allows us to identify distinct sub-regions of different symbiotic relations between the two populations, that is, to define the impact of the competitor on the fitness of a population. These sub-regions are defined by the following expressions

$$\begin{aligned} \frac{\partial }{\partial p_{B}} (f_{A}+\phi _{A}) = -(a-\xi ) N, \qquad \frac{\partial }{\partial p_{A}} (f_{B}+\phi _{B}) = -(a+\kappa ) N, \end{aligned}$$

(4)

where $f_i$ and $\phi _i$, $i=A,B$, are the fitness in the fixed/averaged environment and the fitness surplus obtained due to the environmental oscillations of each population, respectively (for details, see the “Methods” section). $p_i$ and N are the fraction of each population and the total abundance of both populations, respectively.

Due to environmental oscillations, the derivatives (4) can change their signs in the case of time-dependent rates, hence, the type of the symbiotic relation between the two populations cannot be inferred without coarse-graining.

Assuming $r_{A}>r_{B}$, we now consider different regions of the parameters $\xi$ and $\kappa$, corresponding to different evolutionary outcomes.

The phase space of interactions between populations

Coexistence without synergy: pure competition and parasitism

This type of coexistence corresponds to region I in Fig. 2a, where $a \left( 1-\frac{r_{A}}{r_{B}}\right)<\xi<0\ \text {and}\ \kappa <a\left( \frac{r_{B}}{r_{A}}-1\right)$. In this region, the fitness of A, the winner in the fixed environment, is lower than it would be in the latter case. Indeed, the total fitness of A is $f_{A}+\phi _{A}$, where $\phi _{A}=N(1-p_{A})\xi <0$ (Fig. 1). Conversely, the population B obtains a fitness boost due to the environmental oscillations $\phi _{B}=-N p_{A} \kappa >0$. In the emerged game, B always has an advantage over A. Thus, A dominates over B in a fixed environment, but B has the advantage in the new game induced by the oscillating environment. In this region, there is coexistence but no synergy between two populations, that is, the total abundance of the two populations is lower than that in the fixed environment with same parameters, $N_{c}^*<N_{f}^*$.

Fig. 2

Phase space for possible outcomes of the coarse-grained description of two competing populations in an oscillating environment (14, 15), for different values of environment-induced fitness terms. a Partitioning of the phase space into distinct regions (see color code). Regions I, II, and III correspond to the coexistence of the two populations. In the regions II and III, there is synergy between the competing populations, that is, in the unique equilibrium state, the total abundance of both populations is greater than that in the fixed environment (for the considered values of the model parameters, A outcompetes B in the fixed environment). In regions IV, A outcompetes B, like in a fixed environment, V is the region of bistability, where either A or B wins depending on the initial state, and in region VI, B outcompetes A, opposite to the outcome in the fixed environment. Regions IV, V, and VI are defined by violations of one or more of the conditions (1, 2). N/A represents the region where the coarse-grained approach fails, that is (1) holds, but (2) is violated. b Detailed structure of the regions I, II, III, and N/A. The hatched areas represent the regions where the fitness of the two populations are differently affected by the presence of the competitor. The solid red curves show the total abundance of both population defined by (3), increasing from I to III. The model parameters are $r_{A}=1.8$, $r_{B}=1$, $\gamma =0$, and $a=0.1$

Region I includes two sub-regions that correspond to pure competition and parasitism. In the competitive sub-region (unhatched area of I in Fig. 2b), the fitness of each population is a decreasing function on the fraction of the competitor, that is both derivatives in (4) are negative. The derivative of the fitness of the B population changes its sign for $\kappa < -a$, so in the corresponding sub-region, B is a parasite of A ( hatched sub-region of I in Fig. 2b).

Synergy between populations: pure competition

This regime corresponds to region II in Fig. 2a, where ${0<\xi<a\ \text {and}\ -a<\kappa <a\left( \frac{r_{B}}{r_{A}}-1\right) }$. Both populations gain fitness due to the environmental variations, ${\phi _A,\phi _B>0}$, in contrast to region I. Thus, the two populations coexist due to the balance between the two games. The fitness of each population is a decreasing function of the competitor’s fraction, that is, both derivatives are negative in (4); hence, the interaction is purely competitive. Nevertheless, the total abundance of both populations is greater than that in the fixed environment $N_{c}^*>N_{f}^*$, thanks to the increased utilization of the available resources.

Synergy: parasitism and mutualism

This regime corresponds to region III (Fig. 2a, b). The two hatched sub-regions represent parasitic symbiosis between the competing populations A and B. In the vertical sub-region, ${\xi > a}$ and ${-a<\kappa <a\left( \frac{r_{B}}{r_{A}}-1\right) }$, the fitness of A and B are positively and negatively affected by the presence of the competitor, respectively, that is ${\frac{\partial }{\partial p_{B}} (f_{A}+\phi _{A})>0 \ \text {and}\ \frac{\partial }{\partial p_{A}} (f_{B}+\phi _{B})<0}$. Thus, A is a parasite of B. Conversely, in the horizontal sub-region, ${0<\xi <a}$ and ${\kappa < -a}$, fitness of A is negatively impacted by B, whereas A positively affects the fitness of B. Accordingly, the host-parasite relationship is reversed.

The symbiosis between the two populations is mutualistic when $\xi >a$ and $\kappa <-a$. In this sub-region (unhatched in Fig. 2b) of region III, the fitness of each population is an increasing function of the fraction of the other population $\frac{\partial }{\partial p_{j}} \left( f_{i}+\phi _{i}\right) >0$, $i\ne j=A, B$. The total abundance of both populations is greater than that in the fixed environment in both parasitic and mutualistic symbiosis, that is, both these types of symbiosis are synergistic.

The total abundance of both populations (3), $N_{c}^*=\textrm{const}$, increases from region I to region III (see Fig. 2b). The total abundance of both populations in the fixed environment corresponds to $\xi =0$ line; thus, by substituting $\xi =0$ in (3), one recovers $N_{c}^*=N_{f}^*=\frac{r_{A}}{a}$. In the fixed environment, only population A survives, so that the total abundance equals to the abundance of A, whereas in the oscillating environment, both A and B are present in equilibrium with fractions $p_{A}^*=\frac{a}{|\kappa |}\left( 1-\frac{r_{B}}{r_{A}}\right)$ and $p_{B}^*=1-p_{A}^*$, respectively.

Fig. 3

Fractions and total abundance of the two populations in the coexisting regime with various symbiotic relationships. a Fraction of the population A for different values of $\xi$ and $\kappa$. b Total abundance of both populations for different values of $\xi$ and $\kappa$. The oscillating curves show the solutions of (9, 10) with time-dependent reproduction and gene transfer balance rates. The smooth curves show coarse grained time variations of the respective quantities obtained through (14, 15). The blue, orange, and green curves show the time variations of the respective quantities for competition with no synergy $(\xi ,\kappa )=(-0.01,-0.08)$ (region I in Fig. 2b), competition with synergy $(\xi ,\kappa )=(0.08,-0.08)$ (region II), and for mutualistic symbiosis $(\xi ,\kappa )=(0.11,-0.11)$ (region III). The time-dependent rates are as follows: $\tilde{\gamma }(\tau )=0.5 \sin \tau$, $\tilde{r}_{A}(\tau )= c_{1} \cos \tau$, and $\tilde{r}_{B}(\tau )= c_{2} \cos \tau$, where $(c_{1},c_{2})= (-1.6,-0.2),~(-1.6,1.6)$ and $(-2.2, 2.2)$ for the blue, orange, and green curves, respectively. The remaining model parameters are $a=0.1$, $\gamma =0$, $\omega =5$, $r_{A}=1.8$ and $r_{B}=1$. The dashed line in (b) represents the total abundance of both populations in the fixed environment $N_{f}^*=\frac{r_{A}}{a}$

Figure 3 illustrates the solutions for the fractions and total abundance of two populations, for both time-dependent rates and the coarse-grained behavior (14, 15). The ordering of fractions and total abundance does not follow a uniform pattern across the different scenarios of interaction between the populations. In the case of pure competition with synergy, the equilibrium fraction of population A surpasses that in the mutualistic scenario. However, the total abundance of both populations is greater in the mutualistic scenario than in the pure competitive case. The scenario with no synergy (region I in Fig. 2) yields the smallest total abundance for both populations as well as, for the fraction of population A. These observations imply complex relationships between competition, synergy, and mutualistic symbiosis, which affect both population fractions and overall abundance.

So far, we addressed the emergence of coexistence between two competing populations due to fast environmental variations. In other parts of the phase space, however, the environmental variations can preserve the outcome of the competition in the fixed environment, that is, A outcompetes B (region IV in Fig. 2), or reverse it (region VI) resulting in the domination of B and extinction of A. The former scenario occurs if $\xi >a\left( 1-\frac{r_{A}}{r_{B}}\right)$ and $\kappa >a\left( \frac{r_{B}}{r_{A}}-1\right)$, that is, the inequality (1) is reversed for $\kappa$. Conversely, B outcompetes A if the inequality is reversed for $\xi$, whereas that for $\kappa$ holds. Finally, the bistable dynamics (region V) is observed if (1) and (2) are simultaneously violated. Here, both populations are engaged in the pure competition symbiosis, as follows from (4).

Unbalanced gene transfer ($\gamma \ne 0$) in a fixed environment can result in stable coexistence of the two populations (13), (for details, see the “Methods” section). In this case, environmental variations can alter the composition and the total abundance of the two populations and even destroy the coexistence of the two competing populations. Indeed, substituting $a_{AA}=a_{BB}=a$, $a_{AB}=a-\xi$ and $a_{BA}=a+\kappa$ in (11, 12), we obtain the conditions for the stable coexistence of both populations with unbalanced gene transfer in the oscillating environment:

$$\begin{aligned}{} & {} \gamma +\kappa <a\left( \frac{r_{B}}{r_{A}}-1\right) ,\quad \gamma +\xi >a\left( 1-\frac{r_{A}}{r_{B}}\right) ,\end{aligned}$$

(5)

$$\begin{aligned}{} & {} a(\xi -\kappa )+(\gamma +\xi )(\gamma +\kappa )>0 \end{aligned}$$

(6)

The conditions (5) and (6) are the counterparts of (1) and (2), respectively, for the case of unbalanced gene transfer in a fixed environment. From (5), (6), and (13), together, it follows that the environmental variations can change the composition and the total abundance of the two populations in the coexisting equilibrium. Environmental oscillations destroy the coexistence of the two populations observed in the fixed environment if any of the conditions (5, 6) is violated, but (13) holds.

The coarse-grained approach fails as one gets closer to the curve $a(\xi -\kappa )+\xi \kappa =0$, where $N_{c}^*\rightarrow \infty$, whereas the time-dependent solution is bounded (N/A regions in Fig. 2). Unbounded growth corresponds to violation of (2) that holds as equality. Below the curve defined by $a(\xi -\kappa )+\xi \kappa =0$, $N_{c}^*<0$. The fraction of each population satisfy $0<p_{A}^*,p_{B}^*<1$. Thus, there is a non-trivial composition of both populations, but there is no physical level of total abundance. This is one of the differences from the well-known replicator dynamics, which deals with the composition of the total population only. The reason behind the failure is the unbounded growth of both populations (see Additional file 1 for details on the limitations of the coarse-grained approach). The limitations of the model are discussed in the Additional file 1. There, we compare the solutions obtained by numerically integrating the system (9, 10) with time-dependent rates and the coarse grained counterparts obtained via (14, 15). The comparison is based on the Euclidean distance of the time-averaged fractions of both populations obtained via (9, 10) and those obtained from (14, 15). The difference is negligible for the regions IV, V, and VI, but not for the coexistence regimes. Nevertheless, even in the worst case scenario, both methods show the coexistence of both populations although the fractions of the two populations at equilibrium differ. We also compare the relative error of total population abundances. The relative error increases as one approaches the boundary $a(\xi -\kappa )+\xi \kappa =0$ where $N_{c}^*\rightarrow \infty$ (N/A), thus increasing the discrepancy between the results of the coarse-grained approach and the time average of the solution obtained via time-dependent rates, which remains bounded in general.

Discussion

In this work, we developed a mathematical model to explore the role of HGT in the interaction between two cohabiting populations in an oscillating vs a fixed environment. We explored what seems to be the simplest possible model in which two populations differed by a single gene allele, such that each entity in the pool carried only one allele of the given gene that is subject to HGT through copying the gene in the donor cell, and then, substituting the existing gene of the recipient cell, which corresponds to homologous recombination [14, 18]. Notwithstanding the ultimate simplicity of this scheme, it is biologically realistic, for example, in the context of acquisition of multidrug resistance [53]. Indeed, bacteria can adapt to changing environment by exogenous gene uptake and replacement of the already existing sequences with homologous ones [54,55,56].

We start with the stochastic description of all possible interactions between two populations, namely, reproduction and death of each population, competition within and between populations, and gene transfer between the populations. Then, we take the continuous limit and obtain the deterministic description of the time variations of the size of each population (see Additional file 1, [57, 58]). The dynamical system obtained by this procedure is a variation of the well-known replicator dynamics [59,60,61,62,63], with varying total abundance of both populations. The composition and the total abundance of both populations in the equilibrium states are found from the rest points of this system.

By examining this mathematical model, we show that stable coexistence of the two populations in a fixed environment is unattainable within the assumptions that the gene allele replacement via HGT only affects the reproduction rate, the competing abilities of both populations are the same within and between populations, and HGT between the populations is balanced, that is, there is no preferred direction in the gene flow between the populations. These results reflect the classic competitive exclusion (Gauze) principle according to which, in a spatially homogeneous environment, a population with an inferior reproduction rate will inevitably go extinct when competing with a fitter population for the same resource, within the same niche [64,65,66]. However, if the HGT is unbalanced, that is, the rates of gene transfer in the two directions are unequal, the exclusion principle does not apply anymore, and the two populations can coexist.

Evolutionary outcomes in a time-varying environment can drastically differ from those in the fixed environment, and despite the simplicity of our model, the emerging dynamics is complex, with the outcomes critically depending on the parameter combination. Environmental oscillations are incorporated into the model by assuming that the reproduction and gene transfer rates are time-periodic functions, with the implicit assumption that these rates are affected by the changes in the environment. The oscillations were set up such that the averages across environmental variations coincided with the constant rates in the fixed environment, providing for a fair comparison of the dynamics in the two types of environment.

The environmental oscillations occur on a fast time scale, whereby the populations are exposed to numerous changes in the environment before approaching any potential equilibrium. We adopted a coarse-graining approach to account for the environmental variations such that the solution for the composition and total abundance of both populations was obtained as a combination of two components, one delineating the slow-time (coarse grained) behavior and the other capturing the oscillating component [44, 45, 67]. By averaging over the period of environmental variations and retaining the first two leading-order terms of the varying quantities, we obtained the replicator dynamics model with varying total abundance that features fitness contributions derived from two distinct games. The first game represents the fixed (averaged) environment whereas the second one arises due to environmental fluctuations, with the payoffs in this game determined by the intersection of the oscillating components of the reproduction and gene transfer equilibrium rates. These new payoffs are such that the fitness of a particular population is defined by the average of the product of the competitor’s reproduction rate and the gene transfer equilibrium rate between the populations throughout environmental oscillations. For the emergence of the new game, reproduction and gene transfer balance rates must be differentially affected by the environment. Notably, variations either in the reproduction rates or in the gene transfer rate alone did not result in new evolutionary scenarios compared to the fixed environment case.

The emerged game maintains a straightforward payoff structure: once the environmental variations of reproduction and gene transfer rates are given, then either one of the strategies consistently dominates, or there exists a unique non-trivial equilibrium for any given total abundance of both populations. With these adjusted fitness terms, stable coexistence between the two populations becomes possible within a broad range of model parameters. Moreover, depending on the combination of the time-dependent reproduction and gene transfer balance rates, the co-existence of the two populations manifests as all major types of symbiotic relationships. Two populations can coexist in a purely competitive symbiosis, where the fitness of each is negatively impacted by the presence of the other population; a host-parasite relationship, whereby the impact of the second population is positive for one and negative for the other population; and in mutualistic symbiosis, where the presence of the other population is reciprocally beneficial.

We further analyzed the behavior of the total abundance of both populations and defined the regions of synergistic interactions between the competing or cooperating populations. In this case, synergy is observed, that is, the total abundance of the two populations at the stable coexistence in the oscillating environment is greater than the total abundance in the fixed environment, that is, the population size of the winner at equilibrium, under the same model parameter values. As could be expected, mutualism necessarily entails synergy and yields the greatest total abundance of both populations among all the regimes by increasing the resource utilization in the environment. However, under certain combinations of the reproduction and HGT rates, the synergistic effect was observed also in the cases of parasitism and even purely competitive symbiosis.

Outside of the coexistence regimes, the outcomes of the competition between two populations in a stable environment can persist in the oscillating environment such that the winner in the fixed environment wins in the oscillating environment as well. However, in another region of the phase space, the outcome can be reversed due to environmental oscillations. Moreover, there is also a bistability regime where one or the other population goes extinct depending on the initial conditions, a scenario precluded in the fixed environment assuming equal competing abilities and balanced gene transfer. Despite our prior discussion on coexistence and potential outcomes in an oscillating environment with balanced gene transfer in the fixed environment, the introduction of oscillations can notably alter the composition at the coexistence equilibrium, even if both populations coexist in the fixed environment due to non-zero gene transfer balance. Environmental variations, in this case, can even destroy the stable coexistence of the two populations, attained in the fixed environment with unbalanced HGT.

From the methods point of view, it is worth pointing out that coarse graining was an essential ingredient of the present analysis because without it, the symbiotic relationships between two populations would be impossible to elucidate because the fitness gradients can change their signs throughout environmental oscillation due to the oscillating rates.

From the biological standpoint, environmental oscillations provide for the synergy between the two populations by lowering the intensity of the inter-population competition below the level of the intra-population competition, thus, providing for stable coexistence of the two populations. This work shows that stabilization of strain diversity and increase of ecological complexity via HGT in an oscillating environment is an intrinsic feature of even the simplest microbiomes that emerges under minimal assumptions on the basic processes occurring within a microbial community. Given that both environmental variation and HGT are ubiquitous phenomena that affect any microbiome [68, 69], these conclusions appear to be broadly applicable. Notably, all emerging coexistence regimes, with different symbiotic relationships, provide for synergy between the populations, at least under certain parameter combinations, indicating that HGT in an oscillating environment is generally favorable for a microbial community perceived as an integral whole.

Evidently, the model used in this work is grossly (and deliberately) over-simplified. Interactions within microbiomes are highly complex even in a fixed environment [29, 70], and understanding the possible mechanisms stabilizing these communities is of great importance [71,72,73,74]. Real microbiomes encompass interactions among thousands of microbial strains and species, HGT of multiple genes via different routes and many other processes, resulting in extremely complex dynamics [14, 35, 36, 50, 69]. Nevertheless, the general principle established here should apply, amplified by the microbiome complexity. Characterization of the conditions for stabilization of microbiome diversity and the factors that can perturb it is crucial for understanding the role of the microbiome in health and diseases as well as the ecology of microbial communities.

Conclusions

In this work, we employed mathematical modeling to investigate how the diversity of microbiomes persists in the face of the exclusions principle that implies elimination of competing populations that share the same niche. We found that a simple mathematical model of competition between two populations that could exchange a single gene allele via HGT suggested potential solutions to this problem. In a fixed environment, with unbiased HGT, exclusion principle applied, and the fitter population drove the less fit one to extinction. By contrast, in an oscillating environment, in large domains of the phase space bounded by the reproduction rate and HGT, the two populations coexist. Depending on the parameter combinations, all three major types of symbiotic relationships were observed, namely, pure competition, parasitism relationship, and mutualism. Furthermore, in each of these regimes, certain parameter combinations resulted in synergy, that is, a greater total abundance of both populations compared to the abundance of the winning population in the fixed environment. Both environmental oscillations and HGT are ubiquitous in microbial communities. Thus, the results of this work suggest that these fundamental phenomena are both necessary and sufficient to ensure stabilization and persistence of microbial diversity and ecological complexity.

Methods

Model of microbial populations dynamics

We explore a dynamic evolutionary scenario where two populations engage in both within- and between -population competition for shared resources. The two populations differ by a single gene allele, resulting in differential fitness. HGT plays a key role in the model, mediating the exchange of genetic material between the two populations by copying the gene from the donor and replacing the corresponding gene in the recipient. If A is the donor and B is the recipient, then the allele from A replaces that in B, and as a result, the recipient changes its type $B\rightarrow A$, and the population sizes of A and B increase and decrease by one, respectively. The rates of gene transfer vary depending on whether A or B is the donor. Inactivation of genes by mutations is disregarded. This is a general approach for modeling HGT that does not depend on a particular mechanism (such as natural transformation, transduction or conjugation) as long as the new gene acquired via HGT replaces the existing counterpart and there is no spatial heterogeneity.

To capture the dynamics of this HGT-mediated interaction mathematically, we derive a system of differential equations that describes the continuous variation of population sizes over time:

$$\begin{aligned}{} & {} \frac{d n_A}{d t} = n_A(r_{A}-a_{AA} n_A -a_{AB} n_B) + \gamma n_A n_B,\end{aligned}$$

(7)

$$\begin{aligned}{} & {} \frac{d n_B}{d t} = n_B(r_{B}-a_{BA} n_A- a_{BB} n_B) -\gamma n_A n_B, \end{aligned}$$

(8)

where $n_A(t)$ and $n_B(t)$ are population sizes of A and B at time t, respectively. $r_{A}>0$ and $r_{B}>0$ are reproduction rates of A and B, respectively. $a_{ii}$ and $a_{ij}$ $i\ne j=A,B$ are within- and between-population competition rates. $\gamma$ reflects the balance of gene transfer between A and B types, that is, $\gamma >0$ corresponds to the positive gene flow from A to B such that, on average, B entities change their type to A, whereas for $\gamma <0$, A entities change to B. The rates describing the competition, both within and between two populations, are assumed to be positive $a_{ij}>0$, $i,j=A,B$, corresponding to a purely competitive ecosystem [65, 75], in the absence of HGT, $\gamma =0$. The detailed derivation of these equations from elementary processes, based on previous work [57, 58], is provided in the Additional file 1. For the purpose of the present work, we describe the system of (7, 8) by the equivalent representation through the total abundance of both populations $N=n_A+n_B$ and fractions of each population $p_{A}=n_A/N$, $p_{B}=n_B/N$, such that $p_{A}+p_{B}=1$. From (7, 8), we obtain the following dynamical system representing the time variations of total abundance of both populations and the fraction of each population:

$$\begin{aligned}{} & {} \frac{d p_{i}}{d t} = p_{i}\left( f_{i}-\sum \limits _{j}p_{j} f_{j}\right) ,\quad i,j=A,B\end{aligned}$$

(9)

$$\begin{aligned}{} & {} \frac{d N}{d t} = N \sum \limits _{j} p_{j} f_{j}. \end{aligned}$$

(10)

Equation (9) represents replicator dynamics well-known in evolutionary game theory [59,60,61,62,63]. However, the crucial difference between (9, 10) and standard replicator dynamics is the dependence of fitness functions $f_i$ on the total population size N. The fitness of each population can be represented as the expected payoff obtained from a $2 \times 2$ symmetric game with the payoff matrix shown in Fig. 1. The payoffs of each strategy/population A and B depend on the total abundance N of both populations. The time variation of the total abundance of both populations is defined by the mean fitness of both populations. In the equilibrium states, when the right-hand sides of both Eqs. (9, 10) nullify, the mean fitness and the fitness of each population represented in the equilibrium state nullify as well. Solving the right hand-side of (9, 10) with respect to the fractions of each population (taking into an account the normalization $p_{A}+p_{B}=1$) and total population abundance, we find 4 possible non-trivial outcomes: $\left(p_A^*,N^*\right) = \left( 1,\frac{r_{A}}{a_{AA}}\right)$ – extinction of population B, $\left( p_A^*,N^*\right) = \left( 0,\frac{r_{B}}{a_{BB}}\right)$ – extinction of population A, $\left( p_A^*,N^*\right) = \left( \frac{r_{A}a_{BB}-r_{B}(a_{AB}-\gamma )}{r_{A}(a_{BB}-a_{BA}-\gamma )+r_{B}(a_{AA}-a_{AB}+\gamma )},\frac{r_{A}(a_{BB}-a_{BA}-\gamma )+r_{B}(a_{AA}-a_{AB}+\gamma )}{a_{AA}{a_{BB}-(a_{AB}-\gamma )(a_{BA}+\gamma )}}\right)$– if $p_{A}^*\in (0,1)$, along with $N^*>0$, then depending on the stability of equilibrium the outcome is either coexistence of both populations or bistability, where extinction of a population depends on the initial condition, and the interior unstable equilibrium $p^*$ separates basins of attraction of two stable equilibrium states representing extinction of either A or B populations. The trivial equilibrium state, corresponding to the extinction of both populations, is always unstable for positive reproduction rates $r_{i}>0,~i=A,B$.

The stability of equilibrium states, where only one of the competing populations is present, is defined by the following conditions:

$$\begin{aligned} \frac{a_{BA}+\gamma }{a_{AA}} \frac{r_{A}}{r_{B}}>1 \quad \textrm{and} \quad \frac{a_{AB}-\gamma }{a_{BB}} \frac{r_{B}}{r_{A}}>1 \end{aligned}$$

(11)

for $\left( p_A^*,N^*\right) = \left( 1,\frac{r_{A}}{a_{AA}}\right)$ and $\left( p_A^*,N^*\right) = \left( 0,\frac{r_{B}}{a_{BB}}\right)$, respectively. For $\gamma =0$, the conditions (11) link between-populations and within-populations competition with the balance of the reproduction rates. Indeed for $r_{A}=r_{B}$, a given population outcompetes the opponent if the competition within the population of the winner is less severe than the competition between-populations.

If both conditions in (11) are satisfied simultaneously, along with $\frac{a_{AB}-\gamma }{a_{BB}}\frac{a_{BA}+\gamma }{a_{AA}}>1$, then one obtains a bistable dynamical system, where the outcome of competition depends on both the initial abundance and fractions of both populations.

Coexistence is obtained, once both conditions in (11) are violated simultaneously along with

$$\begin{aligned} \frac{a_{AB}-\gamma }{a_{BB}}\frac{a_{BA}+\gamma }{a_{AA}}<1 \end{aligned}$$

(12)

So far, we analyzed the model for general competition and reproduction rates, without imposing any other condition on the rates except for being positive. Therefore, our initial assumption is that populations A and B differ only by a single gene allele, and we further assume that the alleles of this gene affect only the reproduction rates $r_{A}$ and $r_{B}$, whereas the rates describing the competition within and between populations are independent of the gene allele variation and are equal, $a_{ij}=a$. Thus, for equal reproduction rates $r_{A}=r_{B}=r$ and balanced gene transfer between two populations $\gamma =0$, the total abundance of both populations is always equal to $N=\frac{r}{a}$. Coexistence of both populations could not be achieved in the case of balanced gene transfer $\gamma =0$ and unequal reproduction rates $r_{A}\ne r_{B}$, inasmuch as it is impossible to violate both conditions in (11) simultaneously. This impossibility of coexistence between populations corresponds to the exclusion principle which states that, if both populations compete for the same resources, then, one of them will be eliminated, under equal competition rates $a_{ij}=a$; in our setting, the exclusion principle applies in the case of balanced gene transfer $\gamma =0$ [64,65,66].

In the case of unbalanced gene transfer ($\gamma \ne 0$, coexistence of the two competing populations is possible if

$$\begin{aligned} a\left( 1-\frac{r_{A}}{r_{B}}\right)<\gamma < a\left( \frac{r_{B}}{r_{A}}-1\right) \end{aligned}$$

(13)

that is, both conditions in (11) are violated and (12) is satisfied simultaneously, under equal competition rates $a_{ij}=a$. Assuming that $r_{A}<r_{B}$, it follows from (13) that $\gamma >0$, that is, the coexistence of both population is possible if the slow reproducing population, A, has higher gene transfer rate, resulting in B entities change their type to A more frequently than A changes to B.

The above conclusions were obtained for a fixed environment. We now consider the same evolutionary process in a fluctuating environment. We assume that the environment oscillates faster than the characteristic timescale of the population dynamics, that is, both populations experience many changes of the environment before reaching the dynamic equilibrium. The environmental variations are incorporated in the model by assuming that the reproduction rates and the gene transfer balance rate are periodical function $r_{i}(\tau + 2\pi )=r_{i}(\tau ), ~i=A, B$ and $\gamma (\tau +2\pi )=\gamma (\tau )$ with period $2\pi /\omega$, where $\tau =\omega t, ~\omega >1$ describes fast variations of the environment. The time-dependent reproduction and gene transfer balance rates are represented as the sum of fixed and oscillatory parts $r_{i}(\tau )=r_{i}+\tilde{r}_{i}$ and $\gamma (\tau )=\gamma + \tilde{\gamma }$, where the time averages of oscillatory parts nullify over a period of environmental variations, that is $\overline{\tilde{r}_i}\equiv \int _0^{2\pi } \frac{d \tau }{2\pi } \tilde{r}_{i}=0 \ \text {and}\ \overline{\tilde{\gamma }}\equiv \int _0^{2\pi } \frac{d \tau }{2\pi } \tilde{\gamma }=0$. That is, $r_{i}$ and $\gamma$ represent the reproduction and gene transfer balance rates, respectively, in the averaged/fixed environment. Considering the evolutionary processes with averaged rates yields the system (9, 10) and the conclusions obtained for the fixed environment. We are interested in the coarse-grained behavior of the fractions and the total abundance of the competing populations in slow-time t. The coarse-graining procedure is based on the elimination of the fast-oscillating terms by identifying the dynamical impact of these terms on the slow-time behavior (the derivation of the slow-time behavior is given in the Additional file 1). The coarse-grained variations of the fraction of each population and total abundance of both populations, again, are given by the replicator dynamics

$$\begin{aligned}{} & {} \frac{d p_{i}}{d t} = p_{i}\left( f_{i}+\phi _{i}-\sum \limits _{j}p_{j} (f_{j}+\phi _{j})\right) ,\quad i,j=A,B\end{aligned}$$

(14)

$$\begin{aligned}{} & {} \frac{d N}{d t} = N \sum \limits _{j} p_{j} \left( f_{j}+\phi _{j}\right) . \end{aligned}$$

(15)

but, in contrast to (9, 10), the fitness of each population obtains an additional term due to the oscillating environment (Fig. 1). Note, that $f_{i}, i=A,B$ is given by the average/fixed rates, as in (9, 10). The new terms $\phi _{i}$ are expressed through the oscillating rates of reproduction and gene transfer balance as follows

$$\begin{aligned}{} & {} \phi _{A} = N(1-p_{A})\xi ,\quad \phi _{B}=-N\kappa p_{A}, \nonumber \\{} & {} \xi =\frac{1}{\omega }\int _0^{2\pi } \frac{d\tau }{2 \pi } \hat{r}_{B}\tilde{\gamma }, \quad \kappa = \frac{1}{\omega } \int _{0}^{2\pi } \frac{d\tau }{2 \pi } \hat{r}_{A}\tilde{\gamma }. \end{aligned}$$

(16)

In (16), $\hat{r}_{A}$ and $\hat{r}_{B}$ denote the primitives of $\tilde{r}_{A}$ and $\tilde{r}_{B}$, respectively, that is $\partial _\tau \hat{r}_i=\tilde{r}_i$. The existence and stability of possible equilibria of (14, 15) are described by the conditions (11, 12), where the between-population competition rates are substituted as follows $a_{AB} \rightarrow a_{AB}-\xi$ and $a_{BA}\rightarrow a_{BA} + \kappa$, for A and B populations, respectively.

The payoff structure of the emerged game is simple for any positive population abundance N (see Fig. 1). These payoffs are fixed once the time dependence of reproduction and gene transfer balance rate is known. Then, either there is a dominant strategy in the game (that is, either A or B always win in the new game)—if $\xi$ and $\kappa$ have the same sign (for example, A dominates B if $\xi , \kappa >0$), or there is a non-trivial equilibrium for any N if $\xi$ and $\kappa$ have different signs [60, 62, 76].