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

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

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

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

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

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

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

HGT promotes the multistability of communities with multiple niches.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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