- Research article
- Open access
- Published:

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

*BMC Biology*
**volume 22**, Article number: 148 (2024)

## 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.

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).

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

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

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}^*\).

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.

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:

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:

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:

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:

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

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

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

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

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].

## Availability of data and materials

This is a theoretical work that did not involve data analysis.

## References

Chen L, Wang D, Garmaeva S, Kurilshikov A, Vich Vila A, Gacesa R, et al. The long-term genetic stability and individual specificity of the human gut microbiome. Cell. 2021;184(9):2302-2315.e12. https://doi.org/10.1016/j.cell.2021.03.024.

Mehta RS, Abu-Ali GS, Drew DA, Lloyd-Price J, Subramanian A, Lochhead P, et al. Stability of the human faecal microbiome in a cohort of adult men. Nat Microbiol. 2018;3(3):347–55. https://doi.org/10.1038/s41564-017-0096-0.

Huttenhower C, Gevers D, Knight R, Abubucker S, Badger JH, Chinwalla AT, et al. Structure, function and diversity of the healthy human microbiome. Nature. 2012;486(7402):207–14. https://doi.org/10.1038/nature11234.

Gilbert JA, Blaser MJ, Caporaso JG, Jansson JK, Lynch SV, Knight R. Current understanding of the human microbiome. Nat Med. 2018;24(4):392–400. https://doi.org/10.1038/nm.4517.

Freilich S, Zarecki R, Eilam O, Segal ES, Henry CS, Kupiec M, et al. Competitive and cooperative metabolic interactions in bacterial communities. Nat Commun. 2011;2(1):589. https://doi.org/10.1038/ncomms1597.

Levy R, Borenstein E. Metabolic modeling of species interaction in the human microbiome elucidates community-level assembly rules. Proc Natl Acad Sci. 2013;110(31):12804–9. https://doi.org/10.1073/pnas.1300926110.

Zelezniak A, Andrejev S, Ponomarova O, Mende DR, Bork P, Patil KR. Metabolic dependencies drive species co-occurrence in diverse microbial communities. Proc Natl Acad Sci. 2015;112(20):6449–54. https://doi.org/10.1073/pnas.1421834112.

Machado D, Maistrenko OM, Andrejev S, Kim Y, Bork P, Patil KR, et al. Polarization of microbial communities between competitive and cooperative metabolism. Nat Ecol Evol. 2021;5(2):195–203. https://doi.org/10.1038/s41559-020-01353-4.

Lo I, Denef VJ, VerBerkmoes NC, Shah MB, Goltsman D, DiBartolo G, et al. Strain-resolved community proteomics reveals recombining genomes of acidophilic bacteria. Nature. 2007;446(7135):537–41. https://doi.org/10.1038/nature05624.

Zhao S, Lieberman TD, Poyet M, Kauffman KM, Gibbons SM, Groussin M, et al. Adaptive evolution within gut microbiomes of healthy people. Cell Host Microbe. 2019;25(5):656-667.e8. https://doi.org/10.1016/j.chom.2019.03.007.

Bickhart DM, Kolmogorov M, Tseng E, Portik DM, Korobeynikov A, Tolstoganov I, et al. Generating lineage-resolved, complete metagenome-assembled genomes from complex microbial communities. Nat Biotechnol. 2022;40(5):711–9. https://doi.org/10.1038/s41587-021-01130-z.

Olm MR, Crits-Christoph A, Bouma-Gregson K, Firek BA, Morowitz MJ, Banfield JF. inStrain profiles population microdiversity from metagenomic data and sensitively detects shared microbial strains. Nat Biotechnol. 2021;39(6):727–36. https://doi.org/10.1038/s41587-020-00797-0.

Lloyd-Price J, Mahurkar A, Rahnavard G, Crabtree J, Orvis J, Hall AB, et al. Strains, functions and dynamics in the expanded Human Microbiome Project. Nature. 2017;550(7674):61–6. https://doi.org/10.1038/nature23889.

Van Rossum T, Ferretti P, Maistrenko OM, Bork P. Diversity within species: interpreting strains in microbiomes. Nat Rev Microbiol. 2020;18(9):491–506. https://doi.org/10.1038/s41579-020-0368-1.

Wolff R, Shoemaker W, Garud N. Ecological stability emerges at the level of strains in the human gut microbiome. Mbio. 2023;14(2):e02502-22.

Martínez I, Maldonado-Gomez MX, Gomes-Neto JC, Kittana H, Ding H, Schmaltz R, et al. Experimental evaluation of the importance of colonization history in early-life gut microbiota assembly. eLife. 2018;7:e36521. https://doi.org/10.7554/eLife.36521.

Segura Munoz RR, Mantz S, Martínez I, Li F, Schmaltz RJ, Pudlo NA, et al. Experimental evaluation of ecological principles to understand and modulate the outcome of bacterial strain competition in gut microbiomes. ISME J. 2022;16(6):1594–604. https://doi.org/10.1038/s41396-022-01208-9.

Arnold BJ, Huang IT, Hanage WP. Horizontal gene transfer and adaptive evolution in bacteria. Nat Rev Microbiol. 2022;20(4):206–18. https://doi.org/10.1038/s41579-021-00650-4.

Koonin EV. Horizontal gene transfer: essentiality and evolvability in prokaryotes, and roles in evolutionary transitions. F1000Res. 2016;5:F1000 Faculty Rev–1805. https://doi.org/10.12688/f1000research.8737.1.

Fournier GP, Gogarten JP. Evolution of acetoclastic methanogenesis in Methanosarcina via horizontal gene transfer from cellulolytic Clostridia. J Bacteriol. 2008;190(3):1124–7. https://doi.org/10.1128/jb.01382-07.

López-García P, Zivanovic Y, Deschamps P, Moreira D. Bacterial gene import and mesophilic adaptation in archaea. Nat Rev Microbiol. 2015;13(7):447–56. https://doi.org/10.1038/nrmicro3485.

Aminov RI. Horizontal gene exchange in environmental microbiota. Front Microbiol. 2011;2:158. https://doi.org/10.3389/fmicb.2011.00158.

Cordero M, García-Fernández J, Acosta IC, Yepes A, Avendano-Ortiz J, Lisowski C, et al. The induction of natural competence adapts staphylococcal metabolism to infection. Nat Commun. 2022;13(1):1525. https://doi.org/10.1038/s41467-022-29206-7.

Thi LTN, Romero VM, Morikawa K. Cell wall-affecting antibiotics modulate natural transformation in SigH-expressing Staphylococcus aureus. J Antibiot. 2016;69(6):464–6. https://doi.org/10.1038/ja.2015.132.

Fraser C, Hanage WP, Spratt BG. Recombination and the nature of bacterial speciation. Science. 2007;315(5811):476–80.

Dixit PD, Pang TY, Maslov S. Recombination-driven genome evolution and stability of bacterial species. Genetics. 2017;207(1):281–95.

Bobay LM, Ochman H. Biological species are universal across life’s domains. Genome Biol Evol. 2017;9(3):491–501.

Bucci V, Bradde S, Biroli G, Xavier JB. Social interaction, noise and antibiotic-mediated switches in the intestinal microbiota. PLoS Comput Biol. 2012;8(4):e1002497.

Coyte KZ, Schluter J, Foster KR. The ecology of the microbiome: networks, competition, and stability. Science. 2015;350(6261):663–6.

Novozhilov AS, Karev GP, Koonin EV. Mathematical modeling of evolution of horizontally transferred genes. Mol Biol Evol. 2005;22(8):1721–32.

Dohi M, Mougi A. A coexistence theory in microbial communities. R Soc Open Sci. 2018;5(9):180476.

Saakian DB, Koonin EV, Cheong KH. Key role of recombination in evolutionary processes with migration between two habitats. Phys Rev E. 2019;100(3):030401.

Takeuchi N, Cordero OX, Koonin EV, Kaneko K. Gene-specific selective sweeps in bacteria and archaea caused by negative frequency-dependent selection. BMC Biol. 2015;13:1–11.

Wang T, Weiss A, Aqeel A, Wu F, Lopatkin AJ, David LA, et al. Horizontal gene transfer enables programmable gene stability in synthetic microbiota. Nat Chem Biol. 2022;18(11):1245–52.

Pompei S, Bella E, Weitz JS, Grilli J, Lagomarsino MC. Metacommunity structure preserves genome diversity in the presence of gene-specific selective sweeps under moderate rates of horizontal gene transfer. PLoS Comput Biol. 2023;19(10):e1011532.

Coyte KZ, Stevenson C, Knight CG, Harrison E, Hall JP, Brockhurst MA. Horizontal gene transfer and ecological interactions jointly control microbiome stability. PLoS Biol. 2022;20(11):e3001847.

Levins R. Coexistence in a variable environment. Am Nat. 1979;114(6):765–83.

Chesson PL, Warner RR. Environmental variability promotes coexistence in lottery competitive systems. Am Nat. 1981;117(6):923–43.

Levins R. Evolution in changing environments: some theoretical explorations. Princeton University Press; 1968.

Burkart T, Willeke J, Frey E. Periodic temporal environmental variations induce coexistence in resource competition models. Phys Rev E. 2023;108(3):034404.

Steinmetz B, Meyer I, Shnerb NM. Evolution in fluctuating environments: a generic modular approach. Evolution. 2022;76(11):2739–57.

Babajanyan S, Koonin EV, Cheong KH. Can environmental manipulation help suppress cancer? Non-linear competition among tumor cells in periodically changing conditions. Adv Sci. 2020;7(16):2000340.

Babajanyan S, Lin W, Cheong KH. Cooperate or not cooperate in predictable but periodically varying situations? Cooperation in fast oscillating environment. Adv Sci. 2020;7(21):2001995.

Allahverdyan AE, Hu CK. Replicators in a fine-grained environment: adaptation and polymorphism. Phys Rev Lett. 2009;102(5):058102.

Allahverdyan AE, Babajanyan SG, Hu CK. Polymorphism in rapidly changing cyclic environment. Phys Rev E. 2019;100(3):032401.

Babajanyan S, Koonin E, Allahverdyan A. Thermodynamic selection: mechanisms and scenarios. New J Phys. 2022;24(5):053006.

Benincà E, Ballantine B, Ellner SP, Huisman J. Species fluctuations sustained by a cyclic succession at the edge of chaos. Proc Natl Acad Sci. 2015;112(20):6389–94.

Roca CP, Cuesta JA, Sánchez A. Time scales in evolutionary dynamics. Phys Rev Lett. 2006;97(15):158701.

Lee IPA, Eldakar OT, Gogarten JP, Andam CP. Bacterial cooperation through horizontal gene transfer. Trends Ecol Evol. 2022;37(3):223–32.

Palmer JD, Foster KR. Bacterial species rarely work together. Science. 2022;376(6593):581–2.

Kehe J, Ortiz A, Kulesa A, Gore J, Blainey PC, Friedman J. Positive interactions are common among culturable bacteria. Sci Adv. 2021;7(45):eabi7159.

Foster KR, Bell T. Competition, not cooperation, dominates interactions among culturable microbial species. Curr Biol. 2012;22(19):1845–50.

Perron GG, Lee AEG, Wang Y, Huang WE, Barraclough TG. Bacterial recombination promotes the evolution of multi-drug-resistance in functionally diverse populations. Proc R Soc B Biol Sci. 2012;279(1733):1477–84. https://doi.org/10.1098/rspb.2011.1933.

Szollosi GJ, Derényi I, Vellai T. The maintenance of sex in bacteria is ensured by its potential to reload genes. Genetics. 2006;174(4):2173–80.

Takeuchi N, Kaneko K, Koonin EV. Horizontal gene transfer can rescue prokaryotes from Muller’s ratchet: benefit of DNA from dead cells and population subdivision. G3: Genes, Genomes, Genet. 2014;4(2):325–39.

Redfield R. Genes for breakfast: the have-your-cake and-eat-it-too of bacterial transformation. J Hered. 1993;84(5):400–4.

McKane AJ, Newman TJ. Predator-prey cycles from resonant amplification of demographic stochasticity. Phys Rev Lett. 2005;94(21):218102.

McKane AJ, Newman TJ. Stochastic models in population biology and their deterministic analogs. Phys Rev E. 2004;70(4):041902.

Taylor PD, Jonker LB. Evolutionary stable strategies and game dynamics. Math Biosci. 1978;40(1–2):145–56.

Hofbauer J, Sigmund K. Evolutionary games and population dynamics. Cambridge University Press; 1998.

Sandholm WH. Population games and evolutionary dynamics. MIT Press; 2010.

Weibull JW. Evolutionary game theory. MIT Press; 1997.

Logofet D. Stability of biological communities. Mir Publishers; 1983.

Gause GF. The struggle for existence. Dover Books on Biology Series. Dover Publications, Incorporated; 2003. https://books.google.com/books?id=v01OToAhJboC.

Murray JD. Mathematical biology I: an introduction. New York: Springer; 2002.

Hardin G. The competitive exclusion principle: an idea that took a century to be born has implications in ecology, economics, and genetics. Science. 1960;131(3409):1292–7.

Abdul-Rahman F, Tranchina D, Gresham D. Fluctuating environments maintain genetic diversity through neutral fitness effects and balancing selection. Mol Biol Evol. 2021;38(10):4362–75.

Ferreiro A, Crook N, Gasparrini AJ, Dantas G. Multiscale evolutionary dynamics of host-associated microbiomes. Cell. 2018;172(6):1216–27.

Sarkar A, McInroy CJ, Harty S, Raulo A, Ibata NG, Valles-Colomer M, et al. Microbial transmission in the social microbiome and host health and disease. Cell. 2024;187(1):17–43.

Wang Y, Yang Y, Li A, Wang L. Stability of multi-layer ecosystems. J R Soc Interface. 2023;20(199):20220752.

Yang Y, Coyte KZ, Foster KR, Li A. Reactivity of complex communities can be more important than stability. Nat Commun. 2023;14(1):7204.

Yang Y, Foster KR, Coyte KZ, Li A. Time delays modulate the stability of complex ecosystems. Nat Ecol Evol. 2023;7(10):1610–9.

Allesina S, Tang S. Stability criteria for complex ecosystems. Nature. 2012;483(7388):205–8.

Mougi A, Kondoh M. Diversity of interaction types and ecological community stability. Science. 2012;337(6092):349–51.

Hallam TG, Levin SA. Mathematical ecology: an introduction, vol. 17. Springer Science & Business Media; 2012.

Robinson D, Goforth D. The topology of the 2x2 games: a new periodic table, vol. 3. Psychology Press; 2005.

## Acknowledgements

The authors are thankful to Armen E. Allahverdyan and Koonin group members for valuable discussions.

## Funding

Open access funding provided by the National Institutes of Health The authors are supported through the Intramural Research Program of the National Library of Medicine, National Institutes of Health.

## Author information

### Authors and Affiliations

### Contributions

SGB initiated the project; SGB, YIW, and EVK designed the model; SGB implemented the model; SGB, SKG, YIW, and EVK analyzed the results; SGB, SKG, and EVK wrote the manuscript that was read, edited, and approved by all authors.

### Corresponding authors

## Ethics declarations

### Ethics approval and consent to participate

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

## Additional information

### Publisher’s Note

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

## Supplementary information

## 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 licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence 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 licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

## About this article

### Cite this article

Babajanyan, S.G., Garushyants, S.K., Wolf, Y.I. *et al.* Microbial diversity and ecological complexity emerging from environmental variation and horizontal gene transfer in a simple mathematical model.
*BMC Biol* **22**, 148 (2024). https://doi.org/10.1186/s12915-024-01937-7

Received:

Accepted:

Published:

DOI: https://doi.org/10.1186/s12915-024-01937-7