The Persistence Potential of Mobile Genetic Elements

Mobile genetic elements (MGEs), such as plasmids, phages, and transposons, play a critical role in mediating the transfer and maintenance of diverse traits and functions in microbial communities. This role depends on the ability of MGEs to persist. For a community consisting of multiple populations transferring multiple MGEs, however, the conditions underlying the persistence of these MGEs are poorly understood. Computationally, this difficulty arises from the combinatorial explosion associated with describing the gene flow in a complex community using the conventional modeling framework. Here, we describe an MGE-centric framework that makes it computationally feasible to analyze such transfer dynamics. Using this framework, we derive the persistence potential: a general, heuristic metric that predicts the persistence and abundance of any MGEs. We validate the metric with engineered microbial consortia transferring mobilizable plasmids and quantitative data available in the literature. Our modeling framework and the resulting metric have implications for developing a quantitative understanding of natural microbial communities and guiding the engineering of microbial consortia.


Introduction
Mobile genetic elements (MGEs), including plasmids, transposons and phages are major components of the metagenome of microbes 1,2 . Hundreds of plasmids have been identified in diverse microbial communities from rat cecum 3 , cow rumen 4 , sludge 5 , and marine water 6 . In the reference genomes of the human gut microbiome, ~16,000 genes are identified as mobile 7 . MGEs encode diverse biological functions including metabolic capabilities 8 , pathogenic virulence 9 , or traits to cope with environmental stresses [10][11][12] . The interplay between MGEs and the core genomes of the host cells shapes the evolution of microbial communities 13 .
The ability for an MGE to persist in a microbial community can influence the dynamics, function, and even survival of the community 14,15 . Promotion or suppression of MGE persistence, depending on the context, has applications in medicine, biosafety, and biotechnology. For instance, eliminating an MGE encoding resistance to an antibiotic will sensitize a community to the antibiotic 16 , and this resistance reversal can enable more effective use of antibiotics 17 . There is always a biosafety concern about the risk associated with the spread of synthetic genetic constructs into the environment 18,19 . To reduce the risk, the use of genetic vectors with restricted capability to be maintained outside of the laboratory has been proposed 18,19 . In biotechnology, plasmid instability is a major impediment to the large-scale production of recombinant protein products 20,21 . The expression of recombinant proteins encoded by the plasmids are usually burdensome to the cell metabolism, and, as a consequence, the plasmid-carrying cells can be outcompeted by the faster-growing plasmid-free populations. A strategy to promote plasmid persistence could overcome this limitation.
The quantitative studies of MGE persistence and abundance in microbial communities is challenging due to the lack of an effective theoretical framework 22,23 . Since the 1970s, populationbiology models have been developed to predict the persistence of a single MGE in a single species 16,[24][25][26][27][28] . However, microbes in nature often live in complex communities consisting of diverse species and mobile elements 7,29,30 . The limited scope of past quantitative analyses is due in part to the challenge associated with modeling complex communities using the conventional modeling framework, which we refer to as 'subpopulation-centric framework' (SCF). In SCF, a population carrying a particular combination of MGEs is considered a unique subpopulation that requires one ordinary differential equation (ODE) to describe. Modelling a community containing two species and two mobile elements requires eight ODEs, each describing one subpopulation. The model complexity increases combinatorially with the number of MGEs ( Fig. 1A and B). For example, the marine microbiome in a bottle of sea water is estimated to contain ~160 species 30 and ~180 plasmids 6 . In SCF, ~2.4 × 10 '( ODEs and ~2.7 × 10 **+ parameters are needed to model the gene flow dynamics of this community. Therefore, it is currently impossible to model a system in which many MGEs are involved.

An MGE-centric framework to model gene flow
To overcome the challenge associated with SCF, we developed an MGE-centric framework (MCF) that focuses on the overall abundance of each MGE in the community by accounting for the average fitness effect of MGEs to the host. To illustrate the key concepts, consider the transfer of n types of mobilizable plasmids in a community consisting of m species. Let -( = 1,2, … , ) be the abundance of i-th species and -4 ( = 1,2, … , ) the abundance of the i-th species carrying the j-th plasmid (Fig. S1A). The dynamics can be approximately described by two groups of ODEs (see Supplementary  Eq.1 describes the collective growth and dilution of -, where -; is the effective growth rate of the species, is the dilution rate, andrepresents the combined growth effect of all the plasmids carried by the species. In Eq.2, the first term describes the growth of -4 . -4 ; is the effective growth rate; it differs from -; due to the growth effect of the j-th plasmid. -4 is the combined growth effect of the other plasmids. In the second term, 4F-is the conjugation efficiency when the plasmid is transferred from the k-th species to the i-th species. -− -4 is the total abundance of subpopulations ofnot carrying the j-th plasmid. The third term describes the plasmid loss due to segregation error (at a rate constant of -4 ) and dilution.
In this formulation, -4 / -indicates the relative abundance of the j-th plasmid in the i-th species; is the relative abundance of the j-th plasmid in the entire community. Because of the formulation of -, -4 , -; and -4 ; , this framework is applicable for describing arbitrary microbial communities transferring multiple plasmids. In particular, the interactions between populations are accounted for by the appropriate formulation of -; and -4 ; ; plasmid incompatibility can be accounted for by adapting the formulations of -4 and the conjugation terms (Fig. S1C 1A and B, Table S2, see Supplementary Information section 2.2). This simplification renders our framework capable of predicting the dynamics of species composition as well as the distribution patterns of each plasmid (Fig. 1C). In particular, we modeled a community consisting of 200 species transferring 200 plasmids. The conventional SCF requires about 3.2 × 10 (N ODEs and 5.2 × 10 *N( parameters, which is impractical to construct. In contrast, this community can be modeled by 4.0 × 10 + ODEs and 8.1 × 10 ( parameters using MCF, which is feasible to both construct and calculate (Fig. S1B). Although acquiring all the parameters in such complex communities remains challenging due to the technical limitations 31 , MCF enables the theoretical analysis of the dynamics and persistence conditions of MGEs.
This drastic reduction in model complexity is made possible by lumping multiple distinct subpopulations into an average one and then accounting for the average fitness effect of each plasmid. In particular, -lumps all subpopulations carrying different combinations of plasmids or no plasmids and -4 lumps all subpopulations carrying the j-th plasmid. If the plasmids do not have any growth effects, the lumping is exact; the MCF is equivalent to SCF. In general, however, plasmids can confer burden or benefit, which will cause deviation between the results computed from these two frameworks. To evaluate this discrepancy, we conducted numerical simulations on communities transferring one and two plasmids. Testing higher number of plasmids in SCF is computationally prohibitive due to the combinatorial explosion. Indeed, simulation results suggest that the fitness effects of the plasmids are the main factors that determine the discrepancy. Within a reasonably wide range of fitness costs, the predictions of these two models match well, and the smaller the fitness effects, the smaller the discrepancy. (Fig. S2, see Supplementary Information section 2.3). In general, we note that both frameworks represent approximations of real biological systems. The ultimate test of each framework should be from experiments.
The same framework is applicable to the analysis of other MGEs, including phages and transposons. Temperate phages exit the host cells by budding or secretion 32 , without killing the host. This transfer process can be described by MCF by interpreting conjugation efficiencies ( 's) as phage transfer rates. Lytic phages, which cause host death and lysis during release 27,33,34

The persistence potential of an MGE
Our framework makes it feasible to develop a criterion of MGE persistence for complex communities. We first considered an idealized community, where each species has the same set of kinetic parameters and equal abundance ( Fig. 2A). Because of its symmetry, this community allows the analytical derivation of a metric ( ) that determines the persistence and steady-state abundance of the plasmids relative to the total cell density ( Fig. 2A, also see Supplementary Information section 2.5) , [3].
We term this metric the persistence potential for the MGE. In this idealized community, a burdensome MGE ( >0) persists if and only if > 1. If the mobile element is beneficial ( <0), which leads to <0, the element will also persist.
In general, the constituent species in a community are not symmetric (Fig. 2B). Their abundances are different from each other; they have different growth rates; and they transfer plasmids at different rates. For such communities, deriving an analytical solution is not possible. We thus took a heuristic approach to generalize the idealized metric. We kept its derived formulation but replaced each parameter with the weighted average of the corresponding parameter in the general community accounting for the relative abundances of different is the total abundance of all the populations: Z = ∑ -I -J* (Fig. 2). Correspondingly, we replace the number of species, m, with the community diversity , which is calculated through . The general form of persistence potential becomes: [4] where ̅ , ̅ , ̅ , ̅ and stand for the weighted averages of species growth rate, plasmid loss rate, fitness cost, horizontal transfer rate, and dilution rate, respectively ( Fig. 3A and B).
To evaluate the predictive power of the generalized persistence potential, we performed numerical simulations with 2,000 sets of randomized parameters in communities consisting of 5-100 species and 1-50 plasmids. Our randomization process produced communities with various settings that mimic the diversity of the natural genetic-exchange communities. To ensure there were adequate numbers of species coexisting at steady states, the community was divided into multiple coexisting niches. Each niche had its own carrying capacity, and different species within the same niche competed with each other (Fig. S1D). The dynamics of each community was simulated for up to ~30,000 hours until it became unchanged, and we treated it as 'steady state'. The fractions of the MGE-carrying cells were then calculated with respect to values. Our simulation results suggest that the general metric, despite its heuristic nature, remains a robust predictor on whether and to what extent an MGE can persist, with a transition at = 1 (Fig. 3C). When 0 < <1, the abundance of the MGE is close to 0; when > 1, the abundance of the MGE increases monotonically with . The data points computed from the randomized parameter sets approximately collapsed into a single curve, underscoring the general predictive power of .
In addition to assuming symmetry, the basic form of (Eq. 4) is derived by assuming that the system has reached steady state. However, this assumption is not required for the approximate predictive power of . In particular, in the numerical simulation described above, the values have similar predictive power for the MGE abundance well before the system has reached steady state (Fig. S3). This result underscores the general predictive power of and its applicability to experimental systems, which may not be at steady state.

Experimental validation of the persistence potential
To test the predictive power of the persistence potential, we engineered eight communities transferring mobilizable plasmids. Community 1 through 7 were constructed from three E. coli strains (denoted X, B, and R). Strain B expresses BFP on the chromosome, R expresses dTomato, and strain X is not fluorescent. The mobilizable plasmid K was transferred among the strains and expresses GFP constitutively. Using flow cytometry, this system allows the simultaneous quantification of plasmid abundance (with GFP) and population compositions (with BFP and dTomato) (Fig. S4A). Community 8 contained two E. coli strains (MG1655 and DH5α) and five conjugative plasmids (F', PCU1, R388, R6K, and RP4). The community composition and plasmid abundance were quantified by selective plating.
After plasmids were introduced into the communities, dilutions of four different ratios ( 10 b × , 10 + ×, 10 ' × , 10 ( × ) were performed every 24 hours to maintain the growth. We monitored the population dynamics daily over the following 15 days and obtained the fractions of each population and the plasmid (Fig. S4B-H, Fig. S5C). We measured the conjugation rates and fitness costs ( Fig. S4I and J, Fig. S5A and B) and also estimated the average dilution rate and plasmid segregation loss rate. With these parameters, we were able to calculate the plasmid persistence potential in each community. The results were well matched to the predicted pattern, suggesting that values determine plasmid abundance in microbial communities (Fig. 4A).
To examine the general applicability of the metric, we reanalyzed the data from 13 previous studies that have provided sufficient measurements on the kinetic parameters and abundance of plasmids 16,[36][37][38][39][40][41][42][43] , transposons 35 and bacteriophages 33,34,44 . The microbial communities analyzed in these studies covered one, two, or three populations, and up to three MGEs. We collected a total of 92 data points (see Supplementary Information section 3) covering a range from 1.46 × 10 O' to 2.62 × 10 ' (Fig. 4B). These data confirm the predictive power of : in general, the MGE persists when > 1, and its relative abundance increases with . In contrast, when < 1, the MGE tend not to persist, with its relative abundance close to 0.

Discussion
Our work addresses two fundamental challenges facing the quantitative analysis of gene flow dynamics in microbial communities. First, it has been impractical to simulate a complex community in which more than dozen MGEs are transferred. This challenge is resolved by using our MGE-central modeling framework through drastic dimension reduction in model formulation. Second, the MCF enables the heuristic derivation of a metric that predicts the persistence and abundance of any MGEs in a microbial community based on its kinetic parameters and the community composition. Both the modeling framework and the derived MGE persistence potential have implications for future efforts to understand, control, and exploit horizontal gene transfer dynamics in microbial communities.
For instance, the persistence potential predicts that the abundance of an MGE in a microbial community is sensitive to the average growth rate of the constituent populations, especially when the growth rate is close to the system dilution rate. Consistent with this notion, nutrient enrichments, which in general increase community average growth rate, have been shown to increase the relative abundance of plasmids 45 and bacteriophages 44 . In the inflamed mouse gut, the growth of commensal Enterobacteriaceae such as E.coli is boosted by the nitrate generated as a by-product of the host inflammatory response 46 ; this transient boom of Enterobacteriaceae was shown to correlate with an increase in transconjugant abundance 47,48 . In both situations, an increase in the growth rate can promote persistence by increasing the persistence potential directly (through a decrease of its denominator term, Eq2) or indirectly (through increasing the conjugation efficiency) 49 . Conversely, our metric suggests that modulating the overall growth rate of gut microbiota by nutrient supply or its dilution by water inflow could be effective strategies to regulate the abundance of the mobile gene pool in the human gut 50,51 .
Our metric also predicts the influence of bacterial community composition on the MGE persistence potential. For example, the presence of bacterial species with strong efficiencies for interspecific MGE transfers is predicted to promote the average persistence potential of the MGEs. Indeed, for bacteria with poor conjugation efficiencies, coculturing them with efficient donors has been shown to enhance plasmid transfer and maintenance 36,52 . The human gut is equipped with an enormous diversity of microbes 53 , and multiple factors such as diet 54 , age 55 , and antibiotic administration 56 can alter the composition of its microbiome. Our metric predicts that these factors potentially affect the ability of the gut to sustain MGEs and change the contents of mobile elements in gut microbiome. This notion provides new insights on reversing antibiotic resistance in complex microbial communities 57 .
Our results also have implications for the engineering of microbial consortia. Engineering complex bacterial communities in useful ways remains challenging due to the lack of understanding of the ecological principles and intercellular metabolic interactions 58 . Mobile genetic elements are potentially powerful tools for function-oriented microbiota engineering 59 . Our modeling framework and the MGE persistence potential can guide such efforts. In particular, the persistence potential reveals how a few key kinetic parameters predict the approximate abundance of the MGE.  Here, the community is fully symmetric in terms of parameters and population sizes. Each colored circle represents a constituent species; the area of the circle is proportional to the species size. Each directed arrow represents the transfer of an MGE, with the arrow thickness indicating the transfer rate. In this case, can be analytically derived.
(B) Generalization of . We generalized the formulation of by replacing each parameter with its weighted average, which factors in the heterogeneity of species size. The weighted averages of , , and were calculated by

Strains and plasmids
The compositions of the eight engineered communities are shown in Table S1. E.coli strain MG1655 without fluorescence markers was denoted as strain X 16 . E.coli strain DA26735 with chromosomal BFP and chloramphenicol resistance (Cm R ) was denoted as strain B 16,60 , and E.coli strain DA32838 with chromosomal dTomato and Cm R was denoted as strain R 16,60 . All three strains carry helper F plasmid FHR that expresses tetracycline resistance (Tet R ) 16 . FHR is not transmissible but encodes the conjugation machinery that mobilize plasmid K. Plasmid K expresses GFP under the control of a strong constitutive PR promoter, and expresses kanamycin resistance (Kan R ) 16,61 . Plasmid K also carries oriT, so it can be transferred through conjugation 16,60 .
The multi-plasmid community was composed of E.coli strain MG1655 and DH5α. These two strains were distinguished from each other via blue-white screening on X-gal plates, where MG1655 and DH5α colonies were blue and white, respectively. These communities transferred five conjugative plasmids: F' (lncF, Tet R ), PCU1 (lncN, Amp R ),R388 (lncW, Tm R ), R6K (lncX, Strp R ), and RP4 (lncP, Kan R ). These five plasmids are compatible with each other and carry different antibiotic resistance markers. The plasmids were distinguished from each other via selective plating.

Long-term dynamics of the engineered microbial communities
Our methods for measuring the long-term plasmid dynamics are based on the protocols established by Lopatkin et al. 16 . Single colonies of three strains (X, B and R) carrying plasmid FHR and K were grown overnight at 37℃ for 16 hours with shaking (250 rpm) in LB culture (LB broth from APEX) containing appropriate antibiotics (100 µg/mL Cm, 50 µg/mL Kan, or 20 µg/mL Tet). The overnight cultures were resuspended in M9 medium (M9CA medium broth powder from Amresco, supplemented with 0.1 mg/mL thiamine, 2 mM MgSO4, 0.1 mM CaCl2, and 0.4% w/v glucose) without antibiotics, and diluted to the initial density of 10 6 cells/mL. We constructed seven communities using the combinations of these three strains: (a) X; (b) B; (c) R; (d) X+B; (e) X+R; (f) B+R; (g) X+R+B. For the communities (d)-(g), the members of each community were mixed in the equal ratio, and the mixtures were diluted to the initial density of 10 6 cells/mL. The cells were then distributed in a 96-well plate to a final volume of 200 µL/well , and each community had 12 replicates. The 96-well plate was covered with an AeraSeal TM film sealant (Sigma-Aldrich, SKU A9224) followed by a Breath-Easy sealing membrane (Sigma-Aldrich, SKU Z380059). Plates were shaken at 250 rpm at 37℃ for 23 hours. This was denoted as day 0. On day 1, the 12 replicates of each community were divided into four groups, each with three replicates. From each well, 2 µL was removed for flow cytometry. The four groups were subjected to four dilution ratios (10 b ×, 10 + ×, 10 ' ×, 10 ( ×). The new plates were sealed using both membranes and placed back into the incubator. The same protocols were repeated daily from day 1 to day 15.
For the community transferring multiple plasmids, we first transformed the five plasmids into MG1655, respectively, and obtained five plasmid-carrying strains (M F' , M PCU1 , M R388 , M R6K , and M RP4 ). Single colonies of the six strains (M F' , M PCU1 , M R388 , M R6K , M RP4 , and DH5α) were grown overnight at 37℃ for 16 hours with shaking (250 rpm) in LB culture containing appropriate antibiotics (20 µg/mL Tet, 100 µg/mL Amp, 10 µg/mL Tm, 50 µg/mL Strp, and 50 µg/mL Kan, respectively). The overnight cultures were resuspended in M9 medium without antibiotics, and diluted to the initial density of 10 6 cells/mL. The cells of these six strains were mixed in a ratio of 1:1:1:1:1:5 (M F' : M PCU1 : M R388 : M R6K : M RP4 : DH5α), and the mixtures were diluted to the initial density of 10 6 cells/mL. The cells were then distributed in a 96-well plate to a final volume of 200 µL/well, with 12 replicates. The 96-well plate was covered with an AeraSeal TM film sealant followed by a Breath-Easy sealing membrane. Plates were shaken at 250 rpm at 37℃ for 23 hours. This was denoted as day 0. At day 1, the 12 replicates of each community were divided into four groups, each with three replicates. From each well, 10 µL was removed for selective plating. The four groups were subjected to four dilution ratios (10 b ×, 10 + ×, 10 ' ×, 10 ( ×), respectively. The new plates were sealed using both membranes and placed back to the incubator. The same protocols were repeated daily from day 1 to day 15. The ratio between MG1655 and DH5α cells was determined by plating on the X-gal plates (100 µg/mL X-gal, 1mM IPTG). The relative abundance of each plasmid was determined by plating on the plates with the corresponding antibiotics (20 µg/mL Tet, 100 µg/mL Amp, 10 µg/mL Tm, 50 µg/mL Strp, and 50 µg/mL Kan).

Flow cytometry
The community composition and plasmid abundance were quantified using a flow cytometer (MACSQuant ® VYB Analyzer). From day 1 to day 15, the overnight cultures were resuspended and diluted to 1: 1000 in 200 µL fresh M9 media before running through the flow. The channels of emission detectors were set as V1 (450/50 nm) for BFP, Y2 (615/20 nm) for dTomato, and B1 (525/50 nm) for GFP. For each sample, 10,000 cells were collected. All data analysis was performed using FlowJo 10.5.3.

Measuring the fitness costs and conjugation efficiencies of the mobilizable plasmids
Singles colonies of cells that carry the plasmid K (denoted X K , B K , and R K ) or do not carry the plasmid K (denoted X 0 , B 0 and R 0 ) were grown overnight at 37℃ for 16 hours with shaking in LB with appropriate antibiotics. The cultures were resuspended in M9 medium without antibiotics and then diluted in 1: 10 + ratio. The growth curves of these six cell types were measured using a plate reader (TECAN infinite M200 PRO). Six replicates per cell type were used for quantification. The growth rate constants were calculated as the effective growth rates in exponential phases (e.g., the first 5 hours). First, we smoothed the growth curves to filter out the random noises from the plate reader. We then plotted the increments / with regards to the cell density . Linear regression was performed, and the slope was obtained as the growth rate constants. The fitness cost was determined by normalizing the growth rates of plasmid-free populations X 0 , B 0 , and R 0 by the growth rates of plasmid-carrying populations X K , B K , and R K , respectively. The values were obtained through = − 1. The fitness costs of the five conjugative plasmids (F', PCU1, R388, R6K, and RP4) in the host strain MG1655 were measured in the same way as described.
Conjugation efficiencies were estimated using the protocols established by Lopatkin et al. with modifications 16 . X K , B K , and R K served as the donors, and X 0 , B 0 , and R 0 served as the transconjugants. To distinguish donors with transconjugants, we transformed the recipients with another non-mobilizable plasmid pJM31, which carried ColE1 origin and expressed ampicillin resistance (Amp R ). Therefore, donors, recipients, and transconjugants can be distinguished by different resistance markers: Kan for donors, Amp for recipients, and Kan+Amp for transconjugants. Overnight cultures of donors and recipients in LB media with appropriate selection agents were resuspended and diluted (1:100) in fresh LB media. Cells were incubated at 37℃ with shaking for 2 hours until they reached exponential phase. The cells were then resuspended in M9 media and mixed in 1:1 ratio with a total volume of 200 µL. Mixtures were incubated at room temperature (25℃) for 1 hour without shaking. The donor, recipient, and transconjugant densities were measured by diluting the mixtures (1:10 6 for donor and recipient, 1:10 4 for transconjugant), and spreading three replicates onto corresponding selective plates. The conjugation efficiency was obtained as = Z s•t•∆v , where T, D, R stand for the cell densities of transconjugant, donor and recipient, respectively. Before being plugged into the model, the measured values of need to be normalized with respect to the maximum carrying capacity I . We estimated I to be 6 × 10 w cells/mL, which corresponds to OD600≈ 1.2. The dilution rates D were calculated from the daily dilution ratio by D = ]M{ N+ | 62 . Therefore, dilution ratios of 10 b ×, 10 + ×, 10 ' ×, 10 ( × are equal to dilution rates of 0.2878 ℎ O* , 0.3838 ℎ O* , 0.4797 ℎ O* , and 0.5756 ℎ O* , respectively. The loss rate of plasmid K is very small compared with dilution rates D. We used the value ( = 0.001 ℎ O* ) measured by Lopatkin et al. 16 .
The conjugation efficiency of the five conjugative plasmids were measured using similar protocols. Since the abundance of DH5α became negligible in the cocultures at the end of the longterm experiments, only the conjugation between MG1655 cells were considered. For plasmid F' (Tet R ), PCU1 (Amp R ),R388 (Tm R ), and R6K (Strp R ), MG1655 transformed with the nonmobilizable plasmid K-(Kan R , same as the plasmid K but without the transfer origin) served as the recipient. For plasmid RP4 (Kan R ), MG1655 transformed with the non-mobilizable plasmid pJM31 (Amp R ) served as the recipient.