Emergence of alternative stable states in microbial communities undergoing horizontal gene transfer

eLife Assessment

This manuscript offers valuable theoretical predictions on how horizontal gene transfer (HGT) can lead to alternative stable states in microbial communities. Using a modeling framework, solid theoretical evidence is provided to support the claimed role of HGT. However, given that the model has many degrees of freedom, a more comprehensive analysis of the role of different parameters could strengthen the study. Additionally, potential interactions between plasmids that carry out HGT are not discussed in the model. This paper would be of interest to researchers in microbiology, ecology, and evolutionary biology.

https://doi.org/10.7554/eLife.99593.3.sa0

Significance of the findings:

Valuable: Findings that have theoretical or practical implications for a subfield

Landmark
Fundamental
Important
Valuable
Useful

Strength of evidence:

Solid: Methods, data and analyses broadly support the claims with only minor weaknesses

Exceptional
Compelling
Convincing
Solid
Incomplete
Inadequate

During the peer-review process the editor and reviewers write an eLife Assessment that summarises the significance of the findings reported in the article (on a scale ranging from landmark to useful) and the strength of the evidence (on a scale ranging from exceptional to inadequate). Learn more about eLife Assessments

Abstract
Introduction
Results
Discussion
Methods
Appendix 1
Appendix 2
Data availability
References
Article and author information
Metrics

Abstract

Microbial communities living in the same environment often display alternative stable states, each characterized by a unique composition of species. Understanding the origin and determinants of microbiome multistability has broad implications in environments, human health, and microbiome engineering. However, despite its conceptual importance, how multistability emerges in complex communities remains largely unknown. Here, we focused on the role of horizontal gene transfer (HGT), one important aspect mostly overlooked in previous studies, on the stability landscape of microbial populations. Combining mathematical modeling and numerical simulations, we demonstrate that, when mobile genetic elements (MGEs) only affect bacterial growth rates, increasing HGT rate in general promotes multistability of complex microbiota. We further extend our analysis to scenarios where HGT changes interspecies interactions, microbial communities are subjected to strong environmental selections and microbes live in metacommunities consisting of multiple local habitats. We also discuss 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 undergoing HGT. These results reveal how different dynamic processes collectively shape community multistability and diversity. Our results provide key insights for the predictive control and engineering of complex microbiota.

Introduction

Microbial communities living in similar environments often display different community types, each characterized by a unique composition of species (Arumugam et al., 2011; Estrela et al., 2022; Gibson et al., 2016; Gonze et al., 2017). Changes in environmental factors, fluctuations of species abundances, or transient invasions can trigger abrupt and irreversible transitions between two types (Amor et al., 2020; Beardmore et al., 2018; Chaparro-Pedraza and de Roos, 2020). 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 (Fuhrman et al., 2015; Hartmann and Six, 2023; Kumar and Hasty, 2023; Relman, 2012; Sun et al., 2023). Understanding the origin of different community types has broad implications in environments, engineering, and health (Gibson et al., 2016; Gonze et al., 2017).

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 (Amor et al., 2020; Dubinkina et al., 2019; Estrela et al., 2022; Fujita et al., 2023; Gonze et al., 2017; Khazaei et al., 2020; Mangal et al., 2023; Zaoli and Grilli, 2021). 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’ (Gibson et al., 2016; Khazaei et al., 2020). Despite its conceptual importance, the determinants underlying microbiota multistability have been poorly understood due to the intrinsic complexity of microbial interactions (Dubinkina et al., 2019).

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 (Fassarella et al., 2021; Relman, 2012). 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 (van den Berg et al., 2022). 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 (Gonze et al., 2017). 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 (Dubinkina et al., 2019). 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 (Soucy et al., 2015).

Indeed, HGT is prevalent in microbial world and mobile genes occupy a substantial proportion in microbial gene pool (Bordenstein and Reznikoff, 2005; Brito et al., 2016; Brockhurst et al., 2019). For instance, among the closely related microbial isolates in the same human body site, HGT occurs in more than 40% of the species pairs (Smillie et al., 2011). In the draft genomes of human gut microbes, tens of thousands of genes are mobilizable (Brito et al., 2016). Depending on the encoded traits, MGEs shape microbial growth rates, host adaptability and even social interactions between microbes (Rankin et al., 2011). The dynamic gain and loss of MGEs are also a major source of evolutionary innovations of prokaryotic genomes (Werren, 2011). 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 (Dubinkina et al., 2019; Gilpin and Case, 1976; Gonze et al., 2017). 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 in 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 ( $μ_{1}$ and $μ_{2}$ ), interspecies interaction strengths ( $γ_{1}$ and $γ_{2}$ ) 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 (Figure 1A; Dubinkina et al., 2019; Gonze et al., 2017). 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 (Figure 1A). In contrast, if monostable, the population will always rest at the same steady state regardless of the initial composition (Figure 1A). 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 $ϕ_{1} = \frac{μ_{1} - D}{μ_{1}}$ and $ϕ_{2} = \frac{μ_{2} - D}{μ_{2}}$ (see Methods for more details). This condition predicted the thresholds on the strengths of interspecies interactions in order for the system to be bistable.

Figure 1 with 3 supplements see all

Download asset Open asset

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. The right panel is a diagram depicting how a microbial community responds to changes on initial species abundances. The red marble stands for the system state, which will always roll ‘downhill’ on the landscape towards stable states at the low points. The x axis represents the community composition and the valleys in the landscape are the stable states (or equilibrium points) where the community composition will stop changing. Depending on parameter values, the landscape can have one valley (called ‘monostable’) or two valleys (called ‘bistable’). When the system is bistable, the marble can reach either valley, depending on where it starts. (B) The phase diagrams of the population with (bottom panel, $η = 0.2 h r^{- 1}$ ) or without HGT (top panel). Different combinations of $μ_{1}$ and $μ_{2}$ between 0 and 1 hr⁻¹ were tested. For each combination, the initial abundance of each species was randomized 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. All combinations of $μ_{1}$ and $μ_{2}$ resulting in bistability were marked in blue in the diagram. Other parameters are $γ = 1.1$ , $κ = 0.005 h r^{- 1}$ , $D = 0.2 h r^{- 1}$ , $μ_{1}^{0} = μ_{2}^{0} = 0.5 h r^{- 1}$ . Two states, $P$ ( $μ_{1} = 0.3 h r^{- 1}$ , $μ_{2} = 0.8 h r^{- 1}$ ) and $Q$ ( $μ_{1} = 0.4 h r^{- 1}$ , $μ_{2} = 0.6 h r^{- 1}$ ), were labeled as examples. (**C and 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⁻¹ 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 the conjugative transfer of plasmids between competing species (Zhu et al., 2024). While in classic LV model $μ_{1}$ and $μ_{2}$ are constants, our model dissected species growth rates into two components: the basal growth rates ( $μ_{1}^{0}$ and $μ_{2}^{0}$ ) determined by chromosomal genes, and the growth rate effects ( $λ_{1}$ and $λ_{2}$ ) of plasmids. The two components were assumed to combine multiplicatively: $μ_{1} = μ_{1}^{0} (1 + λ_{1})$ , $μ_{2} = μ_{2}^{0} (1 + λ_{2})$ . HGT creates subpopulations (denoted as $p_{1}$ and $p_{2}$ ) within each species that acquire the mobilizable genes from its competitor. The dynamics of $p_{1}$ and $p_{2}$ are governed by cell growth, plasmid transfer (the rate denoted as $η$ ) and plasmid loss (the rate denoted as $κ$ ). The gain and loss of $p_{1}$ and $p_{2}$ in turn lead to the temporal change of the effective growth rates ( $μ_{1}^{e}$ and $μ_{2}^{e}$ ) of each species: $μ_{1}^{e} = μ_{1} (1 + λ_{2} \frac{p_{1}}{s_{1}})$ , $μ_{2}^{e} = μ_{2} (1 + λ_{1} \frac{p_{2}}{s_{2}})$ (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 (Gilpin and Case, 1976). 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. 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 (Figure 1A). We calculated the number of stable states for a wide range of $μ_{1}$ and $μ_{2}$ values and drew the phase diagram (Figure 1B). Whether it is easy for the system to achieve bistability is reflected by the total area of bistability region in the phase diagram. Our numerical results suggested that increasing HGT enlarged the region, creating bistability in many systems which would otherwise be monostable (Figure 1B and C and D). Here, the number of initializations was not critical for the conclusion (Figure 1—figure supplement 1A). The results were not fundamentally changed when the interspecies interaction strengths or gene transfer rates were altered (Figure 1—figure supplement 1B-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 (Figure 1C and D). For two species with large growth rate difference, the system might remain monostable when gene transfer rate increases (Figure 1C). In this case, the strong competitor will always win despite the gene flow. The effect of HGT in promoting bistability is more evident between species with smaller growth rate difference (Figure 1B and D).

Previous studies showed that the cell death rate determines the bistability of two interacting species (Abreu et al., 2019). 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 (Figure 1—figure supplement 2). Our results suggested that varying death rate indeed changed the bistability of the system. When the death rate equaled zero, MGEs 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 interactions also shaped the community bistability. By calculating the bistability probability under different values of $γ$ (Figure 1—figure supplement 3), 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 populations 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. $γ_{i j}$ represents the interspecies interaction strength that the $i$ -th species imposes on the j-th species ( $i, j = 1,2, . . ., m$ ). $η_{j k i}$ is the transfer rate of the $s_{j}$ -originated plasmid from species $k$ to species $i$ . $κ_{i j}$ is the loss rate of the $j$ -th plasmid in the $i$ -th species. The growth rate effect of the $j$ -th plasmid on the $i$ -th species is denoted as $λ_{i j}$ . $p_{i j}$ 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: $μ_{i}^{e} = μ_{i} \prod_{j \neq i} (1 + λ_{i j} \frac{p_{i j}}{s_{i}})$ . The population dynamics can then be simulated by $m + m^{2}$ equations that describe the temporal change of $s_{i}$ and $p_{i j}$ , respectively (see Methods for more details).

To understand how changing HGT 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 hrs till the systems reached equilibrium and calculated the steady-state abundances of different species (Figure 2—figure supplement 1). Limit cycles were not observed in our simulations (Figure 2—figure supplement 2).

To 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 (Figure 2—figure supplement 3). 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 ( $\sum_{i} x_{i} = 1$ , $i = 1,2, \dots, φ$ ). $χ$ is calculated as $χ = e x p (- \sum_{i} x_{i} l o g x_{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 (Figure 2—figure supplement 4).

As shown in Figure 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 (Figure 2A and B). Each stable state is characterized by a different dominant species (Figure 2A). The number of stable states is positively associated with HGT rate (Figure 2C; Figure 2—figure supplements 5 and 6). This prediction is generally applicable to communities consisting of different number of species (Figure 2—figure supplements 5 and 6). For communities of two species, bistability is not possible when $γ_{i j}$ ’s are smaller than 1. However, for communities of multiple species, multistability becomes possible even if all $γ_{i j}$ values are smaller than 1 (Figure 2—figure supplement 7). The effects of HGT on multistability held for different strengths of interspecies interactions (Figure 2—figure supplement 7). Initializing the $p_{i j}$ subpopulations with non-zero abundances did not change the conclusion, either (Figure 2—figure supplement 8).

Figure 2 with 11 supplements see all

Download asset Open asset

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.4 h r^{- 1}$ ), four different attractors appear. Here, populations reaching different stable states were marked with different colors. Other parameters are $μ_{1}^{0} = μ_{2}^{0} = μ_{3}^{0} = μ_{4}^{0} = μ_{5}^{0} = 0.5 h r^{- 1}$ , $μ_{1} = 0.52 h r^{- 1}$ , $μ_{2} = 0.48 h r^{- 1}$ , $μ_{3} = 0.44 h r^{- 1}$ , $μ_{4} = 0.60 h r^{- 1}$ , $μ_{5} = 0.31 h r^{- 1}$ , $γ = 1.05$ , $κ = 0.005 h r^{- 1}$ , $D = 0.2 h r^{- 1}$ . (C) The number of stable states (top panel) and the multistability coefficient $χ$ increase with HGT rate in communities consisting of eight 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.

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 (Figure 2—figure supplement 9). In particular, we randomly sampled $γ_{i j}$ values from uniform distributions. Our results suggested the mean value and variance of $γ_{i j}$ played a role in shaping multistability. The effects of HGT on community multistability became stronger when the mean value of $γ_{i j}$ got larger than 1 and the variance of $γ_{i j}$ was small (Figure 2—figure supplement 9). Nevertheless, the heterogeneous distribution of $γ_{i j}$ 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 $λ_{i j}$ magnitude, we also analyzed different ranges of growth rates effects of MGEs, by sampling $λ_{i j}$ values from uniform distributions with given widths (Figure 2—figure supplement 10). 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 (Figure 2—figure supplement 10B). The qualitative relationship between HGT and community was not dependent on the range of growth rate effects (Figure 2—figure supplement 10A).

So far, our analysis has 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 many species. While this analysis provides the basic proof of principle, natural environments often consist of many niches (Finlay and Medzhitov, 2007). Species living in the same niche compete while species from different niches can coexist (Figure 2—figure supplement 11A). 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 communities (Figure 2—figure supplement 11B-D).

The role of HGT in community multistability when MGEs 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 MGEs.

In nature, the horizontal transfer of some genes can change the interspecies interaction strengths. For instance, the sharing of many mobile genes can promote niche overlapping, leading to an increase of competition strength (Bonham et al., 2017; Caro-Quintero and Konstantinidis, 2015). 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 Appendix 1 for more details). For two-species populations, we divided the interspecies interactions 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 $γ_{1} + δ_{1} \frac{p_{1}}{s_{1}}$ and $γ_{2} + δ_{2} \frac{p_{2}}{s_{2}}$ , respectively. Positive $δ$ ’s mean HGT promotes the interspecies competition, while negative $δ$ ’s mean HGT reduces the strength of interspecies competition (Figure 3A and 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 (Figure 3C and D). We further extended our analysis to multispecies communities and calculated the number of stable states. As shown in Figure 3E and F, the number of stable states increased with HGT rate when MGEs promoted competition, while decreased with HGT rate when MGEs reduced competition.

Figure 3

Download asset Open asset

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

(**A and B**) The schematics of HGT promoting or reducing competition. (C) For populations of two species, when MGEs promote competition, increasing HGT rate enlarges 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 $μ_{1}^{0}$ and $μ_{2}^{0}$ 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$ , $κ = 0.005 h r^{- 1}$ , $D = 0.2 h r^{- 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 five 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$ , $κ = 0.005 h r^{- 1}$ , $D = 0.2 {h r}^{- 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. (G) For populations of two species without HGT, reducing the magnitude of $λ$ (promoting ecological neutrality) enlarges the area of bistability. When gene transfer rate equals zero, changing $δ$ does not influence bistability. When calculating bistability probability, we randomized $λ_{1}$ and $λ_{2}$ 500 times between $- α$ and $α$ following uniform distributions. $α$ represents the magnitude of $λ$ variations. $μ_{1}$ and $μ_{2}$ were calculated as $μ_{1}^{0} (1 + λ_{1})$ and $μ_{2}^{0} (1 + λ_{2})$ , respectively. The overall competition strength was calculated as $γ_{1} + δ_{1} \frac{p_{1}}{s_{1}}$ and $γ_{2} + δ_{2} \frac{p_{2}}{s_{2}}$ . We tested different combinations of $α$ (the magnitude of $λ$ variations) and $δ$ values. We then calculated the fraction ( $f_{0}$ ) of growth rate combinations that generated bistability out of the 500 random combinations. Other parameters are $γ_{1} = γ_{2} = 1.1$ , $κ = 0.005 h^{- 1}$ , $D = 0.2 h^{- 1}$ . (H) In populations of two species with HGT, reducing the magnitude of $λ$ (promoting ecological neutrality) or promoting $δ$ value enlarges the area of bistability ( $f_{1}$ ). $η$ being 0.4 was tested as an example here. (I) The effect of HGT on community bistability depends on how mobile genes modify growth rates and competition. Compared with the population without HGT, gene transfer promotes bistability when $δ$ is zero or positive, while reduces bistability when $δ$ is largely negative.

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 modeling 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 $λ$ , and the influence on competition is described by the parameter $δ$ (Figure 3G, H and I). By varying these two parameters, we can evaluate how the modeling 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 (Figure 3G and H); (2) With HGT, increasing $δ$ promoted bistability (Figure 3H); (2) Compared with the population without HGT, gene transfer promoted bistability when $δ$ was zero or positive, while reduces bistability when $δ$ was largely negative (Figure 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 in community multistability when MGE epitasis effects are present

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 growth rate effects in different genomic backgrounds, a phenomenon called epistasis (Acar Kirit et al., 2020; Gama et al., 2020; San Millan et al., 2014; Silva et al., 2011). 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 Appendix 2 for more details). In particular, we dissected the growth rates $μ_{1}$ and $μ_{2}$ into two components: the basal growth rate ( $μ_{1}^{0}$ and $μ_{2}^{0}$ ), and the growth rate effects ( $λ_{11}$ and $λ_{22}$ ) of the MGEs. The growth rates of $p_{1}$ and $p_{2}$ were calculated as $μ_{1} (1 + λ_{12})$ and $μ_{2} (1 + λ_{21})$ , respectively. Here $λ_{i j}$ ( $i, j = 1,2$ ) stands for the growth rate effects of the $j$ -th MGE in the $i$ -th species. Without epistasis, the growth rate 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 influenced the role of HGT in community stability. Our results suggested that with magnitude epistasis, HGT still enlarged the area of bistability region (Figure 4A & 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 (Figure 4C & D). These results suggest that MGE epistasis might add another layer of complexities into the interplay between HGT and community stability.

Figure 4

Download asset Open asset

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

(**A and B**) With magnitude epistasis, increasing HGT rate still promotes bistability in populations of two competing species. (**C and 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 growth rate 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$ , $κ = 0.005 h r^{- 1}$ , $D = 0.2 h r^{- 1}$ .

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 Methods for 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 (Figure 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 (Figure 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 (Figure 5C and D). Here, without loss of generality we assumed that the MGE was initially carried only by species 1. Bistability was more feasible when $μ_{2} > μ_{1}$ (Figure 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 Figure 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 (Figure 5—figure supplement 1), suggesting the general applicability of HGT’s role in promoting multistability for communities under strong selection.

Figure 5 with 4 supplements see all

Download asset Open asset

HGT creates chances of multistability for communities under strong environmental selection.

(**A and B**) The schematics of communities under strong selection. Without HGT, only the donor species carrying the MGE can survive. HGT allows the other species to acquire MGE from the donor, creating opportunities for the other species to survive under strong selections. (C) The phase diagram of two-species populations transferring an MGE. Different combinations of $μ_{1}$ and $μ_{2}$ between 0 and 1 hr⁻¹ were tested. All combinations of $μ_{1}$ and $μ_{2}$ that led to bistability were marked in blue in the diagram. Other parameters are $γ = 1.1$ , $κ = 0.005 h r^{- 1}$ , $D = 0.2 h r^{- 1}$ , $μ_{1}^{0} = μ_{2}^{0} = 0.5 h r^{- 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 MGE. Three different competition strengths were tested (marked in different colors). (**E and F**) For communities of 5 or 8 species under strong selection, increasing HGT rate 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 carriage of the mobile gene changed the species growth rate from 0 to a positive value $μ_{i}$ . 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⁻¹. 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$ , $κ = 0.005 h r^{- 1}$ , $D = 0.2 h r^{- 1}$ .

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 (Leibold et al., 2004; Miller et al., 2018). 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 matrix of local patches and assumed that the dispersal only occurred between adjacent patches (Figure 5—figure supplement 2A, see Methods for more details). Each local patch carried 10 competing species that exchanged MGEs 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 (Figure 5—figure supplement 2B). 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 (Figure 5—figure supplement 2B). 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 (Figure 5—figure supplement 2B and C). The conclusion is applicable for different dispersal rates (Figure 5—figure supplements 3 and 4). 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 (Grilli et al., 2017; Serván et al., 2018). 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 specific community might also be context dependent. The effect of HGT in promoting multistability might be stronger between species with intense competition (Figure 1—figure supplement 3) and similar growth rates (Figure 2—figure supplement 10). The relative importance of HGT is also dependent on the maximum carrying capacity of the population. Simulations by scanning the $η$ value from 10^–7 to 1 suggested that increasing HGT rates started to promote multistability when $η$ value exceeded 10^–2 per hour (Figure 2—figure supplement 6). This corresponds to a conjugation efficiency of 10^–11 cells. hr⁻¹.mL when the maximum carrying capacity equals 10⁹ cells. mL, or a conjugation efficiency of 10^–14 cells. 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 (Liu et al., 2012). 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, that is 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 (Igler et al., 2022; Lipsitch et al., 2009). 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 (Fujita et al., 2023), human vaginal, gut, and oral microbiomes (Arumugam et al., 2011; Ding and Schloss, 2014; Ravel et al., 2011). In particular, multiple groups of bacteria in human gut exhibit robust bistable abundance distributions, resulting in over 32 different stable state combinations (Lahti et al., 2014). Our results provide a plausible explanation for the large number of community types in these populations, even though our interpretation does not exclude other mechanisms including gradients of environmental factors (Gonze et al., 2017). 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 (Canfora et al., 2019; Kootte et al., 2017; Moustafa et al., 2018). Strategies have been developed to overcome the resilience of unhealthy states, including fecal transplantation or supplementation with probiotics (Fassarella et al., 2021). 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 (Liu et al., 2012). Previous studies suggested that some small molecules such as unsaturated fatty acids can inhibit the conjugative transfer of plasmids (Fernandez-Lopez et al., 2005; Getino and de la Cruz, 2018; Getino et al., 2015). By binding the type IV secretion traffic ATPase TrwD, these compounds limit the pilus biogenesis and DNA translocation (Ripoll-Rozada et al., 2016). 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 (Fernandez-Lopez et al., 2005; Palencia-Gándara et al., 2021).

Methods

Key resources table

Reagent type (species) or resource	Designation	Source or reference	Identifiers	Additional information
Software, algorithm	Codes for all numerical simulations	Github repository Hong et al., 2025	https://doi.org/10.5281/zenodo.14928352

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 (Zhu et al., 2024). Briefly, the model includes four ordinary differential equations (ODEs):

\frac{d S_{1}}{d t} = μ_{1}^{e} S_{1} (1 - \frac{S_{1} + R_{2} S_{2}}{N_{m, 1}}) - D S_{1},

\frac{d S_{2}}{d t} = μ_{2}^{e} S_{2} (1 - \frac{S_{2} + R_{1} S_{1}}{N_{m, 2}}) - D S_{2},

\frac{d P_{1}}{d t} = μ_{1} (1 + λ_{2}) P_{1} (1 - \frac{S_{1} + R_{2} S_{2}}{N_{m, 1}}) + η_{1}^{c} (S_{2} + P_{1}) (S_{1} - P_{1}) - (D + κ_{1}) P_{1},

\frac{d P_{2}}{d t} = μ_{2} (1 + λ_{1}) P_{2} (1 - \frac{S_{2} + R_{1} S_{1}}{N_{m, 2}}) + η_{2}^{c} (S_{1} + P_{2}) (S_{2} - P_{2}) - (D + κ_{2}) P_{2} .

$S_{1}$ and $S_{2}$ are the abundances of the two species. $R_{1}$ and $R_{2}$ 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 $μ_{i} = μ_{i}^{0} (1 + λ_{i})$ where $μ_{i}^{0}$ 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 $μ_{i}^{e}$ , which is calculated as: $μ_{1}^{e} = μ_{1} (1 + λ_{2} \frac{P_{1}}{S_{1}})$ and $μ_{2}^{e} = μ_{2} (1 + λ_{1} \frac{P_{2}}{S_{2}})$ . For instance, all cells of species 1 carry the plasmid 1. Therefore, their growth rate can be obtained as $μ_{1}^{0} (1 + λ_{1})$ . 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 $λ_{2}$ term in $μ_{1}^{e} = μ_{1} (1 + λ_{2} \frac{P_{1}}{S_{1}})$ . $η_{1}^{c}$ and $η_{2}^{c}$ 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. $κ_{1}$ and $κ_{2}$ represent the loss rates of the mobile genetic elements. For plasmids, $κ_{1}$ and $κ_{2}$ describe the loss rate by segregation errors.

Let $s_{1} = S_{1} / N_{m, 1},$ , $s_{2} = S_{2} / N_{m, 2}$ , $p_{1} = P_{1} / N_{m, 1}$ , $p_{2} = P_{2} / N_{m, 2}$ , $γ_{2} = R_{2} \frac{N_{m, 2}}{N_{m, 1}}$ , $γ_{1} = R_{1} \frac{N_{m, 1}}{N_{m, 2}}$ . The model can be further simplified as

\frac{d s_{1}}{d t} = μ_{1}^{e} s_{1} (1 - s_{1} - γ_{2} s_{2}) - D s_{1},

\frac{d s_{2}}{d t} = μ_{2}^{e} s_{2} (1 - s_{2} - γ_{1} s_{1}) - D s_{2},

\frac{d p_{1}}{d t} = μ_{1} (1 + λ_{2}) p_{1} (1 - s_{1} - γ_{2} s_{2}) + η_{1} (ϱ s_{2} + p_{1}) (s_{1} - p_{1}) - (D + κ_{1}) p_{1},

\frac{d p_{2}}{d t} = μ_{2} (1 + λ_{1}) p_{2} (1 - s_{2} - γ_{1} s_{1}) + η_{2} (ϱ^{- 1} s_{1} + p_{2}) (s_{2} - p_{2}) - (D + κ_{2}) p_{2} .

Here, $s_{1}$ , $s_{2}$ , $p_{1}$ and $p_{2}$ are dimensionless. $η_{i}$ is related with $η_{i}^{c}$ by $η_{i} = N_{m, i} η_{i}^{c}$ . $ϱ$ 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, that is $ϱ = 1$ . The unit of $η_{1}$ and $η_{2}$ is hour⁻¹. Without HGT, $p_{1}$ and $p_{2}$ become consistently 0, and the model is equivalent to the classic Lotka-Volterra (LV) framework:

\frac{d s_{1}}{d t} = μ_{1} s_{1} (1 - s_{1} - γ_{2} s_{2}) - D s_{1},

\frac{d s_{2}}{d t} = μ_{2} s_{2} (1 - s_{2} - γ_{1} s_{1}) - D s_{2} .

Stability analysis of the classic LV model

The classic LV model has four fixed points: $g_{1} = [0, 0]$ , $g_{2} = [0, 1 - \frac{D}{μ_{2}}]$ , $g_{3} = [1 - \frac{D}{μ_{1}}, 0]$ , $g_{4} = [\frac{γ_{2} μ_{1} μ_{2} - D γ_{2} μ_{1} - μ_{1} μ_{2} + D μ_{2}}{μ_{1} μ_{2} (γ_{1} γ_{2} - 1)}, \frac{γ_{1} μ_{1} μ_{2} - D γ_{1} μ_{2} - μ_{1} μ_{2} + D μ_{1}}{μ_{1} μ_{2} (γ_{1} γ_{2} - 1)}]$ . The Jacobian matrix of each fixed point can be calculated as $J = [\begin{array}{cc} μ_{1} - 2 μ_{1} s_{1}^{*} - μ_{1} γ_{2} s_{2}^{*} D & - μ_{1} γ_{2} s_{1}^{*} \\ - μ_{2} γ_{1} s_{2}^{*} & μ_{2} - 2 μ_{2} s_{2}^{*} - μ_{2} γ_{1} s_{1}^{*} - D \end{array}]$ , where $s_{1}^{*}$ and $s_{2}^{*}$ are the species abundances at each fixed point.

The system is bistable only when $g_{2}$ and $g_{3}$ are both stable. The Jacobian matrix of $g_{2}$ is $J_{2} = [\begin{array}{cc} μ_{1} - D - μ_{1} γ_{2} + \frac{μ_{1} γ_{2} D}{μ_{2}} & 0 \\ - μ_{2} γ_{1} + D γ_{1} & D - μ_{2} \end{array}]$ . The eigenvalues of $J_{2}$ can be obtained as $μ_{1} - D - μ_{1} γ_{2} + \frac{μ_{1} γ_{2} D}{μ_{2}}$ and $D - μ_{2}$ . $g_{2}$ is stable when both eigenvalues are negative, which leads to $γ_{2} > ϕ_{1} / ϕ_{2}$ where $ϕ_{1} = \frac{μ_{1} - D}{μ_{1}}$ and $ϕ_{2} = \frac{μ_{2} - D}{μ_{2}}$ . Similarly, $g_{3}$ is stable when $γ_{1} > ϕ_{2} / ϕ_{1}$ . Therefore, the condition for the system being bistable is $γ_{1} > ϕ_{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 $d_{i j} = \sqrt{\sum_{k = 1}^{m} {(s_{i, k} - s_{j, k})}^{2}}$ where $k$ is the species index and $s_{i, 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:

\frac{d s_{i}}{d t} = μ_{i}^{e} s_{i} (1 - s_{i} - \sum_{j \neq i} γ_{j i} s_{j}) - D s_{i},

\begin{array}{ll} \frac{d p_{i j}}{d t} = μ_{i} (1 + λ_{i j}) [\prod_{k \neq i, j} (1 + λ_{i k} \frac{p_{i k}}{s_{i}})] p_{i j} (1 - s_{i} - \sum_{j \neq i} γ_{j i} s_{j}) + (s_{i} - p_{i j}) \sum_{k = 1}^{n} η_{j k i} p_{k j} \\ - (D + κ_{i j}) p_{i j} . (i \neq j) . \end{array}

$s_{i}$ is the abundance of the $i$ -th species. $p_{i j}$ is the size of the subpopulation in the $i$ -th species that acquires $s_{j}$ -originated MGEs. $γ_{j i}$ is the negative interaction that $s_{j}$ imposes on the $i$ -th species. $D$ is the dilution rate. $μ_{i}^{e}$ , the effective growth rate of $s_{i}$ , is calculated by $μ_{i}^{e} = μ_{i} \prod_{j \neq i} (1 + λ_{i j} \frac{p_{i j}}{s_{i}})$ . $μ_{i}$ is defined as $μ_{i} = μ_{i}^{0} (1 + λ_{i i})$ where $μ_{i}^{0}$ is the basal growth rate, $λ_{i j}$ is the growth rate effect of the $s_{j}$ -originated MGEs in the $i$ -th species. $η_{j k i}$ is the transfer rate of the $s_{j}$ -originated plasmids from species $k$ to species $i$ . $κ_{i j}$ 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_{i j}$ ’s are $p_{i i} = s_{i}$ and $p_{i j} = 0$ ( $i \neq j$ ). Then the dynamics of each population was simulated for 2000 hr 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, \dots, m$ ):

\frac{d s_{i}}{d t} = μ_{i} p_{i} (1 - s_{i} - \sum_{j \neq i} γ_{j i} s_{j}) - D s_{i},

\frac{d p_{i}}{d t} = μ_{i} p_{i} (1 - s_{i} - \sum_{j \neq i} γ_{j i} s_{j}) + (s_{i} - p_{i}) \sum_{j = 1}^{m} η_{j i} p_{j} - κ_{i} p_{i} - D p_{i} .

$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. $η_{j i}$ is the MGE transfer rate from the $j$ -th to the $i$ -th species. $κ_{i}$ is MGE loss rate in the $i$ -th species. $γ_{j i}$ describes the inter-species competition strength.

Without loss of generality, we assumed the species 1 was the MGE donor and $p_{1}$ was equal with $s_{1}$ . 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 $N_{m}^{(i)}$ ( $i = 1,2, \dots, l$ ). 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:

\frac{d s_{i}}{d t} = μ_{i}^{e} s_{i} (1 - \frac{s_{i} + \sum_{{j : M_{j} = M_{i}}} γ_{j i} s_{j}}{N_{m}^{(M_{i})}}) - D s_{i},

\begin{array}{ll} \frac{d p_{i j}}{d t} = μ_{i} (1 + λ_{i j}) [\prod_{k \neq i, j} (1 + λ_{i k} \frac{p_{i k}}{s_{i}})] P_{i j} (1 - \frac{s_{i} + \sum_{{j : M_{j} = M_{i}}} γ_{j i} s_{j}}{N_{m}^{(M_{i})}}) + (s_{i} - p_{i j}) \sum_{k = 1}^{n} η_{j k i} p_{k j} \\ - (D + κ_{i j}) p_{i j} . (i \neq j) . \end{array}

Here, $s_{i}$ is the abundance of the $i$ -th species. $p_{i j}$ is the size of the subpopulation in the $i$ -th species that acquires $s_{j}$ -originated mobile genes. $γ_{j i}$ is the negative interaction that $s_{j}$ imposes on the $i$ -th species. $η_{j k i}$ 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 \times v$ local patches. We assumed a grid-like arrangement of patched connected by dispersal terms $H$ . The location of each patch can be described by two coordinates $[a, b]$ ( $a = 1,2, \dots, u$ and $b = 1,2, \dots, 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:

\frac{d s_{a, b, i}}{d t} = μ_{a, b, i}^{e} s_{a, b, i} (1 - s_{a, b, i} - \sum_{j \neq i} γ_{j i} s_{a, b, j}) - D s_{a, b, i} + H_{s}^{(a, b)},

\begin{array}{ll} \frac{d p_{a, b, i, j}}{d t} = μ_{i} (1 + λ_{i j}) [\prod_{k \neq i, j} (1 + λ_{i k} \frac{p_{a, b, i, k}}{s_{a, b, i}})] p_{a, b, i, j} (1 - s_{a, b, i} - \sum_{j \neq i} γ_{j i} s_{a, b, j}) \\ + (s_{a, b, i} - p_{a, b, i, j}) \sum_{k = 1}^{n} η_{j k i} p_{a, b, k, j} - (D + κ_{i j}) p_{a, b, i, j} + H_{p}^{(a, b)} . (i \neq j) . \end{array}

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 $μ_{a, b, i}^{e}$ is calculated as $μ_{a, b, i}^{e} = μ_{i} \prod_{j \neq i} (1 + λ_{i j} \frac{p_{a, b, i, j}}{s_{a, b, i}})$ .

$H_{s}^{(a, b)}$ and $H_{p}^{(a, b)}$ 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, $H_{s}^{(a, b)}$ was calculated by summing the dispersal terms at four directions: left, right, up, and down. The left component can be calculated as

H_{s, l e f t}^{(a, b)} = {\begin{cases} ω (s_{a - 1, b, i} - s_{a, b, i}), i f a - 1 > 0 \\ 0, i f a - 1 \leq 0 \end{cases} .

Similarly, $H_{s, r i g h t}^{(a, b)}$ , $H_{s, u p}^{(a, b)}$ and $H_{s, d o w n}^{(a, b)}$ is calculated as:

H_{s, r i g h t}^{(a, b)} = {\begin{cases} ω (s_{a + 1, b, i} - s_{a, b, i}), i f a + 1 < u + 1 \\ 0, i f a + 1 \geq u + 1 \end{cases},

H_{s, u p}^{(a, b)} = {\begin{cases} ω (s_{a, b + 1, i} - s_{a, b, i}), i f b + 1 < v + 1 \\ 0, i f b + 1 \geq v + 1, \end{cases}

H_{s, d o w n}^{(a, b)} = {\begin{cases} ω (s_{a, b - 1, i} - s_{a, b, i}), i f b - 1 > 0 \\ 0, i f b - 1 \leq 0 \end{cases} .

$H_{s}^{(a, b)}$ can then be obtained as $H_{s, l e f t}^{(a, b)} + H_{s, r i g h t}^{(a, b)} + H_{s, u p}^{(a, b)} + H_{s, d o w n}^{(a, b)}$ . $H_{p}^{(a, b)}$ 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} = \sum_{a = 1,2, \dots, u} \sum_{b = 1,2, \dots, v} s_{a, b, i}$ . Then we measured the diversity using Shannon index, a common metric used in ecology. Shannon index can be calculated as $e x p [- \sum_{i = 1}^{m} \frac{s_{i}}{s_{T}} l o g \frac{s_{i}}{s_{T}}]$ . $s_{T}$ is the total species abundance: $s_{T} = \sum_{i = 1}^{m} s_{i}$ .

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 (Kosterlitz et al., 2022; Lopatkin et al., 2017). In this work, we focused on HGT rate in the range of $0 < η_{i} < 0.4 h o u r^{- 1}$ . 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 (Kosterlitz et al., 2022), HGT rates in our analysis are relevant in various environments from human colon (~10¹² cells per gram; Dieterich et al., 2018) to deep oceanic subsurface (10⁶~10⁸ cells per gram; Flemming and Wuertz, 2019).

Appendix 1

The effects of HGT on population multistability when MGEs change interspecies competitions

In nature, the HGT of some genes can change the competition strengths between different species. To understand how the transfer of these genes would influence community multistability, we adapted the main model, by accounting for the dynamic change of competition strength during gene transfer. For a population of two competing species, the model consists of four ODEs:

\frac{d s_{1}}{d t} = μ_{1} s_{1} [1 - s_{1} - (γ_{2} + δ_{2} \frac{p_{2}}{s_{2}}) s_{2}] - D s_{1},

\frac{{d s}_{2}}{d t} = μ_{2} s_{2} [1 - s_{2} - (γ_{1} + δ_{1} \frac{p_{1}}{s_{1}}) s_{1}] - D s_{2},

\frac{d p_{1}}{d t} = μ_{1} p_{1} [1 - s_{1} - (γ_{2} + δ_{2} \frac{p_{2}}{s_{2}}) s_{2}] + η_{1} (s_{2} + p_{1}) (s_{1} - p_{1}) - (D + κ_{1}) p_{1},

\frac{{d p}_{2}}{d t} = μ_{2} p_{2} [1 - s_{2} - (γ_{1} + δ_{1} \frac{p_{1}}{s_{1}}) s_{1}] + η_{2} (s_{1} + p_{2}) (s_{2} - p_{2}) - (D + κ_{2}) p_{2} .

Here, we assumed that mobile genes didn’t cause growth burden or benefits on the species growth rates. The interspecies competitions are determined by two components: the basal competition strengths ( $γ_{1}$ and $γ_{2}$ ) and the added parts by HGT ( $δ_{1}$ and $δ_{2}$ ). The overall competition strengths are calculated as $γ_{1} + δ_{1} \frac{p_{1}}{s_{1}}$ and $γ_{2} + δ_{2} \frac{p_{2}}{s_{2}}$ , respectively. Positive $δ$ values represent HGT promoting competition, while negative $δ$ values represent HGT reducing competition.

This framework can be extended to populations composed of multiple species. For a community of $m$ species. The population dynamics can be described by $2 m$ ODEs:

\frac{d s_{i}}{d t} = μ_{i} s_{i} [1 - s_{i} - \sum_{j \neq 1} (γ_{j i} + δ_{j i} \frac{p_{j i}}{s_{j}} s_{j}) - D s_{i}],

\frac{d p_{i j}}{d t} = μ_{i} p_{i j} [1 - s_{i} - \sum_{j \neq i} (γ_{j i} + δ_{j i} \frac{p_{j i}}{s_{j}}) s_{j}] + (s_{i} - p_{i j}) \sum_{k = 1}^{m} η_{j k i} p_{k j} - (D + κ_{i j}) p_{i j} . (i \neq j) .

$s_{i}$ represents the abundance of the $i$ -th species, and $p_{i j}$ represents the abundance of cells in the $i$ -th species that acquires $s_{j}$ -originated mobile genetic elements. We assumed $p_{i i} = s_{i}$ . $μ_{i}$ is the maximum growth rate of $s_{i}$ . $γ_{j i}$ stands for the basal interaction that $s_{j}$ imposes on the $i$ -th species without HGT. $δ_{j i}$ describes the change of inter-species competition due to gene transfer. The overall interaction that $s_{j}$ imposes on the $i$ -th species becomes $γ_{j i} + δ_{j i} \frac{p_{j i}}{s_{j}}$ . $η_{j k i}$ is the transfer rate of the $s_{j}$ -originated MGEs from species $k$ to species $i$ . $D$ and $κ_{i j}$ are the dilution and MGE loss rate, respectively.

How HGT affects bistability depends on the influences of mobile genes on 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 MGEs simultaneously modify growth rates and competitions:

\frac{{d s}_{1}}{d t} = μ_{1}^{e} s_{1} [1 - s_{1} - (γ_{2} + δ_{2} \frac{p_{2}}{s_{2}}) s_{2}] - D s_{1},

\frac{{d s}_{2}}{d t} = μ_{2}^{e} s_{2} [1 - s_{2} - (γ_{1} + δ_{1} \frac{p_{1}}{s_{1}}) s_{1}] - D s_{2},

\frac{{d p}_{1}}{d t} = μ_{1} (1 + λ_{2}) p_{1} [1 - s_{1} - (γ_{2} + δ_{2} \frac{p_{2}}{s_{2}}) s_{2}] + η_{1} (s_{2} + p_{1}) (s_{1} - p_{1}) - (D + κ_{1}) p_{1},

\frac{d p_{2}}{d t} = μ_{2} (1 + λ_{1}) p_{2} [1 - s_{2} - (γ_{1} + δ_{1} \frac{p_{1}}{s_{1}}) s_{1}] + η_{2} (s_{1} + p_{2}) (s_{2} - p_{2}) - (D + κ_{2}) p_{2} .

Here, $μ_{1}$ and $μ_{2}$ are calculated as $μ_{1} = μ_{1}^{0} (1 + λ_{1})$ and $μ_{2} = μ_{2}^{0} (1 + λ_{2})$ , respectively. $μ_{1}^{0}$ and $μ_{2}^{0}$ are the basal growth rates determined by the non-mobilizable genes. $λ_{1}$ and $λ_{2}$ are the growth rate effects of the MGEs. The effective species growth rates are calculated as: $μ_{1}^{e} = μ_{1} (1 + λ_{2} \frac{P_{1}}{S_{1}})$ and $μ_{2}^{e} = μ_{2} (1 + λ_{1} \frac{P_{2}}{S_{2}})$ . By changing the $δ$ values and the magnitudes of $λ$ ’s, this model structure allowed us to evaluate how the interplay among HGT, species growth and interactions affected the baseline expectation.

Appendix 2

Modeling the epistasis of mobile genetic elements

To evaluate the influence of MGE epistasis on the role of HGT in mediating population bistability, we dissected the growth rates $μ_{1}$ and $μ_{2}$ into two components: the basal growth rate ( $μ_{1}^{0}$ and $μ_{2}^{0}$ ), and the growth rate effects ( $λ_{11}$ and $λ_{22}$ ) of the mobilizable genes. The growth rates of $p_{1}$ and $p_{2}$ were calculated as $μ_{1} (1 + λ_{12})$ and $μ_{2} (1 + λ_{21})$ , respectively. Here, $λ_{i j}$ ( $i, j = 1,2$ ) describes the growth rate effects of the $j$ -th mobile gene in the $i$ -th species. Without epistasis ( $λ_{11} = λ_{21}$ and $λ_{22} = λ_{12}$ ), the growth rate effects are independent of the host species. Epistasis leads to the difference between $λ_{11}$ and $λ_{21}$ (or between $λ_{22}$ and $λ_{12}$ ). The population dynamics can be described by the following ODEs:

\frac{{d s}_{1}}{d t} = μ_{1} (1 + λ_{12} \frac{p_{1}}{s_{1}}) s_{1} (1 - s_{1} - γ_{2} s_{2}) - D s_{1},

\frac{{d s}_{2}}{d t} = μ_{2} (1 + λ_{21} \frac{p_{2}}{s_{2}}) s_{2} (1 - s_{2} - γ_{1} s_{1}) - D s_{2},

\frac{{d p}_{1}}{d t} = μ_{1} (1 + λ_{12}) p_{1} (1 - s_{1} - γ_{2} s_{2}) + η_{1} (s_{2} + p_{1}) (s_{1} - p_{1}) - (D + κ_{1}) p_{1},

\frac{{d p}_{2}}{d t} = μ_{2} (1 + λ_{21}) p_{2} (1 - s_{2} - γ_{1} s_{1}) + η_{2} (s_{1} + p_{2}) (s_{2} - p_{2}) - (D + κ_{2}) p_{2} .

Here, $γ_{1}$ and $γ_{2}$ are the competition strengths. $D$ is the dilution rate. $η_{1}$ and $η_{2}$ are transfer rates. $κ_{1}$ and $κ_{2}$ are the loss rates of the mobilizable genes.

We quantified the epistasis using by two ratios: $ξ_{1} = λ_{21} / λ_{11}$ and $ξ_{2} = λ_{12} / λ_{22}$ . $ξ_{1} > 0$ and $ξ_{2} > 0$ represent magnitude epistasis, while $ξ_{1} < 0$ or $ξ_{2} < 0$ represent the sign epistasis. In numerical simulations, we tested both types of epistasis. Specifically, for systems of two competing species, we approximated the areas of bistability regions under different $ξ_{1}$ and $ξ_{2}$ values, by numerical simulations with randomized parameters. Given competition strength, we randomized $μ_{1}$ and $μ_{2}$ 500 times between 0 and 1 hr^–1 following uniform distributions while keeping $μ_{1}^{0}$ and $μ_{2}^{0}$ 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 generate bistability out of the 500 random combinations.

Data availability

All the codes and source data associated with the numerical simulations and analysis of this manuscript are available on GitHub (copy archived at Hong et al., 2025).

The following data sets were generated

1. Hong J
2. Xue W
3. Wang T
(2025) Zenodo
MultistabilityByHGT.

https://doi.org/10.5281/zenodo.14928352

References

(2019) Mortality causes universal changes in microbial community composition
Nature Communications 10:2120.

https://doi.org/10.1038/s41467-019-09925-0
- PubMed
- Google Scholar
(2020) Experimental determination of evolutionary barriers to horizontal gene transfer
BMC Microbiology 20:326.

https://doi.org/10.1186/s12866-020-01983-5
- PubMed
- Google Scholar
1. Amor DR
2. Ratzke C
3. Gore J
(2020) Transient invaders can induce shifts between alternative stable states of microbial communities
Science Advances 6:eaay8676.

https://doi.org/10.1126/sciadv.aay8676
- PubMed
- Google Scholar
1. Arumugam M
2. Raes J
3. Pelletier E
4. Le Paslier D
5. Yamada T
6. Mende DR
7. Fernandes GR
8. Tap J
9. Bruls T
10. Batto J-M
11. Bertalan M
12. Borruel N
13. Casellas F
14. Fernandez L
15. Gautier L
16. Hansen T
17. Hattori M
18. Hayashi T
19. Kleerebezem M
20. Kurokawa K
21. Leclerc M
22. Levenez F
23. Manichanh C
24. Nielsen HB
25. Nielsen T
26. Pons N
27. Poulain J
28. Qin J
29. Sicheritz-Ponten T
30. Tims S
31. Torrents D
32. Ugarte E
33. Zoetendal EG
34. Wang J
35. Guarner F
36. Pedersen O
37. de Vos WM
38. Brunak S
39. Doré J
40. MetaHIT Consortium
41. Antolín M
42. Artiguenave F
43. Blottiere HM
44. Almeida M
45. Brechot C
46. Cara C
47. Chervaux C
48. Cultrone A
49. Delorme C
50. Denariaz G
51. Dervyn R
52. Foerstner KU
53. Friss C
54. van de Guchte M
55. Guedon E
56. Haimet F
57. Huber W
58. van Hylckama-Vlieg J
59. Jamet A
60. Juste C
61. Kaci G
62. Knol J
63. Lakhdari O
64. Layec S
65. Le Roux K
66. Maguin E
67. Mérieux A
68. Melo Minardi R
69. M’rini C
70. Muller J
71. Oozeer R
72. Parkhill J
73. Renault P
74. Rescigno M
75. Sanchez N
76. Sunagawa S
77. Torrejon A
78. Turner K
79. Vandemeulebrouck G
80. Varela E
81. Winogradsky Y
82. Zeller G
83. Weissenbach J
84. Ehrlich SD
85. Bork P
(2011) Enterotypes of the human gut microbiome
Nature 473:174–180.

https://doi.org/10.1038/nature09944
- PubMed
- Google Scholar
1. Beardmore RE
2. Cook E
3. Nilsson S
4. Smith AR
5. Tillmann A
6. Esquivel BD
7. Haynes K
8. Gow NAR
9. Brown AJP
10. White TC
11. Gudelj I
(2018) Drug-mediated metabolic tipping between antibiotic resistant states in a mixed-species community
Nature Ecology & Evolution 2:1312–1320.

https://doi.org/10.1038/s41559-018-0582-7
- PubMed
- Google Scholar
(2017) Extensive horizontal gene transfer in cheese-associated bacteria
eLife 6:e22144.

https://doi.org/10.7554/eLife.22144
- PubMed
- Google Scholar
1. Bordenstein SR
2. Reznikoff WS
(2005) Mobile DNA in obligate intracellular bacteria
Nature Reviews Microbiology 3:688–699.

https://doi.org/10.1038/nrmicro1233
- Google Scholar
1. Brito IL
2. Yilmaz S
3. Huang K
4. Xu L
5. Jupiter SD
6. Jenkins AP
7. Naisilisili W
8. Tamminen M
9. Smillie CS
10. Wortman JR
11. Birren BW
12. Xavier RJ
13. Blainey PC
14. Singh AK
15. Gevers D
16. Alm EJ
(2016) Mobile genes in the human microbiome are structured from global to individual scales
Nature 535:435–439.

https://doi.org/10.1038/nature18927
- PubMed
- Google Scholar
(2019) The ecology and evolution of pangenomes
Current Biology 29:R1094–R1103.

https://doi.org/10.1016/j.cub.2019.08.012
- PubMed
- Google Scholar
1. Canfora EE
2. Meex RCR
3. Venema K
4. Blaak EE
(2019) Gut microbial metabolites in obesity, NAFLD and T2DM
Nature Reviews. Endocrinology 15:261–273.

https://doi.org/10.1038/s41574-019-0156-z
- PubMed
- Google Scholar
1. Caro-Quintero A
2. Konstantinidis KT
(2015) Inter-phylum HGT has shaped the metabolism of many mesophilic and anaerobic bacteria
The ISME Journal 9:958–967.

https://doi.org/10.1038/ismej.2014.193
- PubMed
- Google Scholar
1. Chaparro-Pedraza PC
2. de Roos AM
(2020) Ecological changes with minor effect initiate evolution to delayed regime shifts
Nature Ecology & Evolution 4:412–418.

https://doi.org/10.1038/s41559-020-1110-0
- PubMed
- Google Scholar
(2018) Microbiota in the gastrointestinal tract
Medical Sciences 6:116.

https://doi.org/10.3390/medsci6040116
- PubMed
- Google Scholar
1. Ding T
2. Schloss PD
(2014) Dynamics and associations of microbial community types across the human body
Nature 509:357–360.

https://doi.org/10.1038/nature13178
- PubMed
- Google Scholar
(2019) Multistability and regime shifts in microbial communities explained by competition for essential nutrients
eLife 8:e49720.

https://doi.org/10.7554/eLife.49720
- PubMed
- Google Scholar
(2022) Functional attractors in microbial community assembly
Cell Systems 13:29–42.

https://doi.org/10.1016/j.cels.2021.09.011
- PubMed
- Google Scholar
1. Fassarella M
2. Blaak EE
3. Penders J
4. Nauta A
5. Smidt H
6. Zoetendal EG
(2021) Gut microbiome stability and resilience: elucidating the response to perturbations in order to modulate gut health
Gut 70:595–605.

https://doi.org/10.1136/gutjnl-2020-321747
- PubMed
- Google Scholar
(2005) Unsaturated fatty acids are inhibitors of bacterial conjugation
Microbiology 151:3517–3526.

https://doi.org/10.1099/mic.0.28216-0
- PubMed
- Google Scholar
1. Finlay BB
2. Medzhitov R
(2007) Host-microbe interactions: fulfilling a niche
Cell Host & Microbe 1:3–4.

https://doi.org/10.1016/j.chom.2007.02.004
- PubMed
- Google Scholar
1. Flemming HC
2. Wuertz S
(2019) Bacteria and archaea on Earth and their abundance in biofilms
Nature Reviews. Microbiology 17:247–260.

https://doi.org/10.1038/s41579-019-0158-9
- PubMed
- Google Scholar
(2015) Marine microbial community dynamics and their ecological interpretation
Nature Reviews. Microbiology 13:133–146.

https://doi.org/10.1038/nrmicro3417
- PubMed
- Google Scholar
1. Fujita H
2. Ushio M
3. Suzuki K
4. Abe MS
5. Yamamichi M
6. Iwayama K
7. Canarini A
8. Hayashi I
9. Fukushima K
10. Fukuda S
11. Kiers ET
12. Toju H
(2023) Alternative stable states, nonlinear behavior, and predictability of microbiome dynamics
Microbiome 11:63.

https://doi.org/10.1186/s40168-023-01474-5
- PubMed
- Google Scholar
(2020) Plasmid interactions can improve plasmid persistence in bacterial populations
Frontiers in Microbiology 11:2033.

https://doi.org/10.3389/fmicb.2020.02033
- PubMed
- Google Scholar
(2015) Synthetic fatty acids prevent plasmid-mediated horizontal gene transfer
mBio 6:e01032.

https://doi.org/10.1128/mBio.01032-15
- PubMed
- Google Scholar
1. Getino M
2. de la Cruz F
(2018) Natural and artificial strategies to control the conjugative transmission of plasmids
Microbiology Spectrum 6:10.

https://doi.org/10.1128/microbiolspec.mtbp-0015-2016
- PubMed
- Google Scholar
1. Gibson TE
2. Bashan A
3. Cao HT
4. Weiss ST
5. Liu YY
(2016) On the origins and control of community types in the human microbiome
PLOS Computational Biology 12:e1004688.

https://doi.org/10.1371/journal.pcbi.1004688
- PubMed
- Google Scholar
1. Gilpin ME
2. Case TJ
(1976) Multiple domains of attraction in competition communities
Nature 261:40–42.

https://doi.org/10.1038/261040a0
- PubMed
- Google Scholar
1. Gonze D
2. Lahti L
3. Raes J
4. Faust K
(2017) Multi-stability and the origin of microbial community types
The ISME Journal 11:2159–2166.

https://doi.org/10.1038/ismej.2017.60
- PubMed
- Google Scholar
1. Grilli J
2. Adorisio M
3. Suweis S
4. Barabás G
5. Banavar JR
6. Allesina S
7. Maritan A
(2017) Feasibility and coexistence of large ecological communities
Nature Communications 8:14389.

https://doi.org/10.1038/ncomms14389
- Google Scholar
1. Hartmann M
2. Six J
(2023) Soil structure and microbiome functions in agroecosystems
Nature Reviews Earth & Environment 4:4–18.

https://doi.org/10.1038/s43017-022-00366-w
- Google Scholar
Software
1. Hong J
2. Xue W
3. Wang T
(2025) MultistabilityByHGT, version 10.5281/zenodo.14928352
Zenodo.

https://doi.org/10.5281/zenodo.14928352
(2022) Plasmid co-infection: linking biological mechanisms to ecological and evolutionary dynamics
Philosophical Transactions of the Royal Society B 377:20200478.

https://doi.org/10.1098/rstb.2020.0478
- Google Scholar
(2020) Metabolic multistability and hysteresis in a model aerobe-anaerobe microbiome community
Science Advances 6:eaba0353.

https://doi.org/10.1126/sciadv.aba0353
- Google Scholar
1. Kootte RS
2. Levin E
3. Salojärvi J
4. Smits LP
5. Hartstra AV
6. Udayappan SD
7. Hermes G
8. Bouter KE
9. Koopen AM
10. Holst JJ
11. Knop FK
12. Blaak EE
13. Zhao J
14. Smidt H
15. Harms AC
16. Hankemeijer T
17. Romijn HA
18. Schaap FG
19. Olde Damink SWM
20. Ackermans MT
21. Dallinga-Thie GM
22. Zoetendal E
23. de Vos WM
24. Serlie MJ
25. Stroes ESG
26. Groen AK
27. Nieuwdorp M
(2017) Improvement of insulin sensitivity after lean donor feces in metabolic syndrome is driven by baseline intestinal microbiota composition
Cell Metabolism 26:611–619.

https://doi.org/10.1016/j.cmet.2017.09.008
- PubMed
- Google Scholar
1. Kosterlitz O
2. Muñiz Tirado A
3. Wate C
4. Elg C
5. Bozic I
6. Top EM
7. Kerr B
(2022) Estimating the transfer rates of bacterial plasmids with an adapted Luria-Delbrück fluctuation analysis
PLOS Biology 20:e3001732.

https://doi.org/10.1371/journal.pbio.3001732
- PubMed
- Google Scholar
1. Kumar S
2. Hasty J
(2023) Stability, robustness, and containment: preparing synthetic biology for real-world deployment
Current Opinion in Biotechnology 79:102880.

https://doi.org/10.1016/j.copbio.2022.102880
- PubMed
- Google Scholar
(2014) Tipping elements in the human intestinal ecosystem
Nature Communications 5:4344.

https://doi.org/10.1038/ncomms5344
- PubMed
- Google Scholar
1. Leibold MA
2. Holyoak M
3. Mouquet N
4. Amarasekare P
5. Chase JM
6. Hoopes MF
7. Holt RD
8. Shurin JB
9. Law R
10. Tilman D
11. Loreau M
12. Gonzalez A
(2004) The metacommunity concept: a framework for multi‐scale community ecology
Ecology Letters 7:601–613.

https://doi.org/10.1111/j.1461-0248.2004.00608.x
- Google Scholar
1. Lipsitch M
2. Colijn C
3. Cohen T
4. Hanage WP
5. Fraser C
(2009) No coexistence for free: neutral null models for multistrain pathogens
Epidemics 1:2–13.

https://doi.org/10.1016/j.epidem.2008.07.001
- PubMed
- Google Scholar
1. Liu L
2. Chen X
3. Skogerbø G
4. Zhang P
5. Chen R
6. He S
7. Huang DW
(2012) The human microbiome: a hot spot of microbial horizontal gene transfer
Genomics 100:265–270.

https://doi.org/10.1016/j.ygeno.2012.07.012
- PubMed
- Google Scholar
1. Lopatkin AJ
2. Meredith HR
3. Srimani JK
4. Pfeiffer C
5. Durrett R
6. You L
(2017) Persistence and reversal of plasmid-mediated antibiotic resistance
Nature Communications 8:1689.

https://doi.org/10.1038/s41467-017-01532-1
- PubMed
- Google Scholar
1. Mangal U
2. Noh K
3. Lee S
4. Cha JK
5. Song JS
6. Cha JY
7. Lee KJ
8. Kim KM
9. Kwon JS
10. Choi SH
(2023) Multistability and hysteresis in states of oral microbiota: Is it impacting the dental clinical cohort studies?
Journal of Periodontal Research 58:381–391.

https://doi.org/10.1111/jre.13098
- Google Scholar
(2018) Microbiomes as metacommunities: understanding host-associated microbes through metacommunity ecology
Trends in Ecology & Evolution 33:926–935.

https://doi.org/10.1016/j.tree.2018.09.002
- Google Scholar
1. Moustafa A
2. Li W
3. Anderson EL
4. Wong EHM
5. Dulai PS
6. Sandborn WJ
7. Biggs W
8. Yooseph S
9. Jones MB
10. Venter JC
11. Nelson KE
12. Chang JT
13. Telenti A
14. Boland BS
(2018) Genetic risk, dysbiosis, and treatment stratification using host genome and gut microbiome in inflammatory bowel disease
Clinical and Translational Gastroenterology 9:e132.

https://doi.org/10.1038/ctg.2017.58
- PubMed
- Google Scholar
(2021) Conjugation inhibitors effectively prevent plasmid transmission in natural environments
mBio 12:e0127721.

https://doi.org/10.1128/mBio.01277-21
- PubMed
- Google Scholar
(2011) What traits are carried on mobile genetic elements, and why?
Heredity 106:1–10.

https://doi.org/10.1038/hdy.2010.24
- PubMed
- Google Scholar
1. Ravel J
2. Gajer P
3. Abdo Z
4. Schneider GM
5. Koenig SSK
6. McCulle SL
7. Karlebach S
8. Gorle R
9. Russell J
10. Tacket CO
11. Brotman RM
12. Davis CC
13. Ault K
14. Peralta L
15. Forney LJ
(2011) Vaginal microbiome of reproductive-age women
PNAS 108:4680–4687.

https://doi.org/10.1073/pnas.1002611107
- Google Scholar
1. Relman DA
(2012) The human microbiome: ecosystem resilience and health
Nutrition Reviews 70 Suppl 1:S2–S9.

https://doi.org/10.1111/j.1753-4887.2012.00489.x
- PubMed
- Google Scholar
(2016) Type IV traffic ATPase TrwD as molecular target to inhibit bacterial conjugation
Molecular Microbiology 100:912–921.

https://doi.org/10.1111/mmi.13359
- PubMed
- Google Scholar
(2014) Positive epistasis between co-infecting plasmids promotes plasmid survival in bacterial populations
The ISME Journal 8:601–612.

https://doi.org/10.1038/ismej.2013.182
- PubMed
- Google Scholar
(2018) Coexistence of many species in random ecosystems
Nature Ecology & Evolution 2:1237–1242.

https://doi.org/10.1038/s41559-018-0603-6
- PubMed
- Google Scholar
1. Silva RF
2. Mendonça SCM
3. Carvalho LM
4. Reis AM
5. Gordo I
6. Trindade S
7. Dionisio F
(2011) Pervasive sign epistasis between conjugative plasmids and drug-resistance chromosomal mutations
PLOS Genetics 7:e1002181.

https://doi.org/10.1371/journal.pgen.1002181
- PubMed
- Google Scholar
1. Smillie CS
2. Smith MB
3. Friedman J
4. Cordero OX
5. David LA
6. Alm EJ
(2011) Ecology drives a global network of gene exchange connecting the human microbiome
Nature 480:241–244.

https://doi.org/10.1038/nature10571
- PubMed
- Google Scholar
(2015) Horizontal gene transfer: building the web of life
Nature Reviews. Genetics 16:472–482.

https://doi.org/10.1038/nrg3962
- PubMed
- Google Scholar
1. Sun L
2. Wang S
3. Narsing Rao MP
4. Shi Y
5. Lian Z-H
6. Jin P-J
7. Wang W
8. Li Y-M
9. Wang K-K
10. Banerjee A
11. Cui X-Y
12. Wei D
(2023) The shift of soil microbial community induced by cropping sequence affect soil properties and crop yield
Frontiers in Microbiology 14:1095688.

https://doi.org/10.3389/fmicb.2023.1095688
- Google Scholar
1. van den Berg NI
2. Machado D
3. Santos S
4. Rocha I
5. Chacón J
6. Harcombe W
7. Mitri S
8. Patil KR
(2022) Ecological modelling approaches for predicting emergent properties in microbial communities
Nature Ecology & Evolution 6:855–865.

https://doi.org/10.1038/s41559-022-01746-7
- PubMed
- Google Scholar
1. Werren JH
(2011) Selfish genetic elements, genetic conflict, and evolutionary innovation
PNAS 108:10863–10870.

https://doi.org/10.1073/pnas.1102343108
- Google Scholar
1. Zaoli S
2. Grilli J
(2021) A macroecological description of alternative stable states reproduces intra- and inter-host variability of gut microbiome
Science Advances 7:eabj2882.

https://doi.org/10.1126/sciadv.abj2882
- PubMed
- Google Scholar
1. Zhu S
2. Hong J
3. Wang T
(2024) Horizontal gene transfer is predicted to overcome the diversity limit of competing microbial species
Nature Communications 15:800.

https://doi.org/10.1038/s41467-024-45154-w
- Google Scholar

Article and author information

Author details

Juken Hong

Key Laboratory of Quantitative Synthetic Biology, Shenzhen Institute of Synthetic Biology, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen, China

Contribution
Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft

Competing interests
No competing interests declared
Wenzhi Xue

Key Laboratory of Quantitative Synthetic Biology, Shenzhen Institute of Synthetic Biology, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen, China

Contribution
Investigation, Visualization, Writing – original draft

Competing interests
No competing interests declared
Teng Wang

Key Laboratory of Quantitative Synthetic Biology, Shenzhen Institute of Synthetic Biology, Shenzhen Institutes of Advanced Technology, Chinese Academy of Sciences, Shenzhen, China

Contribution
Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing

For correspondence
t.wang1@siat.ac.cn

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0009-0000-5176-5095

Funding

National Natural Science Foundation of China (12401660)

Teng Wang

Shenzhen Institue of Synthetic Biology Scientific Research Program (HSE499011086)

Teng Wang

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This study was partially supported the National Natural Science Foundation of China (12401660 to TW) and Shenzhen Institute of Synthetic Biology Scientific Research Program (HSE499011086 to TW).

Version history

Sent for peer review: June 18, 2024
Preprint posted: June 22, 2024
Reviewed Preprint version 1: November 21, 2024
Reviewed Preprint version 2: February 4, 2025
Version of Record published: March 3, 2025

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.99593. This DOI represents all versions, and will always resolve to the latest one.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

200

views
12

downloads
0

citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Juken Hong
Wenzhi Xue
Teng Wang

(2025)

Emergence of alternative stable states in microbial communities undergoing horizontal gene transfer

eLife 13:RP99593.

https://doi.org/10.7554/eLife.99593.3

Share this article

Cite this article

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

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

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

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

HGT creates chances of multistability for communities under strong environmental selection.

Author details

Juken Hong

Contribution

Competing interests

Wenzhi Xue

Contribution

Competing interests

Teng Wang

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags

Further reading