Introduction

Microbial communities living in similar environments often display different community types, each characterized by a unique composition of species^1-4. Changes in environmental factors, fluctuations of species abundances, or transient invasions can trigger abrupt and irreversible transitions between two types^5-7. Such transitions, also known as regime shifts, lead to the dramatic alterations of microbiome functions and are tightly associated with host diseases, soil fertility, marine biogeochemical flux, or bioproduction yield^8-12. Understanding the origin of different community types has broad implications in environments, engineering and health^1,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 microbiomes^1,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 interactions¹⁴.

Understanding the origin and the underlying determinants of multistability in complex microbiota is an important challenge in microbial ecology and evolution. In host-associated microbiota, regime shifts between alternative stable states can trigger the critical transitions between healthy microbiome and dysbiosis^8,18. Unravelling how different biological processes drive community multistability is important for the predictive control of complex microbiota, and can inspire novel treatments to modulate the communities towards the desired stable states or away from the unhealthy states.

Mathematical models are indispensable to decipher the emergent properties of complex populations¹⁹. Different ecological models have been developed to analyze the origin of multistability in microbial communities. For instance, based on the generalized Lotka-Volterra (LV) model, Gonze D et al. showed that multistability could emerge in communities where different species mutually inhibited each other¹. Dubinkina V et al. applied a consumer-resource model in which microbes competed for different metabolites, and showed that multistability required microbial species to have different stoichiometries of essential nutrients¹⁴. These works provide critical insights on the determinants of community multistability, yet the roles of many other mechanisms remain to be understood. Chief among these is horizontal gene transfer (HGT), a process where microbes share mobile genetic elements (MGEs) with their neighbors²⁰.

Indeed, HGT is prevalent in microbial world and mobile genes occupy a substantial proportion in microbial gene pool^21-23. For instance, among the closely related microbial isolates in the same human body site, HGT occurs in more than 40% of the species pairs²⁴. In the draft genomes of human gut microbes, tens of thousands of genes were mobilizable²². Depending on the encoded traits, MGEs shape microbial growth rates, host adaptability and even social interactions between microbes²⁵. The dynamic gain and loss of MGEs are also a major source of evolutionary innovations of prokaryotic genomes²⁶. 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 multistability^1,14,27. By combining mathematical modelling and numerical simulations, we demonstrated that the horizontal transfer of plasmids, a major type of MGEs, could promote the emergence of alternative stable states which otherwise would not exist. We further extended our analysis to scenarios where HGT changed interspecies interactions, where microbial communities were subjected to strong environmental selections and where microbes lived in metacommunities consisting of multiple local habitats. We also analyzed the role of different mechanisms, including interspecies interaction strength, the growth rate effects of MGEs, MGE epistasis and microbial death rates in shaping the multistability of microbial communities. These results created a comprehensive framework to understand how different dynamic processes, including but not limited to HGT rates, collectively shaped community multistability and diversity. 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. Without HGT, the community dynamics can be captured by the classic LV model that accounts for species growth rates (µ₁ and µ₂), interspecies interaction strengths (γ₁ and γ₂) and the dilution or death rate (D). When all the parameters are given, one can simulate the system dynamics from the initial abundances of the two species. The system will reach an equilibrium, characterized by the steady-state abundances of the two species. Bistability means the coexistence of two distinct stable states (Fig. 1a)^1,14. Depending on the initial abundances, either one of the two species can dominate the community. With bistability, even when the model parameters are not altered, a change in initial species abundances is sufficient to drive the system to a different stable state (Fig. 1a). In contrast, if monostable, the population will always rest at the same steady state regardless of the initial composition (Fig. 1a). Given the simplicity of this model, the condition of the system being bistable can be analytically solved as γ₁ >ϕ₂/ϕ₁ and γ₂ >ϕ₁/ϕ₂ where ϕ and (see Methods for more details). This condition predicted the thresholds on the strengths of interspecies interactions in order for the system to be bistable.

Fig. 1.

To analyze the effect of HGT on bistability, we applied a mathematical model that we previously developed to describe the conjugative transfer of plasmids between competing species²⁸. While in classic LV model µ₁ and µ₂ are constants, our model dissected species growth rates into two components: the basal growth rates ( and ) determined by chromosomal genes, and the growth rate effects (λ₁ and λ₂) of plasmids. The two components were assumed to combine multiplicatively: . HGT creates subpopulations (denoted as p₁ and p₂) within each species that acquire the mobilizable genes from its competitor. The dynamics of p₁ and p₂ were governed by cell growth, plasmid transfer (the rate denoted as η) and plasmid loss (the rate denoted as κ). The gain and loss of p₁ and p₂ in turn lead to the temporal change of the effective growth rates ( and ) of each species: , (see Methods for more details). Therefore, plasmid transfer allows the two competitors to partially exchange their growth benefits or disadvantages. Here, we assumed that the species internal inhibition was independent of the acquired MGEs. We sought to understand whether the process of HGT would reshape the stability landscape of the system.

After HGT being introduced, analytical solution becomes infeasible due to the increase of modeling complexity. Instead, we determined whether the system was bistable by numerical simulations with randomized initial compositions²⁷. 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 µ₁ and µ₂ 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 many systems which would otherwise be monostable (Fig. 1b-d). Here, the number of initializations was not critical for the conclusion (Fig. S1a). The results were not fundamentally changed when the interspecies interaction strengths or gene transfer rates were altered (Fig. S1b-d), suggesting the robustness of our conclusion.

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). For two species with large growth rate 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 growth rate difference (Fig. 1b and d).

Previous studies showed that the cell death rate determines the bistability of two interacting species²⁹. To evaluate the impact of death rate D on the interplay between HGT and system bistability, we performed additional analysis by calculating the bistability probability under different D values (Fig. S2). Our results suggested that varying death rate indeed changed the bistability of the system. When the death rate equaled zero, mobile genetic elements that only modified growth rates would have no effects on population bistability. These results highlighted the importance of added death rate in driving multistability.

The strength of interspecies interaction also shaped the community bistability. By calculating the bistability probability under different values of γ (Fig. S3), we showed that while changing γ didn’t fundamentally alter our conclusion, the influence of HGT on population bistability became weak when interspecies interaction strength got smaller than 1. These results suggested that the role of HGT was more significant in population with strong competitions.

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. Consider a community of m species. Let s_i be the abundance of the i-th species. γ_ij represents the interspecies interaction strength that the i-th species imposes on the j-th species (i, j = 1,2, …, m). η_jki is the transfer rate of the s_j -originated plasmid from species k to species i. κ_ij is the loss rate of the j-th mobile genetic element in the i-th species. The growth rate effect of the j-th plasmid on the i-th species is denoted as λ_ij. p_ij describes the abundance of the subpopulation in the i-th species that carries the plasmid from the j-th species. The effective growth rate of the i-th species can be calculated based on all the plasmids that it carries: . The population dynamics can then be simulated by m + m² equations that describe the temporal change of s_i and p_ij, respectively (see Methods for more details).

To understand how changing 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. S4). Limit cycles were not observed in our simulations (Fig. S5).

The calculate the absolute number of stable states, we grouped the steady states into different attractors using a given threshold. Steady states with their Euclidean distances smaller than the threshold belonged to the same group (Fig. S6). 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 x_i be the relative size of each domain of attraction (∑_i x_i = 1, i = 1,2, …, φ). χ is calculated as χ = exp (− ∑_i x_i logx_i), 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. S7).

As shown in Fig. 2, when HGT is absent, the number of stable states was usually small in a multi-species population. Increasing gene transfer rate reshaped the stability landscape of a microbiome by creating additional stable states (Fig. 2a and b). Each stable state is characterized by a different dominant species (Fig. 2a). The number of stable states is positively associated with HGT rate (Fig. 2c, Fig. S8 and S9). This prediction is generally applicable to communities consisting of different number of species (Fig. S8 and S9). For communities of two species, bistability is not possible when γ_ij’s are smaller than 1. However, for communities of multiple species, multistability becomes possible even if all γ_ij values are smaller than 1 (Fig. S10). The effects of HGT on multistability held for different strengths of interspecies interactions (Fig. S10). Initializing the p_ij subpopulations with non-zero abundances didn’t change the conclusion, either (Fig. S11).

Fig. 2.

Here we assumed the interactions among different species were homogeneous. To evaluate whether this assumption was critical, we extended our analysis by considering heterogenous interaction strengths in multispecies communities (Fig. S12). In particular, we randomly sampled γ_ij values from uniform distributions. Our results suggested the mean value and variance of γ_ij played a role in shaping multistability. The effects of HGT on community multistability becomes stronger when the mean value of γ_ij gets larger than 1 and the variance of γ_ij is small (Fig. S12). Nevertheless, the heterogeneous distribution of γ_ij didn’t fundamentally change our conclusion.

MGEs can drive large growth rate differences when they encode adaptative traits like antibiotic resistance. In many cases, however, the growth rate effects of MGEs can be small. To evaluate the influence of λ_ij magnitude, we also analyzed different ranges of growth rates effects of mobile genetic elements, by sampling λ_ij values from uniform distributions with given widths (Fig. S13). Greater width led to larger magnitude of growth rate effects. We used five-species populations as an example and tested different ranges. Our results suggested that multistability was more feasible when the growth rate effects of MGEs were small (Fig. S13b). The qualitative relationship between HGT and community was not dependent on the range of growth rate effects (Fig. S13a).

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 niches³⁰. Species living in the same niche competes while species from different niches can coexist (Fig. S14a). 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 suggested that increasing HGT rate still promoted multistability, regardless of the niche structure of the competing communities (Fig. S14b-d).

The role of HGT on community multistability when mobile genes modify interspecies interactions

Our previous analysis has focused on the role of HGT in communities where MGEs only affect species growth rates and have no influences on interspecies interactions. To evaluate how this modeling structure and the underlying assumptions affected the prediction, we extended our framework by accounting for the modifications of interspecies interactions by mobile genes.

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

Fig. 3.

These results suggested that the effects of HGT on bistability depended on how mobile genes shaped species growth rates and competitions. To place these two factors in a comprehensive framework, we further generalized the model structure, by accounting for the scenario where mobile genes modified growth rates and competitions at the same time. The effect of mobile genes on growth rates was represented by the magnitude of λ’s, and the influence on competition is described by the parameter δ (Fig. 3g to i). By varying these two parameters, we can evaluate how the model structure and the underlying assumptions affected the baseline expectation. We performed additional simulations with broad ranges of λ and δ values. In particular, we analyzed whether HGT would promote the likelihood of bistability in two-species communities compared with the scenario without gene transfer. Our results suggested that: (1) With or without HGT, reducing λ (increasing ecological neutrality) promoted bistability (Fig. 3g and h); (2) With HGT, increasing δ promoted bistability (Fig. 3h); (2) Compared with the population without HGT, gene transfer promoted bistability when δ was zero or positive, while reduces bistability when δ was largely negative (Fig. 3i). These results suggested that the interplay among HGT, species growth and interactions added new insights into the fundamental question of how microbes competing for a limited number of resources stably coexisted.

The role of HGT on community multistability when MGE epitasis effects are present

In previous analysis, we also assumed that the effects of MGEs on host growth (λ₁ and λ₂) were independent of the host species. In nature, however, the same MGE can have different growth rate effects in different genomic backgrounds, a phenomenon called epistasis^33-36. 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 µ₁ and µ₂ into two components: the basal growth rate ( and ), and the growth rate effects (λ₁₁ and λ₁₂) of the mobile genetic elements. The growth rates of p₁ and p₂ were calculated as µ₁(1 + λ₁₂) and µ₂(1 + λ₂₁), respectively. Here λ_ij (i, j = 1,2) stands for the growth rate effects of the j-th mobile genetic element in the i-th species. Without epistasis, the growth rate effects are independent of the host species (λ₁₁ = λ₂₁ and λ₂₂ = λ₁₂). The epistasis can be quantified by two ratios: ξ₁ = λ₂₁/λ₁₁ and ξ₂ = λ₁₂/λ₂₂. ξ₁ >0 and ξ₂ >0 represent magnitude epistasis, whereas ξ₁ < 0 or ξ₂ < 0 represent the sign epistasis. We carried out additional numerical simulations to evaluate how different epistasis types influenced the role of HGT in community stability. Our results suggested that with magnitude epistasis, HGT still enlarged the area of bistability region (Fig. 4a 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. 4c and d). These results suggest that MGE epistasis might add another layer of complexities into the interplay between HGT and community stability.

Fig. 4.

Interplay between HGT and the multistability of communities under strong selections

The growth rate effect of an MGE can be extreme in some cases. For instance, under strong antibiotic selection, only cells carrying the antibiotic resistance 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 MGE in a population under strong environmental selection (see Supplementary Materials for more details). Specifically, we assumed that the MGE was initially carried by one species in the populations. Without HGT, only the donor species carrying the MGE can survive due to selection (Fig. 5a). In this case, the system will always be monostable regardless of the initial species compositions. With HGT, the MGE can be transferred to other species, creating chances for other species to survive (Fig. 5b). To understand how HGT would affect 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. 5c and d). Here, without loss of generality we assumed that the MGE was initially carried only by species 1.

Fig. 5.

Bistability was more feasible when µ₂ >µ₁ (Fig. 5c), suggesting that the effect of HGT on promoting bistability would be stronger when the weak competitor acted as the MGE donor. We further analyzed the number of stable states in multispecies communities. As shown in Fig. 5e and f, for these populations under antibiotic selection, increasing transfer rate of the MGE also promoted the emergence of multistability. This conclusion equally held under different interaction strengths or growth rate ranges (Fig. S15), suggesting the general applicability of HGT’s role in promoting multistability for communities under strong selection.

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 in similar environments, also known as ‘patches’, form a metacommunity^37,38. When multiple stable states coexist, the species compositions of different patches may differ from each other due to the variability of their initial configurations. 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. S16a, see Methods for more details). Each local patch carried 10 competing species that exchanged mobile gene elements 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 the steady state.

Without HGT, the metacommunity exhibited strong spatial homogeneity at regional scale (Fig. S16b). 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. S16b). 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. S16b and c). The conclusion is applicable for different dispersal rates (Figs. S17 and S18). These results highlight the fundamental role of HGT in community stability and diversity at local and regional scale.

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 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 puzzle^39,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 (Fig. S3) and similar growth rates (Fig. S13). The relative importance of HGT in real environment is also dependent on the maximum carrying capacity of the population. Simulations by scanning the η value from 10⁻⁷ to 1 suggested that increasing HGT rates started to promote multistability when η value exceeded 10⁻² per hour (Fig. S9). This corresponds to a conjugation efficiency of 10⁻¹¹ cell⁻¹ · hr⁻¹ · mL when the maximum carrying capacity equals 10⁹ cells · mL⁻¹, or a conjugation efficiency of 10⁻¹⁴ cell⁻¹ · hr⁻¹· mL when the maximum carrying capacity equals 10¹² cells · mL⁻¹. Therefore, in environments with high cell density and abundant MGEs, such as mammal gut, the role of HGT will be more prominent⁴¹. In contrast, the influence of HGT on community multistability might become less important in populations with poor cell density and slow gene transfer.

Our work provided a comprehensive framework to explore how different parameters, including but not limited to HGT rates, collectively shaped community multistability. For instance, our simulations suggested that nonzero cell death or dilution rate was essential for the bistability of two competing species. Promoting competition strength in general enlarged multistability region. Besides, multistability was also dependent on community neutrality, i.e. the growth rate similarities among different competitors. Multistability is more feasible for communities with high neutrality. Which mechanisms contribute most to the multistability of a given community might be context-specific, depending on the values of different parameters in the population.

Our model allowed multiple plasmids to infect a single cell, but ignored the ecological interactions between different plasmids. These complex interactions are governed by many different processes: from systems that exclude the entry of other plasmids, incompatibility system that eliminates within-cell coexistence, to interference in the regulation of plasmid copy number^42,43. These interactions can further change the distribution of plasmids among different species, alter the structure of gene transfer network, and reshape the stability landscape of a microbial community. However, the general conclusion regarding how these interactions shape population multistability remains to be established in future studies.

Previous studies documented the presence of multiple alternative community types in soil¹⁵, human vaginal, gut and oral microbiomes^2,44,45. In particular, multiple groups of bacteria in human gut exhibit robust bistable abundance distributions, resulting in over 32 different stable state combinations⁴⁶. 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 factors¹. 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 disorders^47-49. Strategies have been developed to overcome the resilience of unhealthy states, including fecal transplantation or supplementation with probiotics¹⁸. 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 flow⁴¹. Previous studies suggested that some small molecules such as unsaturated fatty acids can inhibit the conjugative transfer of plasmids^50-52. By binding the type IV secretion traffic ATPase TrwD, these compounds limit the pilus biogenesis and DNA translocation⁵³. Therefore, modulating MGE spread using these molecules might offer new opportunities to reshape the stability landscape and narrow down the attraction domains of the disease states^51,54.

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 mobilizable plasmids with growth rate effects²⁸. Briefly, the model includes four ordinary differential equations (ODEs):

S₁ and S₂ are the abundances of the two species. R₁ and R₂ describe the strengths of interspecies competition. N_m,1 and N_m,2 are the maximum carrying capacities for species 1 and 2, respectively. D is the dilution rate. P_i (i = 1, 2) represents the abundance of the subpopulation in the i-th species that acquires the mobilizable genes from its competitor. Here, we assumed that the populations would not lose their own native plasmids completely. The internal inhibition within each species was also assumed to be independent of the acquired MGEs. µ_i represents the maximum ‘per capita’ 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 growth rate effect of the plasmids. Interspecies gene flow allows each species to acquire plasmids from its competitor, resulting in the dynamic change of the effective species growth rates , which is calculated as: and . For instance, all cells of species 1 carry the plasmid 1. Therefore, their growth rate can be obtained as . After HGT, the overall growth rate of species 1 changes because some of the individuals get the plasmid from the second species. Such a change is reflected by the λ₂ term in . and stand for HGT rates. Here, we described the HGT process by mass action kinetics, an assumption commonly made for plasmid transfer dynamics. We further assumed that intraspecies and interspecies transfer occurred at the same rates. κ₁ and κ₂ represent the loss rates of the mobile genetic elements. For plasmids, κ₁ and κ₂ describe the loss rate by segregation errors.

Let . The model can be further simplified as

Here, s₁, s₂, p₁ and p₂ are dimensionless. η_i is related with describes the ratio between N_m,2 and N_m,1: ϱ = N_m,2/N_m,1. For simplicity, we assumed that the two species had the same carrying capacity, i.e. ϱ = 1. The unit of η₁ and η₂ is hour⁻¹. Without HGT, p₁ and p₂ 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 ℊ₂ and ℊ₃ are both stable. The Jacobian matrix of ℊ₂ is . The eigenvalues of J₂ can be obtained as and D – µ₂. ℊ₂ is stable when both eigenvalues are negative, which leads to γ₂ >ϕ_2/ϕ₁ where and . Similarly, ℊ₃ is stable when γ₁ >ϕ_2/ϕ₁. Therefore, the condition for the system being bistable is γ₂ >ϕ₁/ϕ₂ and γ₁ >ϕ_2/ϕ₁.

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 s_j,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:

s_i is the abundance of the i-th species. p_ij is the size of the subpopulation in the i-th species that acquires s_i -originated MGEs. γ_ji is the negative interaction that s_j imposes on the i-th species. D is the dilution rate., the effective growth rate of s_i, is calculated by is defined as where is the basal growth rate, λ_ij is the growth rate effect of the s_j-originated MGEs in the i-th species. η_jki is the transfer rate of the s_j-originated plasmids from species k to species i. κ_ij is the loss rate of the j-th mobile genetic element in the i-th species.

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 < s_i < 1 following uniform distributions. The initial conditions for p_ij’s are p_ii = s_i and p_ij = 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.

Horizontal gene transfer within communities under strong environmental selection

Under strong antibiotic selections, only cells carrying antibiotic resistant MGEs can survive. To examine whether our conclusion is still applicable in this scenario, we generalized our model by considering the transfer of an antibiotic-resistant gene (ARG) in a population of m species under strong antibiotic selection. Within each species, only cells carrying the MGE can grow, with µ_i being their maximum growth rate. The cells without the MGE will be eventually depleted by dilution D. The community dynamics can be described by the following ODEs (i = 1,2, …, m):

s_i is the total abundance of the i-th species, and p_i is the abundance of cells carrying the MGE in s_i. In the first equation, the first term describes the overall growth of the i-th species, which is only contributed by p_i since the MGE-free cells are not able to grow. η_ji is the MGE transfer rate from the j-th to the i-th species. κ_i is MGE loss rate in the i-th species. γ_ji describes the inter-species competition strength.

Without loss of generality, we assumed the species 1 was the MGE donor and p₁ was equal with s₁. Without HGT, the MGE only resides in the donor species, while HGT allows the MGE to be shared with other species. To understand how HGT rate influences community multistability, we calculated the phase diagram of populations of two competing species. Increasing HGT rate enlarged the bistability region. We further analyzed the number of alternative stable states in communities of 5 or 8 species. In particular, we randomly initialized the species abundances 500 times, simulated the steady states and clustered them into different attractors. When calculating the number of stable states in multi-species communities, we randomly drew the µ_i values from a uniform distribution between 0.3 and 0.7 hr⁻¹. Our results suggested that the number of stable states in general increased with HGT rate.

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. M_i = k means that the i-th species lives in the k-th niche. For species i and j, if M_i = M_j, 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:

Here, s_i is the abundance of the i-th species. p_ij is the size of the subpopulation in the i-th species that acquires s_j-originated mobile genes. γ_ji is the negative interaction that s_j imposes on the i-th species. η_jki is the transfer rate of the s_j-originated genes from species k to species i. 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. We assumed a grid-like arrangement of patched connected by dispersal terms ℋ. 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], s_a,b,i is the abundance of the i-th species, and p_a,b,i,j is the size of the subpopulation in the i-th species that acquires s_j-originated mobile genes. The effective growth rate is calculated as .

and are the gains of s_a,b,i and p_a,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: s_i = ∑_a=1,2,…,u ∑_b=1,2,…,v s_a,b,j. Then we measured the diversity using Shannon index, a common metric used in ecology. Shannon index can be calculated as . s_p is the total species abundance: .

Biological relevance of our analysis

The empirically measured gene transfer rates (denoted as η^c, with the unit of cells⁻¹· mL · hour⁻¹) should be multiplied by N_m, the maximum carrying capacity of the population, before being used in our model. Therefore, the transfer rate η in our analysis are several orders of magnitude higher than the experimentally measured rates^55,56. In this work, we focused on HGT rate in the range of 0 < η_i < 0.4 hour⁻¹. Whether this range is biologically relevant depends on N_m of the population. Considering the typical estimates of plasmid transfer rates cross species (10⁻¹³ cells⁻¹ · mL · hour⁻¹, from Klebsiella pneumoniae to Escherichia coli) and within species (10⁻⁷ cells⁻¹ · mL · hour⁻¹, between E. coli strains) provided by a previous study⁵⁵, HGT rates in our analysis are relevant in various environments from human colon (∼10¹² cells per gram) ⁵⁷ to deep oceanic subsurface (10⁶∼10⁸ cells per gram)⁵⁸.

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.

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