Introduction

Microbial communities living in similar environments often display different community types, each characterized by a unique composition of species1-4. Changes in environmental factors, fluctuations of species abundances, or transient invasions can trigger abrupt and irreversible transitions between two types5-7. Such transitions, also known as regime shifts, lead to the dramatic alterations of microbiome functions and are tightly associated with many host diseases, soil fertility, marine biogeochemical flux, or bioproduction yield8-12. Understanding the origin of different community types has broad implications in environments, engineering and health1,3.

Multistability, a phenomenon where multiple stable states coexist for the same set of system parameters, is one of the most important mechanisms behind alternative community types and has been widely observed in a variety of natural or experimental microbiomes1,4,7,13-17. Under perturbations, a multistable system can be driven away from a former stable state, cross the tipping point and reach a new steady state. The population can persist in the new state even after the perturbations are withdrawn, exhibiting a behavior of ‘hysteresis’3,13. Despite its conceptual importance, the determinants underlying microbiota multistability have been poorly understood due to the intrinsic complexity of microbial interactions14.

Unraveling how multistability emerges in complex microbiota can provide crucial insights in different scenarios. For instance, the existence of alternative ‘hidden’ stable states in gut microbiome undermines the community resilience against external perturbations such as pathogenic infections or dietary changes18. Understanding which mechanisms contribute to multistability will help the design of modulatory strategies that prevent the developments of unhealthy stable states8,18. Multistability also renders the catastrophic shift of microbial populations under antibiotic exposure from being drug-sensitive to drug-resistant6. Inhibiting the processes underlying multistability can be useful to reverse antibiotic resistance and prolong the efficacy of current antibiotics19. In engineered microbiomes, the drastic compositional changes can result in the unexpected loss of the desired functions. The ability to predict and control community multistability is essential to better harness the bioproduction potential of synthetic microbial consortia15.

Mathematical models are indispensable to decipher the emergent properties of complex populations20. While different ecological models have been developed to analyze the influences of some factors, such as microbial mutual inhibition and nutrient stoichiometries, in mediating community multistability, the roles of many critical mechanisms have been overlooked in existing theories1,14. Chief among these is horizontal gene transfer (HGT), a process where microbes share mobile genetic elements (MGEs) with their neighbors21. Indeed, HGT is prevalent in microbial world and mobile genes occupy a substantial proportion in microbial gene pool22-24. For instance, among the closely related microbial isolates in the same human body site, HGT occurs in more than 40% of the species pairs25. In the draft genomes of human gut microbes, tens of thousands of genes were mobilizable23. However, despite its functional and ecological importance, the interplay between genetic exchange and community multistability has not been rigorously examined in previous studies and a general conclusion remains lacking.

In this work, we focus on communities of competing microbes, a model system for microbiome multistability1,14,26. By combining theoretical analysis and numerical simulations, we demonstrated that in many scenarios gene transfer could promote the emergence of alternative stable states which otherwise would not exist. We further extended our analysis to metacommunities consisting of multiple local habitats. Our results suggested that while multistability enabled by HGT didn’t affect local diversity, it enabled the emergence of species diversity at regional scale. We also discussed the influence of different confounding factors on the applicability of the predictions. Our results provide valuable insights for microbial ecology, health and microbiome engineering.

Results

The role of HGT on bistability of two-species populations

To illustrate the basic concepts, let’s first consider a population of two competing species. The mutual inhibition between the two members creates a positive feedback, making bistability possible1,14. Without HGT, the community dynamics can be simulated by the classic Lotka-Volterra (LV) model that accounts for species growth rates (μ1 and μ2), interspecies competition (γ1 and γ2) and the dilution rate (D). Given the simplicity of this model, the condition of the system being bistable can be analytically solved as γ1 > ϕ2/ϕ1 and γ2 > ϕ1/ϕ2 where and (see Methods for more details). When this condition is satisfied, either species can dominate the population at steady state, depending on their initial abundances. Otherwise, the population will always rest at the same steady state regardless of the initial composition (Fig. 1a).

Horizontal gene transfer (HGT) promotes bistability of communities consisting of two competing microbes.

(a) Schematics of monostable and bistable systems. A monostable system will always rest at the same steady state no matter where the system starts from. In contrast, a bistable system can reach two distinct steady states, depending on the initial ratio between the two competitors.

(b) The phase diagram of the population with (bottom panel, η = 0.2 hr−1) or without HGT (top panel). Different combinations of μ1 and μ2 between 0 and 1 hr−1 were tested. For each combination, the initial abundance of each species was randomized 200 times between 0 and 1 following uniform distribution. The system was monostable if all initializations led to the same steady state. Otherwise, the system was bistable. All combinations of μ1 and μ2 that led to bistability were marked in blue in the diagram. Other parameters are γ = 1.1, k = 0.005 hr−1, D = 0.2 hr−1, . Two states, P (μ1 = 0.3 hr−1, μ2 = 0.8 hr−1) and Q (μ1 = 0.4 hr−1, μ2 = 0.6 hr−1), were labeled as examples. (c-d) Population dynamics of P and Q under 100 different initializations. Here, the dynamic changes of species 1’s abundance were shown. Three different HGT rates (0, 0.1, 0.2 hr−1 from top to bottom) were tested. When gene transfer rate increases, the system at Q changed from monostable to bistable, while the system at P remained monostable.

To analyze the effect of HGT on bistability, we applied a mathematical model that we previously developed to describe gene transfer between competing species27. Compared with the classic LV model, the new framework considered the growth rate of each species as a combination of two components contributed by non-mobilizable genes and mobile genes, respectively (see Methods for more details). Due to the fitness effects of mobile genes (denoted as λ1 and λ2, respectively), their gain and loss within each species caused the dynamic change of species growth rates. Gene transfer (with the rate denoted as η) allows the two competitors to partially exchange their fitness benefits or disadvantages, and we sought to understand whether this process would reshape the stability landscape of the system.

Due to the increase of modeling complexity, analytical solution becomes infeasible for the new framework. Instead, we determined whether the system was bistable by numerical simulations with randomized initial compositions26. In general, the resolution of this approach increases with the number of randomizations. Without loss of generality, here we randomly initialized the species abundances for 200 times between 0 and 1 following uniform distributions (Fig. S1a). Then we simulated the population dynamics for each initialization until steady state and obtained their species abundances. A system is called monostable if all populations starting from different compositions converge into the same state. In contrast, if two attractors exist, the system is called bistable (Fig. 1a). We calculated the number of stable states for a wide range of μ1 and μ2 values and drew the phase diagram (Fig. 1b). Whether it is easy for the system to achieve bistability is reflected by the total area of bistability region on the phase diagram. Our numerical results suggested that increasing HGT enlarged the region, creating bistability in a lot of systems which would otherwise be monostable (Fig. 1b-d). Here, the number of initializations was not critical for the conclusion (Fig. S1a). The results also hold for populations with a variety of interspecies competitions or gene transfer rates (Fig. S1b-d), suggesting the robustness of our conclusion against a variety of confounding factors.

Our results suggested that increasing HGT made it more feasible for the system to become bistable. However, this should not be taken as an assertion that the system will always be bistable when HGT exists (Fig. 1c and d). The effect of HGT in creating bistability in a specific community depends on the competition strength, gene transfer rate and the fitness difference between species (Fig. 1b-d). Bistability will not appear when γ1γ2 < 1, regardless of HGT rate (Methods). In addition, for two species with large fitness difference, the system might remain monostable when gene transfer rate increases (Fig. 1c). In this case, the strong competitor will always win despite the gene flow. The effect of HGT in promoting bistability is more applicable between species with smaller fitness difference and stronger competition (Fig. 1b-d, Fig. S1b).

HGT shapes the stability landscape of multi-species communities

The modeling framework can be readily generalized to more complex communities consisting of multiple competing species by taking into account the gene transfer between every pair of species (see Methods for more details). To understand how gene transfer rate affects the stability landscape of complex microbiota, we estimated the number of alternative stable states by randomly initializing the abundance of each species 500 times between 0 and 1 following uniform distributions. For each initialization, we simulated the population dynamics for 2000 hours till the systems reached equilibrium and calculated the steady-state abundances of different species (Fig. S2). Limit cycles were not observed in our simulations (Fig. S3).

The stability landscape of a complex community can be characterized by the absolute number of stable states. To this end, we grouped the steady states into different attractors using a given threshold. Two steady states with their Euclidean distance smaller than the threshold belonged to the same group (Fig. S4). However, the estimation of this number can be biased by the intrinsic limitations of the numerical approach. For instance, the limited number of initializations can underestimate this number. In contrast, the inability of a system to reach equilibrium within the limited time window of simulations can cause its overestimation. To overcome such limitations, we developed an alternative metric χ to measure the multistability of a system by borrowing the concept of ‘entropy’. Let φ be the number of stable states and xi be the relative size of each domain of attraction (Σi xi = 1, i = 1,2, …, φ). χ is calculated as χ = exp (− Σi xilogxi), which accounts for the absolute number of stable states as well as the size of each basin. χ equals 1 in monostable systems. χ gets larger when more stable states emerge or the distribution of the basin sizes become more even. Numerical tests suggested that χ was robust to the variations of distance threshold, simulation timespan or number of initializations (Fig. S5).

As shown in Fig. 2, when HGT is absent, the number of stable states is usually small in a multi-species population. Increasing gene transfer rate reshaped the stability landscape of a microbiome by creating additional stable states. Each stable state is characterized by a different dominant species (Fig. 2a). This prediction is generally applicable to communities consisting of different number of species (Fig. S6).

The number of alternative stable states in multi-species communities increases with HGT rate.

(a), A schematic of the stability landscape of a population composed of multiple competing species with or without HGT. Each alternative stable state is characterized by a different dominant species (shown in filled circles).

(b), The emergence of alternative stable states in a five-species community with gene transfer. We simulated the dynamics of 500 parallel communities, each with the same kinetic parameters but different initial compositions. The species compositions of these communities at three different timepoints were shown by principal component analysis from left to right. Without HGT (η = 0), all populations converged into a single stable state. In contrast, with HGT (η = 0.2 hr−1), four different dynamic attractors coexist. Here, populations reaching different stable states were marked with different colors. Other parameters are , μ1 = 0.52 hr−1, μ2 = 0.48 hr−1, μz = 0.44 hr−1, μ4 = 0.60 hr−1, μ5 = 0.31 hr−1, γ = 1.05, k = 0.005 hr−1, D = 0.2 hr−1.

(c), The number of stable states (top panel) and the multistability coefficient χ increase with HGT rate in communities consisting of 8 competing species. We calculated the number of stable states by randomly initializing the species abundances 500 times between 0 and 1 following uniform distributions. Then we simulated the steady states and clustered them into different attractors by applying a threshold of 0.05 on their Euclidean distances. The data were presented as mean ± standard deviation of 10 replicates. Each replicate corresponded to a different combination of randomized species growth rates.

So far, our analysis has been focused on communities where the interaction network is fully connected, and each species competes with every other member. In each of such populations, it is highly infeasible to achieve the stable coexistence of different species. While this analysis provides the basic proof of principle, natural environments often consist of many niches28. Species living in the same niche compete with each other while species from different niches can coexist. To understand whether our conclusion is still applicable when niche effects are considered, we extended our analysis to communities carrying a random number of niches. Each species was randomly allocated into one of the niches, and the stability landscape was analyzed in the same way by randomly initializing the systems (see Methods for more details). Our results suggest that increasing HGT rate still promotes multistability, regardless of the niche structure of the competing communities (Fig. 3).

HGT promotes the multistability of communities with multiple niches.

(a), A schematic of a microbial population consisting of multiple niches. Species within each niche share the same carrying capacity and compete with each other, while species from different niches can coexist.

(b-c), Multiple alternative stable states emerge in niche communities with gene exchange. Here, a system consisting of 10 species and 4 niches were shown as examples. The steady-state compositions of 100 parallel populations, each with different initializations, were calculated with (η = 0.1 hr−1) or without HGT. The heights of the colored bars represent the relative abundances of different species. Without HGT, 2 stable states existed for the 100 populations, while with HGT, 6 different stable states existed. The carrying capacity of each niche was randomized between 0 and 1. Other parameters are , k = 0.005 hr−1, D = 0.2 hr−1. Interspecies competition strength equals 1.05 within each niche.

(d), The multistability χ in niche communities increases with HGT rate. Communities of 8 species were considered as an example. The data were presented as mean ± standard deviation of 10 replicates. For each replicate, the species growth rates were randomized between 0.3 and 0.7 hr−1 following uniform distributions. The number of niches were randomized as 2, 3 or 4.

HGT-enabled multistability promotes the regional coexistence of competing species

Our results highlight the role of HGT on shaping the stability landscape of closed local communities. Natural microbiomes, however, exist at scales broader than a localized population. A set of local communities subject to similar environments, also known as ‘patches’, form a metacommunity29,30. When multiple stable states coexist, the compositions of different patches may differ from each other due to the variability of their initial species occupation. Collectively, the multistability in the local patches might give rise to the regional species diversity in the metacommunity. This intuition leads to our hypothesis that HGT allows the stable coexistence of multiple competing species at regional scale, even when species exclusively outcompete each other in every local population.

To examine this hypothesis, we constructed metacommunities of 100 patches connected by dispersal of microbes. Without loss of generality, we considered a 2-D array of local patches and assumed that the dispersal only occurred between adjacent patches (Fig. 4a, see Methods for more details). Each local patch carried 10 competing species that exchanged mobile genes with each other. The dynamics of each patch was governed by the same set of parameters but initialized with different species compositions. We then numerically simulated the dynamics of the metacommunity till it reached steady state.

HGT promotes the regional diversity in metacommunities of local populations connected by dispersal.

(a), A schematic of microbial dispersal within the metacommunity. The filled circles represent different local patches.

(b), A metacommunity composed of 100 connected patches was shown as an example. Each local population contained 10 competing species and was randomly initialized with different species compositions. Dispersal occurred between neighboring patches with a rate of 0.001 hr−1. Without HGT, the systems were monostable, resulting in low regional diversity. HGT promotes the regional diversity by the emergence of multistability in each local patch. The collective diversity of the metacommunity was calculated by pooling all the local populations together. Other parameters are , γ = 1.05, k = 0.005 hr−1, D = 0.2 hr−1.

(c), The regional diversity of a metacommunity increases with HGT rate. We first pooled all local populations together and calculated the total abundance of every species in the metacommunity. Then we quantified the regional diversity using Shannon index. The data were presented as mean ± standard deviation of 10 replicates.

Without HGT, the metacommunity exhibited strong spatial homogeneity at regional scale (Fig. 4b). Different patches rested at similar species compositions. In contrast, HGT created a spatial pattern of species distributions across different patches, due to the emergent multistability (Fig. 4b). To analyze the influence of HGT on regional diversity, we pooled all local populations together and calculated the collective species diversity in the pool using Shannon index. Our results suggested that while local diversity remained low regardless of HGT rate, the regional diversity of microbial species increased with HGT rate (Fig. 4b and c). The conclusion is applicable for different dispersal rates (Figs. S7 and S8). These results highlight the fundamental role of HGT in community stability and diversity at local and regional scale.

The influences of different confounding factors on the role of HGT in mediating population multistability

Our previous analysis has focused on the role HGT in communities where: (1) MGEs only affect species growth rates and have no influences on inter-species competition; (2) the fitness effects of MGEs (λ1 and λ2) are independent of the host species. To understand how these assumptions affected the prediction and clarify the potential caveats of our conclusion, we further extended our modeling framework by accounting for different confounding factors.

In nature, the horizontal transfer of some genes can change the competition strengths between different species. For instance, the sharing of many mobile genes can promote niche overlapping, leading to an increase of competition strength31,32. Some mobile genes encode public goods that benefit the host’s neighboring species33,34. HGT of these genes will reduce the strength of interspecies competition. To understand how the transfer of these genes would influence community bistability, we adapted the previous model by considering the dynamic change of competition strength during HGT (see Supplementary Materials for more details). In particular, for two-species populations, we divided the interspecies competitions into two components: the basal competition strengths (γ1 and γ2) and the added parts by HGT (δ1 and δ2). The overall competition strength was calculated as and , respectively. Positive δ ‘s mean HGT promotes the inter-species competition, while negative δ ‘s mean HGT reduces competition. Numerical simulations with randomized parameters suggested that when δ’s were positive, HGT still promoted the bistability, whereas if δ’s were negative, gene transfer would reduce the chance of bistability (Fig. 5a-d). We further extended our analysis to multispecies communities and calculated the number of stable states. As shown in Fig. 5e and f, the number of stable states increases with HGT rate when MGEs promotes competition, while decreases with HGT rate when MGEs reduces competition.

The effects of HGT on population multistability when MGEs promote or reduce competition.

(a-b), The schematics of HGT promoting or reducing competition.

(c), For populations of two species, when MGEs promote competition, increasing HGT rate enlarged the area of bistability region in the phase diagram. Here, δ describes the effect of mobile genes on the competition strength. Positive δ represents HGT promoting competition. In numerical simulations, we tested three different δ values (marked in different colors). When calculating the area of bistability region, we randomized μ1 and μ2 500 times between 0 and 1 hr-1 following uniform distributions while keeping and constants. For each pair of growth rates, we randomized the initial abundance of each species 200 times between 0 and 1 following uniform distributions. The system was monostable if all initializations led to the same steady state. Otherwise, the system was bistable. Then we calculated the fraction of growth rate combinations that generated bistability out of the 500 random combinations. Other parameters are γ1 = γ2 = 1.1, k = 0.005 h−1, D = 0.2 h−1.

(d), When MGEs reduce competition, the area of bistability region decreases with HGT rate. Three negative δ values were tested and shown as examples here.

(e), For populations of 5 species, when MGEs promote competition, the number of stable states increases with HGT rate. We calculated the number of stable states by randomly initializing the species abundances 500 times. Then we simulated the steady states and clustered them into different attractors by applying a threshold of 0.05 on their Euclidean distances. The data were presented as mean ± standard deviation of 10 replicates. Each replicate corresponds to a different combination of randomized species growth rates. Other parameters are γi,j = 1.1, k = 0.005 h−1, D = 0.2 h−1, δ = 0.5.

(f), When MGEs reduce competition, the number of stable states decrease with HGT rate. The data were presented as mean ± standard deviation of 10 replicates. δ = −0.5 was used in the simulation.

In previous analysis, we also assumed that the effects of MGEs on host growth (λ1 and λ2) were independent of the host species. In nature, however, the same MGE can have different fitness effects in different genomic backgrounds, a phenomenon called epistasis35-38. There are two types of epistasis: magnitude epistasis, where the host genomic background only affects the magnitude but not the sign of λ value, and sign epistasis, where the same MGE is burdensome in one species while beneficial in the other host. The model can be readily extended to account for these two types of epistasis (see Supplementary Materials for more details). In particular, we dissected the growth rates μ1 and μ2 into two components: the basal growth rate ( and ), and the fitness effects (λ11 and λ22) of the mobilizable genes. The growth rates of p1 and p2 were calculated as μ1(1 + λ12) and μ2(1 + λ21), respectively. Here λij (i, j = 1,2) stands for the fitness effects of the j-th mobile gene in the i-th species. Without epistasis, the fitness effects are independent of the host species (λ11 = λ21 and λ22 = λ12). The epistasis can be quantified by two ratios: ξ1 = λ21/λ11 and ξ2 = λ12/λ22. ξ1 > 0 and ξ2 > 0 represent magnitude epistasis, whereas ξ1 < 0 or ξ2 < 0 represent the sign epistasis. We carried out additional numerical simulations to evaluate how different epistasis types influence the role of HGT in community stability. Our results suggested that with magnitude epistasis, HGT still enlarged the area of bistability region (Fig. 6a and b). In contrast, sign epistasis overturned HGT’s role on community stability: with sign epistasis, increasing gene transfer rate reduced the area of bistability region (Fig. 6c and d). These results suggest that MGE epistasis might add another layer of complexities into the interplay between HGT and community stability.

The influence of epistasis on the role of HGT in mediating population bistability.

(a-b) With magnitude epistasis, increasing HGT rate still promotes bistability in populations of two competing species.

(c-d) With sign epistasis, increasing HGT rate reduces the area of bistability region. For each type of epistasis, we considered two scenarios. In a and c, host genetic background influences the fitness effect of only one MGE, while in b and d, host genetic background influences both MGEs. ξ1 and ξ2 are defined as λ21/λ11 and λ12/λ22, respectively. ξ1 > 0 and ξ2 > 0 represent magnitude epistasis, while ξ1 < 0 or ξ2 < 0 represents the sign epistasis. Other parameters are γ1 = γ2 = 1.1, k = 0.005 h−1, D = 0.2 h−1.

The fitness effect of an MGE can be extreme in some cases. For instance, under strong antibiotic selection, only cells carrying the antibiotic resistant genes (ARGs) can survive. To examine whether our conclusion is still applicable in this scenario, we generalized our model by considering the transfer of an ARG in a population under strong antibiotic selection (see Supplementary Materials for more details). Specifically, we assumed that the ARG was initially carried by one species in the populations. Without HGT, only the donor species carrying the ARG can survive due to antibiotic selection (Fig. 7a). In this case, the system will always be monostable regardless of the initial species compositions. With HGT, the ARGs can be transferred to other species, creating chances for other species to survive (Fig. 7b). To understand how HGT affects community stability under this condition, we first analyzed the phase diagram of the two-species population. Our simulation results suggested that bistability became possible when there existed gene transfer and increasing HGT rate in general enlarged the bistability area (Fig. 7c and d). Here, without loss of generality we assumed that the ARG was initially carried only by species 1. Bistability was more feasible when μ2 > μ1 (Fig. 7c), suggesting that the effect of HGT on promoting bistability will be stronger when the weak competitor acts as the ARG donor. We further analyzed the number of stable states in complex communities transferring an ARG. As shown in Fig. 7e and f, for these populations under antibiotic selection, increasing transfer rate of the ARG also promoted the emergence of multistability, suggesting the general applicability of our prediction in this scenario.

The horizontal transfer of ARGs creates chances of multistability for communities under strong antibiotic selection.

(a-b), The schematics of communities under strong antibiotic selection. Without HGT, only the donor species carrying the ARG can survive. HGT allows the other species to acquire ARG from the donor, creating opportunities for the other species to survive under antibiotic treatments.

(c), The phase diagram of two-species populations transferring an ARG. Different combinations of μ1 and μ2 between 0 and 1 hr−1 were tested. All combinations of μ1 and μ2 that led to bistability were marked in blue in the diagram. Other parameters are γ = 1.1, k = 0.005 hr−1, D = 0.2 hr−1, . Here, we assumed that the ARG was initially carried only by species 1. Bistability was more feasible when μ2 > μ1.

(d), Increasing HGT rate enlarges the area of bistability region in the phase diagram of two competing species transferring an ARG. Three different competition strengths were tested (marked in different colors).

(e-f), For communities of 5 or 8 species under antibiotic selection, increasing HGT rate of an ARG promotes the emergence of alternative stable states. We calculated the number of stable states by randomly initializing the species abundances 500 times. Then we simulated the steady states and clustered them into different attractors. The data were presented as mean ± standard deviation of 10 replicates. Each replicate corresponds to a different combination of randomized species growth rates. Other parameters are γi,j = 1.5, k = 0.005 h−1, D = 0.2 h−1.

Discussion

Our work predicts that in many scenarios horizontal gene flow promotes multistability in communities of competing microbes. At the scale of local populations, the emergence of alternative stable states undermines the global stability of the community, making the system less resilient to perturbations on species abundances and creating greater opportunities for regime shifts. In metacommunities, gene flow creates spatial heterogeneity of species compositions in different patches. Despite the low diversity in every local patch, HGT-driven multistability creates substantial diversity at regional scale. In microbial ecology, how microbes competing for a limited number of resources stably coexist has been a long-standing puzzle39,40. Our results highlight the importance of considering gene exchange when addressing this fundamental challenge.

While this work predicts the potential role of HGT in microbiota diversity and stability, several caveats need to be considered when applying this prediction. For MGEs with sign epistasis, or MGEs that relax interspecies competitions, HGT might reduce the number of stable states. The relative importance of gene flow in a real community might also be context dependent. The effect of HGT in promoting multistability might be stronger between species with intense competition and similar fitness. In environments with high cell density and abundant MGEs, such as mammal gut, the role of HGT will also be more prominent41. In contrast, the influence of HGT on community multistability might become less important in populations with poor cell density and slow gene transfer. Nevertheless, by comparison with empirical estimates of plasmid conjugation rates from a previous study42, the HGT rates in our analysis are biologically relevant in a variety of natural environments (Methods).

Previous studies documented the presence of multiple alternative community types in soil15, human vaginal, gut and oral microbiomes2,43,44. In particular, multiple groups of bacteria in human gut exhibit robust bistable abundance distributions, resulting in over 32 different stable state combinations45. Our results provide a plausible explanation for the large number of community types in these populations, even though our interpretation doesn’t exclude other mechanisms including gradients of environmental factors1. Disentangling the contributions of HGT from the other factors in diverse natural microbiomes is a promising topic for future research.

Controlling the multistability of complex microbiota has important applications in human health. For instance, in gut microbiome, the development of unhealthy stable states is closely related to many diseases such as obesity, type 2 diabetes mellitus and gastrointestinal disorders46-48. Strategies have been developed to overcome the resilience of unhealthy states, including fecal transplantation or supplementation with probiotics18. Our work predicts that the horizontal transfer of MGEs such as plasmids might be an additional mechanism for the emergence of the unhealthy states. Indeed, human gut is a ‘hotspot’ of gene flow41. Therefore, inhibiting the MGE spread using small molecules might offer new opportunities to reshape the stability landscape and narrow down the attraction domains of the disease states49,50.

Methods

Modeling framework of two competing species sharing mobile genes

We used a mathematical model developed in a previous work to analyze the population dynamics of two competing species transferring mobile genes with fitness effects27. Briefly, the model includes four ordinary differential equations (ODEs):

S1 and S2 are the abundances of the two species. γ1 and γ2 describe the strengths of interspecies competition. Nm is the maximum carrying capacity, and D is the dilution rate. Pi (i = 1, 2) represents the abundance of the subpopulation in the i-th species that acquires the mobilizable genes from its competitor. μi represents the maximum growth rate of the i-th species when interspecies HGT is absent. μi is calculated as where is the basal growth rate determined by the non-mobilizable genes and λi is the fitness effect of the mobile genes. Interspecies gene flow allows each species to acquire mobile genes from its competitor, resulting in the dynamic change of the effective species growth rates , which is calculated as: and . and stand for HGT rates. Here, we assumed that intraspecies and interspecies transfer occurred at the same rates. k1 and k2 are the loss rates of the mobile genes.

Let s1 = S1/Nm, s2 = S2/Nm, p1 = P1/Nm and p2 = P2/Nm. The model can be further simplified as

Here, s1, s2, p1 and p2 are dimensionless. ηi is related with by . Therefore, the unit of η1 and η2 is hour−1. Without HGT, p1 and p2 become consistently 0, and the model is equivalent to the classic Lotka-Volterra (LV) framework:

Stability analysis of the classic LV model

The classic LV model has four fixed points:. The Jacobian matrix of each fixed point can be calculated as, where and are the species abundances at each fixed point.

The system is bistable only when ℊ2 and ℊ3 are both stable. The Jacobian matrix of ℊ2 is . The eigenvalues of J2 can be obtained as and D – μ2. ℊ2 is stable when both eigenvalues are negative, which leads to γ2 > ϕ1/ϕ2 where and . Similarly, ℊ3 is stable when γ1 > ϕ2/ϕ1.

Therefore, the condition for the system being bistable is γ2 > ϕ1/ϕ2 and γ1 > ϕ2/ϕ1.

Bistability analysis of two-species populations by numerical simulations

We performed numerical simulations to determine the bistability of a population consisting of two competing species transferring mobile genes. Specifically, for a system with given parameters, we randomly initialized the species abundances between 0 and 1 following uniform distributions for 200 times. Then we simulated the dynamics of each population till the abundance of every species reached equilibrium, and we treated it as a steady state. We grouped the 200 steady states into different attractors by applying the threshold of 0.01. Two steady states with their Euclidean distance smaller than the threshold belonged to the same group. The Euclidean distance was calculated as where k is the species index and si,k represents the abundance of the k-th species in the i-th steady state.

Theoretical framework of multiple competing species transferring mobile genes

For a community of m competing species that share mobile genes with each other, the modeling framework consists of two groups of ODEs:

si is the abundance of the i-th species. pij is the size of the subpopulation in the i-th species that acquires sj-originated mobile genes. γij is the negative interaction that sj imposes on the i-th species. D is the dilution rate. , the effective growth rate of si, is calculated by is defined as where is the basal growth rate, λij is the fitness effect of the sj-originated mobile genes in the i-th species. ηjki is the transfer rate of the sj-originated genes from species k to species i. kij is the gene loss rate.

Analysis of multistability of complex communities

For a population of m species with given niches and kinetic parameters, we initialized the system 500 times in the ranges of 0 < si < 1 following uniform distributions. The initial conditions for pij’s are pii = si and pij = 0 (i ≠ j). Then the dynamics of each population was simulated for 2000 hours till the system reached the steady state. For each steady state we obtained the corresponding species composition and grouped the 500 solutions into different stable states based on their Euclidean distances.

Niche-based model of multi-species communities

For a community consisting of m competing species, let the number of niches be l. Each niche has a maximum carrying capacity . Each species lives in one of the niches. Let M be the map from species index to the index of its niche. Mi = k means that the i- th species lives in the k-th niche. For species i and j, if Mi = Mj, these two species belong to the same niche and compete with each other. Species from different niche do not directly interact. Therefore, the model becomes:

In simulations, the number of niches was randomized between 2 and m/2. The maximum carrying capacity of each niche was randomized between 0 and 1 following uniform distributions. Then each species was randomly allocated into one of the niches.

Theoretical framework of metacommunity dynamics

Consider a metacommunity composed of u × v local patches. The location of each patch can be described by two coordinates [a, b] (a = 1,2, …, u and b = 1,2, …, v). Let ω be the dispersal rate between the adjacent patches. The dynamics of the metacommunity can then be described by generalizing equations (11) and (12) to each local patch:

For the patch located at [a, b], sa,b,i is the abundance of the i-th species, and pa,b,i,j is the size of the subpopulation in the i-th species that acquires sj-originated mobile genes. The effective growth rate is calculated as.

and are the gains of sa,b,i and pa,b,i,j by dispersal, respectively. We assumed that microbial dispersal only occurred between adjacent patches. Therefore, was calculated by summing the dispersal terms at four directions: left, right, up, and down. The left component can be calculated as

Similarly, and is calculated as:

can then be obtained as can be formulated in the similar way.

The collective species diversity in a metacommunity

To analyze the total species diversity in a metacommunity, we first calculated the total abundance of each species across all patches: si = Σa=1,2,…,u Σb=1,2,…,v sa,b,i. Then we measured the diversity using Shannon index, a common metric used in ecology. Shannon index can be calculated as. sT is the total species abundance: .

Biological relevance of our analysis

The empirically measured gene transfer rates (denoted as ηc, with the unit of cells−1 · mL · hour−1) should be multiplied by Nm, the maximum carrying capacity of the population, before being used our model. Therefore, the transfer rate η in our analysis are several orders of magnitude higher than the experimentally measured rates51,52. In this work, we focused on HGT rate in the range of 0 < ηi < 0.4 hour−1. Whether this range is biologically relevant depends on Nm of the population. Considering the typical estimates of plasmid transfer rates cross species (10−13 cells−1 · mL · hour−1, from Klebsiella pneumoniae to Escherichia coli) and within species (10−7 cells−1 · mL · hour−1, between E. coli strains) provided by a previous study51, HGT rates in our analysis are relevant in various environments from human colon (∼1012 cells per gram) 53 to deep oceanic subsurface (106∼108 cells per gram)54.

Code Availability

All the codes and source data associated with the numerical simulations and analysis of this manuscript are available at the Github repository: (https://github.com/twang1993/MultistabilityByHGT).

Acknowledgements

This work is supported by Shenzhen Institute of Synthetic Biology Scientific Research Program (Grant No. HSE199011086). The funder had no role in study design, data collection and analysis, the decision to publish, or the preparation of the manuscript.

Declaration of interests

The authors declare no competing interests.