# Plasmid Interactions Can Improve Plasmid Persistence in Bacterial Populations

^{1}Department of Pharmacy, Faculty of Health Sciences, UiT The Arctic University of Norway, Tromsø, Norway^{2}Department of Plant Biology, Faculty of Sciences, University of Lisbon, Lisbon, Portugal^{3}cE3c – Centre for Ecology, Evolution and Environmental Changes, Faculty of Sciences, University of Lisbon, Lisbon, Portugal

It is difficult to understand plasmid maintenance in the absence of selection and theoretical models predict the conditions for plasmid persistence to be limited. Plasmid-associated fitness costs decrease bacterial competitivity, while imperfect partition allows the emergence of plasmid-free cells during cell division. Although plasmid conjugative transfer allows mobility into plasmid-free cells, the rate of such events is generally not high enough to ensure plasmid persistence. Experimental data suggest several factors that may expand the conditions favorable for plasmid maintenance, such as compensatory mutations and accessory genes that allow positive selection. Most of the previous studies focus on bacteria that carry a single plasmid. However, there is increasing evidence that multiple plasmids inhabit the same bacterial population and that interactions between them affect their transmission and persistence. Here, we adapt previous mathematical models to include multiple plasmids and perform computer simulations to study how interactions among them affect plasmid maintenance. We tested the contribution of different plasmid interaction parameters that impact three biological features: host fitness, conjugative transfer and plasmid loss – which affect plasmid persistence. The interaction affecting conjugation was studied in the contexts of intracellular and intercellular interactions, i.e., the plasmids interact when present in the same cell or when in different cells, respectively. First, we tested the effect of each type of interaction alone and concluded that only interactions affecting fitness (epistasis) prevented plasmid extinction. Although not allowing plasmid maintenance, intracellular interactions increasing conjugative efficiencies had a more determinant impact in delaying extinction than the remaining parameters. Then, we allowed multiple interactions between plasmids and concluded that, in a few cases, a combined effect of (intracellular) interactions increasing conjugation and fitness lead to plasmid maintenance. Our results show a hierarchy among these interaction parameters. Those affecting fitness favor plasmid persistence more than those affecting conjugative transfer and lastly plasmid loss. These results suggest that interactions between different plasmids can favor their persistence in bacterial communities.

## Introduction

Mobile plasmids carry genes essential for their replication, conjugative transfer, and stability in the host population (Toussaint and Merlin, 2002). Additionally, they may encode traits useful for their hosts, such as resistance to antibiotics and heavy metals (Foster, 1983), virulence factors (Gyles and Boerlin, 2014) or catabolism of xenobiotics (Nojiri et al., 2004). This genetic cargo promotes bacterial diversification (Heuer and Smalla, 2012) reshaping lifestyles (Kaper et al., 2004; Howard et al., 2009; Hadi, 2020) that can become critical to human populations, as, for instance, the antibiotic resistance crisis (San Millan, 2018).

These traits are beneficial only under specific circumstances and cannot explain plasmid persistence (Carroll and Wong, 2018). Upon cell division, imperfect plasmid segregation allows plasmid-free cells to emerge in the population, and plasmid loss tends to provide a growth advantage in a non-selective environment by removing the plasmid-imposed fitness cost. Moreover, gene migration into the chromosome renders the plasmid obsolete, because plasmid loss would no longer be detrimental to the host even under conditions selecting for beneficial genes since the beneficial gene is kept (Bergstrom et al., 2000).

Although plasmids can enforce their persistence by conjugating back into plasmid-free cells, these rates of transfer seem to be too low to allow long-term plasmid maintenance. The low transfer rates result, at least in part, from plasmid-encoded mechanisms that control the expression of the genes involved in the conjugative process (Frost and Koraimann, 2010). By repressing these genes, for example when there are no plasmid-free cells nearby (Koraimann and Wagner, 2014), the plasmid reduces its conjugation and consequently imposes lower fitness costs on the host (Haft et al., 2009). Furthermore, repression of the conjugative machinery may prevent infection by bacteriophages that use conjugative pili as receptors (Anderson, 1968, but see an alternative explanation in Dionisio, 2005). However, in newly formed transconjugants there is a time-delay until the repression of conjugative genes is re-established. Therefore, plasmids become transitorily derepressed in a community due to these cells, consequently displaying much higher conjugative efficiencies. This effect of transitory derepression that leads to epidemic spread could, in theory, explain plasmid maintenance (Lundquist and Levin, 1986), but the equilibrium density of donor and recipient populations would have to be unrealistically high (Simonsen, 1991).

Mutational events can also increase the odds of plasmid persistence, by reducing the plasmid fitness cost (Dahlberg and Chao, 2003; San Millan et al., 2014b; Harrison et al., 2015; Porse et al., 2016; Stalder et al., 2017) such that sometimes plasmids even become beneficial (Bouma and Lenski, 1988; Dionisio et al., 2005; Starikova et al., 2013; Loftie-Eaton et al., 2017), or improve its stability in the host (Sota et al., 2010; Loftie-Eaton et al., 2016; Yano et al., 2016). Amelioration can occur in a diversity of environments ranging from antagonistic to mutualistic (Harrison et al., 2015), and pulses of positive selection (San Millan et al., 2014b; Stevenson et al., 2018) facilitates plasmid maintenance. Likewise, distinct plasmids affect the host differently (Hall et al., 2015; Di Luca et al., 2017; San Millan et al., 2018) and vice versa, meaning that the same plasmid behaves differently in response to different hosts, varying its fitness/physiological effect (Harr and Schlötterer, 2006; Kottara et al., 2018), conjugative efficiency (Dionisio et al., 2002; De Gelder et al., 2005; Sheppard et al., 2020) and stability (De Gelder et al., 2007). Therefore, plasmid behavior may also differ in populations composed of multiple hosts and/or plasmids (De Gelder et al., 2008; Hall et al., 2016; Kottara et al., 2016; Harrison et al., 2018; Jordt et al., 2020). Indeed, mathematical models that incorporate more complex factors, including fluctuating selective pressures (Peña-Miller et al., 2015), plasmid compensation (Harrison et al., 2016; Hall et al., 2017; Zwanzig et al., 2019), population diversity (i.e., multiple strains) (Dionisio, 2005), and pleiotropy (Jordt et al., 2020) expand the range of conditions that allow plasmid maintenance relatively to earlier models that focused primarily on essential plasmid features (fitness costs, conjugative efficiency, and loss rate) (Stewart and Levin, 1977; Simonsen, 1991; Bergstrom et al., 2000).

Besides experiencing changing environments or the presence of additional strains, bacterial communities may also face several plasmids simultaneously. Indeed, bacteria commonly carry more than one plasmid (San Millan et al., 2014a), which allows interactions among them (reviewed in Dionisio et al., 2019). For example, interactions between plasmids that alter fitness costs (epistasis) (Silva et al., 2011; San Millan et al., 2014a) affect the total cost of harboring the two plasmids simultaneously such that it becomes higher (or lower) than simply the sum of the cost imposed by each plasmid individually. Plasmid behavior may also change whether they occupy the same cells or if they are present in different cells of the same population (Sagai et al., 1977; Gama et al., 2017a). Recently it has also been shown that upon host-plasmid co-evolution under antibiotic selection, the stability of two distinct plasmids in communities consisting of two different bacterial species increased once antibiotics were removed (Jordt et al., 2020). Interestingly, pleiotropic effects resulted in greater plasmid persistence in both novel host-plasmid combinations and, in some cases, multi-plasmid hosts (Jordt et al., 2020).

In this work, we modify earlier mathematical models devoted to studying plasmid maintenance (Stewart and Levin, 1977; Levin et al., 1979; Simonsen, 1991) as to incorporate multiple plasmids and interactions among them. Such interactions affect the three main features of plasmid biology: fitness cost, conjugative efficiency, and rate of loss. Therefore, we aim to understand how interactions between plasmids affect their maintenance in bacterial populations.

## Materials and Methods

### Parameters

To model the ecology of plasmid-harboring bacteria we consider three main parameters of plasmid biology: the fitness effect of plasmid carriage on the host cell, the rate of conjugative transfer and the rate of plasmid loss due to missegregation. These models are adapted from mass-action models by Stewart and Levin (1977), Levin et al. (1979), Simonsen (1991), but we now consider two plasmids X (focal) and Y (competitor). We assume that all cells are isogenic, only differing in plasmid content. We also assume that there is no evolution of hosts or of plasmids, otherwise pleiotropic effects between mutations in different replicons should be considered (Jordt et al., 2020). For language simplification, X may refer to plasmid X, to a cell harboring plasmid X, or even the density of cells harboring plasmid X, depending on the context. The same definitions apply to Y. In a bacterial population where cells carry each of two different conjugative plasmids, there are four possible types of cells: X and Y carrying, respectively, each of the plasmids, XY carrying simultaneously both plasmids as a result of conjugation, and plasmid-free cells Ø produced by plasmid loss. Following the simplification mentioned above, XY may refer to cells harboring both the X and Y plasmids or the density of cells harboring both plasmids, while Ø may refer to plasmid-free cells or their density. The parameters and their values are presented in Table 1.

#### Fitness

*Escherichia coli* K-12 has a doubling time of 20 min in Lysogeny Broth (Sezonov et al., 2007), equivalent to three generations per hour, which translates to a standard growth rate ψ = 3⋅ln(2) h^{–1}. We assume this to be the maximum growth rate of plasmid-free strains. Since plasmid carriage commonly entails a fitness cost, plasmid-carrying strains will grow at a rate ψ weighed by their fitness ω (relative to the plasmid-free strain), such that the growth rate of a strain carrying plasmid X is ψ_{X} = ψ⋅ω_{X}. The fitness cost of different conjugative plasmids has been determined in the ranges 3.9–8% and 0–14.3%, respectively, in *E. coli* (Silva et al., 2011) and *P. aeruginosa* (San Millan et al., 2014a). Accordingly, we employed fitness costs of 2.5%, 5%, 10%, and 15% on the models, translating into ω_{X} (and ω_{Y}) ∈ {0.85; 0.9; 0.95; 0.975}. Although other works also report fitness costs for conjugative plasmids, we focused solely on the two above as they also provide another important information for this study, the fitness effects of the interactions (epistasis) between plasmid pairs.

Among conjugative plasmid combinations, positive and negative epistatic interactions (ε) are equally frequent and within a 95% confidence interval of [−0.035; 0.06] (Silva et al., 2011). Epistasis between a conjugative (or a mobilizable) and a non-conjugative plasmid tended to be positive, varying between 0.005 and 0.159 (San Millan et al., 2014a). Therefore, we considered values of epistasis ε∈ {−0.05; 0; 0.05; 0.1}, such that the fitness of strains carrying both plasmids X and Y is ω_{XY} = ω_{X}⋅ω_{Y} + ε.

#### Loss

Plasmid-carrying cells can become plasmid-free due to plasmid loss during cell division. This rate of plasmid loss can vary substantially, at least ranging from 10^{–9} to 10^{–3} (Lau et al., 2013; Loftie-Eaton et al., 2017) and we considered it as δ∈ {10^{–8}; 10^{–6}; 10^{–4}}. We considered that plasmids X and Y always have the same loss rate, to reduce the number of possible combinations that would otherwise be too high. Therefore, a proportion δ⋅X of cells carrying plasmid X becomes plasmid-free every hour.

The rate of plasmid loss can be affected by incompatible plasmids that interfere with replication or partition mechanisms (Novick, 1987; Nordström and Austin, 1989), but also disturbed by compatible plasmids that may perturb stability mechanisms such as toxin–antitoxin (also known as post-segregation killing or addiction) systems (Kamruzzaman et al., 2017). Interference between addiction systems can increase plasmid loss at least one order of magnitude (Radnedge et al., 1997), thus we regarded interference on loss rate as σ∈ {1; 10}, such that when the two plasmids are present in the same cell the loss rate for either of them is δ_{XY} = δ⋅σ.

#### Conjugation

The remaining parameter is plasmid transfer which occurs at a specific conjugation rate γ and depends on the proportion of donor (*D*) and recipient (*R*) cells such that transconjugants (*T*) emerge as $\frac{dT}{dt}=\mathrm{\gamma}\cdot D\cdot R$. Values for plasmid conjugation rates were based on our previous publication (Gama et al., 2017a). These conjugation efficiencies were converted by the end-point method (Simonsen et al., 1990) to reflect conjugation rates per hour (they are an approximation as our experiments lasted 90 min). Conjugation rates varied between 10^{–17} to 10^{–8}. Transfer of plasmids R477-1 (temperature-sensitive transfer) and R6K was near the detection limit as transconjugant colonies were only found occasionally; exclusion of these plasmids revealed a lower bound of 5 × 10^{–13}. By contrast, plasmids R1drd19 and F (which are de-repressed for conjugation) and R124 displayed conjugation efficiencies above 10^{–10} characteristic of de-repressed plasmids (Gordon, 1992; Dionisio et al., 2002) and were thus ignored such that their exclusion revealed a higher bound of 5 × 10^{–11}. Therefore, we considered conjugation rates γ_{X} (and γ_{Y}) ∈ {10^{–13}; 10^{–12}; 10^{–11}} which are within limits previously reported for natural plasmids (Gordon, 1992).

Conjugation rates can be affected by plasmid interactions either if the plasmids reside in the same cell (intracellular effect α) or if one is present in the donor and the other in the recipient cell (intercellular effect ξ). Each of these effects can be negative (inhibiting plasmid transfer) or positive (facilitating plasmid transfer). We calculated these effects by dividing the endpoint conjugation rate of a plasmid when in combination by the respective rate when alone. The median positive intracellular effect was 21.9 and the median negative intracellular effect was 0.006. Therefore, we considered α∈ {10^{–3}; 1; 10} as the conjugation rate of X in the presence of Y in the same cell [γ_{X(Y}_{)}] and the conjugation rate of Y in the presence of X [γ_{Y}_{(}_{X}_{)}]. The conjugation rate of plasmid X from XY cells to Ø cells (plasmid-free cells) becomes γ_{X(Y)} = γ_{X}⋅α_{X} while the conjugation rate of plasmid Y from XY cells to Ø cells becomes γ_{Y(X)} = γ_{Y}⋅α_{Y}, such that α_{X} and α_{Y} are the respective intracellular effects on X and Y. The median positive intercellular effect was 8.3 while the median negative intercellular effect was 0.04, such that we considered ξ∈ {10^{–2}; 1; 10}. The conjugation rate from X cells to Y cells thus becomes γ_{X}_{Y} = γ_{X}⋅ξ_{X} and from Y into X becomes γ_{Y}_{X} = γ_{Y}⋅ξ_{Y}, where ξ_{X} and ξ_{Y} are the respective intercellular effects on X and Y.

Three types of transconjugants can develop from matings between cells XY with Ø that may acquire either plasmid X or Y but also acquire both plasmids. Simultaneous transfer of plasmids has been shown to occur at an identical rate as the transfer of the least efficient plasmid in the combination (Gama et al., 2017b), thus we estimated co-transfer γ_{XY} as the minimum between γ_{X(Y)} and γ_{Y(X)}. As a consequence of these matings, Ø cells receive both plasmids at a rate γ_{XY} but plasmid X and Y individually at rates γ_{X(Y)} − γ_{XY} and γ_{Y(X)} − γ_{XY}, respectively.

Additionally, strains carrying a single plasmid can receive the second one from XY cells. Transconjugants from matings between XY and X cells and between XY and Y cells should occur, respectively, as γ_{X(Y)}_{Y} = γ_{X(Y)} and γ_{Y(X)}_{X} = γ_{Y(X)} since they depend on the intracellular interactions between plasmids co-residing in the donor XY cells. However, we consider this to be true only when the intracellular effect decreases the rate of transfer (α ≤ 1). The rationale for this exception is the following. Facilitation occurs because the conjugative pili expressed by the second plasmid help to stabilize the mating pair (Gama et al., 2017a). However, in matings between XY and X cells or between XY and Y cells, one of the plasmids is present in both donor and recipient. Therefore, surface/entry exclusion mechanisms expressed by the plasmid in the recipient cell prevent the same plasmid in the donor cells to transfer efficiently (Garcillan-Barcia and de la Cruz, 2008). Exclusion tends to have a stronger effect on mating efficiency than facilitation (Garcillan-Barcia and de la Cruz, 2008), thus neutralizing it. Therefore, we consider that γ_{X(Y)}_{Y} = γ_{X} and γ_{Y(X)}_{X} = γ_{Y} when α > 1.

### General Model

The models employed in this study were adapted from previous works (Stewart and Levin, 1977; Levin et al., 1979; Simonsen, 1991) that describe chemostat cultures. The differential equations 1–5 define the rate of change (derivative in time) of the concentration of nutrients and the density of the cell types.

The rate of change of the density of plasmid-free cells at a specific time is given by Eq. (1). The first term defines population growth as the product of population density of plasmid-free cells (Ø) by its growth rate (ψ) and the Monod function $\frac{R}{R+Q}$, which models bacterial growth in a liquid environment relative to the concentration of a limiting nutrient. The second term defines the number of cells that exit the chemostat per volume unit of the chemostat at a washout rate Ω. The third term defines the density of the population of single-plasmid-carrying cells (X or Y) that lose the plasmid at a rate δ. The following terms define the proportion of the density of the population of plasmid-free cells that become transconjugants on different matings, which depends on the rate of conjugation and the population densities of both donor and recipient cells (in this case Ø). The two first conjugation terms describe matings between the plasmid-free recipient cells with donor cells carrying either plasmid X or Y, respectively. The last conjugation term describes matings with donor cells carrying simultaneously both plasmids X and Y, such that three types of transconjugates occur at different rates: γ_{XY} for transconjugants receiving both plasmids, γ_{X(Y)} − γ_{XY} and γ_{Y(X)} − γ_{XY} for those receiving, respectively, only plasmid X and Y.

Equation 2 relates to the population density of cells carrying plasmid X only. The first term corresponding to cell growth is similar to the one in Eq. (1), but takes the product by the strain fitness (ω_{X}) to weigh the effect of carrying the plasmid. Second and third terms reflect, respectively, cell washout and loss of plasmid X becoming plasmid-free. Additionally, cells carrying both plasmids can lose plasmid Y at a rate δ_{XY}, thus creating cells carrying only plasmid X. The amount of X cells also increases due to cells Ø that receive only plasmid X from matings with either cells X or XY, but decreases due to cells X that receive plasmid Y from matings with either cells Y or XY at rates γ_{Y}_{X} and γ_{Y(X)}_{X}, respectively.

Equation 3 is identical to Eq. (2), but refers to cells carrying only plasmid Y instead.

Equation 4 describes the kinetics of cells carrying simultaneously both plasmids X and Y and resembles Eqs (2 and 3), except the loss term is doubled (as cells can lose either of the plasmids) and that all conjugation terms lead to an increase in the population density of cells with both plasmids. Note that for simplicity we consider that cells can lose one plasmid at a time, but never simultaneously.

Lastly, Eq. (5) describes the kinetics of nutrients in the chemostat, such that at a specific time point their concentration depends on the balance between the amount of nutrients (per volume unit) that enters (Ω⋅*R*_{0}) and exits (Ω⋅*R*) the chemostat minus the nutrients exhausted for cell growth. The latter being proportional to the amount of resources required per cell division concentration (Π) and the cumulative growth of each of the different strains.

### Model Implementation

We performed simulations under different models. In model 1, plasmid X was the only plasmid in the population, while in models 2–6 plasmid Y is also present. In model 2, the two plasmids do not interact, while in models 3–6 the plasmids interact according to a single interaction parameter (ε, α, ξ, or σ).

The models were implemented in R (R Core Team, 2018) version 3.4.4 and ran with package rootSolve (Soetaert, 2009) to analyze the steady-state of ordinary differential equations and with package deSolve (Soetaert et al., 2010) for solving the equations across time to find the time required for extinction. Packages doParallel (Microsoft Corporation and Weston, 2019) and foreach (Microsoft and Weston, 2019) were used to run multiple simulations in parallel.

### Analysis

The results obtained for the steady-state are on a continuous scale, such that the decreasing population density can become < 1 without reaching 0 when the population met the stable equilibrium. Thus, we considered the plasmid to go to extinction if there is less than one cell per volume unit (mL, see Table 1) carrying plasmid X (i.e., X < 1 and XY < 1) at the steady-state, otherwise the plasmid is stably maintained which we define here as survival. When the plasmid did not survive, we estimated the time of extinction (ToE) as the first time point (non-negative integer scale) where X < 1 and XY < 1. For model 1, where the populations only carried plasmid X but not Y, these conditions simply become X < 1.

We analyzed the effect of each plasmid interaction (models 3–6) relatively to the “null” models 1 and 2. The analysis against model 1 allows us to understand if the presence of cells carrying a second plasmid Y affects the persistence of the focal plasmid X. For this, we normalized the time of extinction obtained from the model under study by the respective value obtained with model 1. That is, we divided the time of extinction under that model by the respective time of extinction under model 1. After normalization, values > 1 indicate that the plasmid increases its time of extinction, and vice-versa. We quantitatively compared (Kruskal–Wallis or Wilcoxon tests) these values in function of the different parameter values of the interaction to understand how the different parameter values affected the outcome of plasmid X. Next, we did a complementary qualitative analysis. To do so, we categorized the normalized times of extinction as the following outcomes: “increase” for values > 1, “decrease” when < 1 and “null” when = 1. Then we applied χ^{2} tests on the tabulated outcomes per interaction value. We performed equivalent analyses against model 2, which allows us to evaluate how favorable (or not) it is for the focal plasmid X to interact with the co-resident plasmid Y.

Data analyses were performed in R (R Core Team, 2018) version 3.5.1 using the following packages: reshape2 (Wickham, 2007) and data.table (Dowle and Srinivasan, 2019) to handle data structures, heplots (Friendly, 2010) to calculate partial η^{2}, FSA (Ogle et al., 2019) and rcompanion (Mangiafico, 2019) for Kruskal–Wallis *post hoc* testing, pscl (Jackman, 2020) and performance (Lüdecke et al., 2020) to calculate logistic regression’s pseudo-*R*^{2}. Figures were created using packages ggplot2 (Wickham, 2016), patchwork (Pedersen, 2020), ggthemes (Arnold, 2019), and ggpubr (Kassambara, 2020). Supplementary tables were created using package openxlsx (Schauberger and Walker, 2019).

## Results

### Model 1: One Plasmid Only

We started the study by evaluating the fate of different plasmids when alone in the bacterial population. In this model, the initial population consists of 5 × 10^{5} X cells (as explained in the “Materials and Methods” section, and an X cell is a cell harboring plasmid X). Plasmid loss creates plasmid-free cells, which can later reacquire it through conjugation.

We studied the fate of 36 plasmids, each consisting of a different combination of features (fitness, conjugation rate, and loss rate). None of these plasmids could be maintained in the chemostat. Therefore, we analyzed their time of extinction (ToE) as a function of the three plasmid features and interactions among them. Time of extinction depends on fitness, conjugation and loss rates but also on the combined effect of fitness and conjugation rate (ANOVA, *df* = 7, *P* < 0.05, Supplementary Table S1). As shown in Figure 1, plasmid-carrying cells take longer to go extinct with increasing fitness and conjugation rates but with decreasing loss rates, as it would be expected. The combined effects of fitness and conjugation rate is illustrated as cells with γ_{X} = 10^{–11} and ω_{X} = 0.975 are not parallel to the others in Figures 1A and B, respectively, while cells with ω_{X} = 0.975 and γ_{X} = 10^{–11} display much higher time of extinction in Figure 1C. This is why, in the ANOVA, the effect size (η^{2}) of the fitness-conjugation rate interaction is high – in fact even higher than that of loss (Supplementary Table S1).

**Figure 1.** Effect of fitness (ω), conjugation (γ), and loss (δ) rates on the plasmid time of extinction (ToE). Results from model 1, comprising a single plasmid. Time of extinction (ToE) is the first time point where there is less than one plasmid-carrying cell in the population. Time of extinction represented in order to fitness **(A)**, conjugation **(B)**, and loss rate **(C)**.

### Model 2: Two Non-interacting Plasmids

From now on, we model populations starting with both X and Y cells such that plasmid-free cells will be created by plasmid loss and XY cells by conjugation. In this model 2, as well as the following models 3 to 6, we simulated populations starting with 5 × 10^{5} X cells and 5 × 10^{5} Y cells, such that the initial population size is 10^{6} cells. To control for this change, we repeated the simulation of a single plasmid (model 1) considering the initial population as 10^{6} X cells instead of 5 × 10^{5} X cells. The result is qualitatively the same since no plasmid survives (i.e., no plasmid was stably maintained), while, quantitatively, the time of extinction only increased by 1 h or remained unchanged.

In model 2, although all mating combinations are possible, the interactions between plasmids have no effect (ε = 0 and α = ξ = σ = 1). Both focal (X) and competitor (Y) plasmids express all 36 feature combinations, which consequently results in 432 combinations of two plasmids. In all cases the focal plasmid X always becomes extinct, as well as plasmid Y. In 36 cases the two plasmids become extinct at the same time, which is when they display the same characteristics, while in 198 cases plasmid X disappears before Y. Among those 198 cases, plasmid X becomes extinct faster than plasmid Y when imposing higher fitness cost or when having the same fitness effect but lower conjugation rate (note that plasmids X and Y always have identical loss rates). In the remaining 198 cases, Y cells disappear before X cells, since plasmid features are inverted.

We compared the times of extinction for plasmid X obtained from models 1 and 2 to understand if the presence of cells carrying a second plasmid Y affects its persistence (Figure 2). When a second plasmid is present the time of extinction decreases (paired Wilcoxon test with continuity correction, *P* < 2.2 × 10^{–16}). This shows that competition with another plasmid-carrying strain accelerates plasmid extinction.

**Figure 2.** A competitor plasmid Y accelerates the extinction of the focal plasmid X. The vertical axis represents the time of extinction (ToE) of plasmid X when alone in the bacterial population (model 1) and the horizontal axis represents ToE of plasmid X when a competitor plasmid Y is present in the same population (model 2). The blackline represents *y* = *x*.

### Model 3: Interactions Affecting Fitness (Epistasis)

In model 3, the two plasmids can interact, affecting the total fitness effect on XY cells. Thus, ω_{XY} = ω_{X}⋅ω_{Y} + ε, where ε∈ {−0.05; 0; 0.05; 0.1}. This data set consists of a total of 1728 cases, 432 per value of ε.

#### Plasmid Survival

The focal plasmid can survive in 128 of these cases, 105 where epistasis is positive with ε = 0.1 and 23 cases if ε = 0.05. These positive values of epistasis imply that, in all these cases, the fitness of XY cells is higher than that of plasmid-free cells and they are always present in the final population at higher frequencies than X or Y cells. Interestingly, in 45 cases there are no plasmid-free cells in the final population. This happens when loss rates are very low, namely δ = 10^{–8}. In the remaining 83 cases, the maximum frequency of plasmid-free cells is 18.66%.

Yet, there are seven more cases where the fitness of cells carrying both plasmids is higher than that of plasmid-free cells (ω_{XY} ≥ 1) but the focal plasmid goes to extinction. All these exceptions have loss rate δ = 10^{–4}, four consisting in positive epistasis of ε = 0.05 and ω_{X} = ω_{Y} = 0.975 and the remaining three in ε = 0.1 and ω_{X} = ω_{Y} = 0.95. Additionally, in these seven cases, γ_{X} and γ_{Y} < 10^{–11} (but not γ_{X} = γ_{Y} = 10^{–12} when ε = 0.1). These results show that the plasmids can survive if there is positive epistasis, elevating the fitness XY cells above that of the plasmid-free cells. Nonetheless, there is no guarantee of plasmid survival as lower conjugation rates preclude the emergence of XY cells and higher loss rates decrease their maintenance.

As just mentioned, in 135 cases carriage of two plasmids confers a fitness benefit to the host cells relatively to plasmid-free cells. Despite that, in model 1, we only considered plasmids that produce fitness costs. As a control, we expanded model 1 with ω_{X} ≥ 1, such that the focal plasmid confers no fitness cost or even a benefit (thus ω_{X} ∈ {1; 1.05}). We considered a fitness benefit of 5% because, when XY cells had a fitness higher than that of plasmid-free cells, the fitness effect ranged between 0.000625 and 0.050625. Among the new 18 simulated cases, there are two where plasmid X goes extinct – when ω_{X} = 1 and γ_{X} < 10^{–11} and δ = 10^{–4}. These two cases of extinction further illustrate the conceivable detrimental effect of conjugation and loss rates. Moreover, even in the cases where the plasmid can be maintained, the population at equilibrium is not pure and plasmid-free cells exist at a frequency varying from 2.96 × 10^{–6} to 1 × 10^{–1}. Overall, this result shows that fitness alone is not sufficient for plasmid survival and depends on the rates of conjugation and loss. It is also interesting to note that plasmid-free cells still persist in these conditions, unlike what was observed for several cases when another plasmid is present in the population.

#### Time of Extinction

Among the 1600 cases where the focal plasmid does not survive there are 736 cases of positive epistasis. In 135 of these cases, XY cells are still less fit than both X and Y cells and are the first to get extinct. In 223 cases XY cells are fitter than both X and Y, but they are the ones displaying the longest time of extinction in only 58 cases. This again shows that XY cells are not necessarily stable even if they are fitter than X and Y. Therefore, even when there is positive epistasis between plasmids, conjugation and loss rates are still determinant for plasmid maintenance.

Among the 432 cases of negative epistasis, XY cells are always less fit than either X or Y cells. Nonetheless, in 12 cases, XY cells are not the first getting extinct – if γ_{X} = γ_{Y} = 10^{–11} and δ > 10^{–4} and either ω_{X} or ω_{Y} (but not both) is 0.975. Thus, provided that both plasmids exhibit high conjugation and low loss rates and one of them imposes low fitness costs XY cells can persist longer that one of the single-plasmid-carrying strains.

We compared these results with those obtained with model 2 to unveil how epistatic interactions affect the fate of the focal plasmid relative to when the two plasmids did not interact (Figure 3A). Therefore, we analyzed the 1600 times of extinction after normalizing them by the results of model 2 – which is equivalent to divide by the respective value when there was no epistasis (ε = 0).

**Figure 3.** Relative time of extinction when plasmid interactions affect fitness (epistasis; ε). **(A)** Relative time of extinction (ToE) in presence and absence of epistasis: the vertical axis represents ToE of plasmid X when there is epistasis (model 3) and the horizontal axis represents the respective ToE when there are no plasmid interactions (ε = 0; model 2). **(B)** Relative ToE in presence and absence of a competitor plasmid Y: the vertical axis represents ToE of plasmid X when Y is present and plasmid interactions affect fitness (model 3) while the horizontal axis represents the respective ToE when plasmid X is alone in the bacterial population (model 1). The analysis is always in function of the focal plasmid X. The blackline represents *y* = *x*.

We analyzed the normalized data in function of the different values of epistasis. This type of interaction has a significant effect (Kruskal–Wallis test, *df* = 3, *P* < 2.2 × 10^{–16}, Supplementary Table S2), and all four epistasis values produce distinct results (Dunn test with Bonferroni correction, Supplementary Table S2). We further analyzed the effects qualitatively, i.e., we checked if times of extinction increased, decreased or did not change. These results show that the four epistasis values produce different results (χ^{2} test with Yates continuity correction, *df* = 6, *P* < 2.2 × 10^{–16} followed by Fisher pairwise tests of independence, Supplementary Table S2), which agrees with the conclusion of the quantitative analysis. Altogether the results correspond to the prediction that the times of extinction increase with the strength of positive epistasis and decrease with negative epistasis. However, it is important to mention that it is not straightforward to conclude the outcome based on epistasis only because the result also depends on the features of the two plasmids, namely fitness impact, conjugative transfer and loss rate. This means that with positive epistasis not all cases reveal an increased time of extinction, and the reverse is also true for negative epistasis. In addition, ε = −0.05 (when normalized relatively to ε = 0) has a much smaller effect (median = 1 and standard deviation = 0.03) than ε = 0.05 (median 1 and standard deviation = 1.25), showing that negative epistasis has a less detrimental effect than positive epistasis has a beneficial effect.

We also analyzed the 1600 cases to understand how the times of extinction differ from when the competitor plasmid is absent (Figure 3B), that is, we compared them with those obtained with model 1. We analyzed this effect by normalizing the time of extinction of the focal plasmids by the respective value obtained when they were alone in the population (model 1). Positive epistasis has a significant effect (Kruskal–Wallis test, *df* = 3, *P* < 2.2 × 10^{–16}, followed by *post hoc* Dunn test with Bonferroni correction, Supplementary Table S3), such that there is a stronger effect for ε = 0.1 and then for ε = 0.05; but ε = −0.05 and ε = 0 do not differ from each other. We obtained the same conclusion when analyzing the effects qualitatively (χ^{2} test with Yates continuity correction, *df* = 3, *P* < 2.2 × 10^{–16} followed by Fisher pairwise tests of independence, Supplementary Table S3), i.e., comparing the number of cases among the categories where times of extinction increase, decrease or do not change. Thus, negative epistasis does not lead to a significantly different outcome than not having epistasis, and the plasmid always gets extinct earlier than when alone. Nevertheless, the number of cases of increased plasmid persistence increases with the strength of positive epistasis (Supplementary Table S3).

### Model 4: Intracellular Interactions Affecting Conjugation

In this model, the two plasmids can interact, such that the conjugation rate of one of the plasmids changes when both plasmids occupy the same cell – an intracellular conjugation effect α. There are two possibilities: either X or Y is the target of the effect, but never both. Thus, either the transfer rate of X when co-inhabiting with Y is the transfer rate of X when alone times α_{X} (γ_{X(Y)} = γ_{X}⋅α_{X}) or the transfer rate of Y when co-inhabiting with X is the transfer rate of Y when alone times α_{Y} γ_{Y(X)} = γ_{Y}⋅α_{Y}), where α_{X} or α_{Y} ∈ {10^{–3}; 1; 10}. Next, we analyze both possibilities, but always in function of plasmid X. Each of these data sets consists of a total of 1296 cases, 432 per value of the intracellular conjugation effect α.

#### Plasmid X as the Target of Interactions

The focal plasmid cannot survive in any of the 1296 cases when it was the target of intracellular conjugation effects. Therefore, we analyzed the times of extinction, comparing among the three values of the intracellular conjugation effect α_{X} (Figure 4A). These values were normalized by the respective values when α_{X} = 1 (which is equivalent to divide by the respective values of model 2). There is a significant effect (Kruskal–Wallis test, *df* = 2, *P* < 2.2 × 10^{–16}, Supplementary Table S4) revealing differences between the three α_{X} parameter values (Dunn test with Bonferroni correction, Supplementary Table S4). Even though conjugation rate increases only 10 times, it has a stronger impact (median = 1, standard deviation = 0.14) than when it decreases 1000 times (median = 1, standard deviation = 0.05). We further analyzed the effects qualitatively, i.e., whether times of extinction increased, decreased or did not change. The three α_{X} parameter values produce different results (χ^{2} test with Yates continuity correction, *df* = 4, *P* < 2.2 × 10^{–16} followed by Fisher pairwise tests of independence, Supplementary Table S4), which agrees with the quantitative analysis. While with α_{X} = 10^{–3} times of extinction either decrease (198 cases) or remain the same, with α_{X} = 10 they increased in 247 cases but never decreased.

**Figure 4.** Relative time of extinction when plasmid intracellular interactions (α) affect conjugation rates. **(A,C)** Relative time of extinction (ToE) in presence and absence of interactions: the vertical axis represents ToE of plasmid X when the plasmids interact (model 4) and the horizontal axis represents the respective ToE when the plasmids do not interact (α = 1; model 2). **(B,D)** Relative ToE in presence and absence of a competitor plasmid Y: the vertical axis represents ToE of plasmid X when Y is present and plasmid intracellular interactions affect conjugation rates (model 4) while the horizontal axis represents the respective ToE when plasmid X is alone in the bacterial population (model 1). In panels **(A,B)**, plasmid X is the target of the effect α_{X}, while in panels **(C,D)** plasmid Y is the target of effect α_{Y} (models 4.1 and 4.2, respectively). The analysis is always in function of the focal plasmid X. The blackline represents *y* = *x*.

Next, we compared how times of extinction differ from those obtained with model 1 when the competitor plasmid was absent (Figure 4B). Only α_{X} = 10 has a significant effect (Kruskal–Wallis test, *df* = 2, *P* < 1.3 × 10^{–6}, followed by *post hoc* Dunn test with Bonferroni correction, Supplementary Table S5), while α_{X} = 10^{–3} and α_{X} = 1 are not significantly different. Then, we analyzed the effects qualitatively. While all cases with α_{X} = 10^{–3} and α_{X} = 1 the time of extinction decreased relatively to model 1, α_{X} = 10 led to increase in 46 cases (χ^{2} test with Yates continuity correction, *df* = 2, *P* < 2.2 × 10^{–16} followed by Fisher pairwise tests of independence, Supplementary Table S5). In all these 46 cases, the focal plasmid had γ_{X} = 10^{–11} and α_{X} = 10 such that in cells carrying both plasmids γ_{X(Y)} = 10^{–10}. Therefore, we tested the outcome of plasmids with γ_{X} = 10^{–10} when alone in a population (model 1). This revealed that all 12 plasmids could be stably maintained in the population, such that the frequencies of plasmid-free cells would vary between 1.15 × 10^{–6} and 7.8 × 10^{–2}.

In conclusion, only an increasing intracellular effect on conjugation affected times of extinction significantly, although this was not sufficient to allow plasmid survival in the presence of a second plasmid.

#### Plasmid Y as the Target of Interactions

When the competitor plasmid was the target of intracellular conjugation effects the focal plasmid could not survive in any of the 1296 cases. We analyzed the time of extinction of the focal plasmid X in terms of the intracellular conjugation effect α_{Y} (Figure 4C) relative to model 2. There is a significant effect and the three α_{Y} parameter values produce significantly different results (Kruskal–Wallis test, *df* = 2, *P* < 2.2 × 10^{–16}, followed by Dunn test with Bonferroni correction, Supplementary Table S6). There is also a significant difference when analyzing the effects qualitatively (χ^{2} test with Yates continuity correction, *df* = 4, *P* < 2.2 × 10^{–16} followed by Fisher pairwise tests of independence, Supplementary Table S6). The time of extinction of plasmid X increased in 211 cases when α_{Y} = 10^{–3} but in 187 when α_{Y} = 10. However, there are also cases in which the time of extinction decreased: 3 cases when α_{Y} = 10^{–3} and 65 when α_{Y} = 10. Interestingly, the quantitative and qualitative analyses reveal different aspects of these results. Both results with α_{Y} = 10^{–3} and α_{Y} = 10 have a median of 1. Though α_{Y} = 10 allows the time of extinction to increase in fewer cases, it has a stronger effect in the sense that it produces a higher standard deviation (0.06 vs. 0.02). Nonetheless, unlike the previous section, the survival of the focal plasmid can increase whether the effect is below or over 1.

There was no significant difference related to the effect α_{Y} (Kruskal–Wallis test, *df* = 2, *P* = 0.06, Supplementary Table S7) when analyzing times of extinction relative to model 1 where the competitor plasmid was absent (Figure 4D). However, times of extinction always decreased with the exception of 21 cases with α_{Y} = 10 which were sufficient to produce a significant difference (χ^{2} test with Yates continuity correction, *df* = 2, *P* = 5.37 × 10^{–10} followed by Fisher pairwise tests of independence, Supplementary Table S7).

### Model 5: Intercellular Interactions Affecting Conjugation

As in model 4, the two plasmids can interact but the conjugation rate of only one of them changes when both donor and recipient cells carry plasmids – an intercellular conjugation effect ξ. Here there are again two non-overlapping possibilities: either X or Y is the target of the effect, that we analyze only in function of the focal plasmid X. Thus, either γ_{X}→_{Y} = γ_{X}⋅ξ_{X} or γ_{Y}→_{X} = γ_{Y}⋅ξ_{Y}, where ξ∈ {10^{–2}; 1; 10}. Each of these data sets consists on a total of 1296 cases, 432 per value of *ξ*.

#### Plasmid X as the Target of Interactions

The focal plasmid does not persist in any of the 1296 cases where it was the target of intercellular conjugation effects. We normalized times of extinction relatively to ξ_{X} = 1 (model 2) to evaluate the effect of the different values of ξ_{X} (Figure 5A). The three values of ξ_{X} produce significantly different results (Kruskal–Wallis test, *df* = 2, *P* < 2.2 × 10^{–16}, followed by Dunn test with Bonferroni correction, Supplementary Table S8). An increase of conjugation rate with ξ_{X} = 10 displays a slightly stronger impact (median = 1, standard deviation = 0.06) than a decrease with ξ_{X} = 0.01 (median = 1, standard deviation = 0.04). Moreover, with ξ_{X} = 0.01 the time of extinction either decreases (313 cases) or remains the same while with ξ_{X} = 10 the time of extinction increased in most cases (408/432) and did not change in the remaining 24 (in all of which γ_{X} = 10^{–13}). This difference in outcomes is significant (χ^{2} test with Yates continuity correction, *df* = 4, *P* < 2.2 × 10^{–16} followed by Fisher pairwise tests of independence, Supplementary Table S8).

**Figure 5.** Relative time of extinction (ToE) when plasmid intercellular interactions (ξ) affect conjugation rates. **(A,C)** Relative time of extinction (ToE) in presence and absence of interactions: the vertical axis represents ToE of plasmid X when the plasmids interact (model 5) and the horizontal axis represents the respective ToE when the plasmids do not interact (ξ = 1; model 2). **(B,D)** Relative ToE in presence and absence of a competitor plasmid Y: the vertical axis represents ToE of plasmid X when Y is present and plasmid intercellular interactions affect conjugation rates (model 5) while the horizontal axis represents the respective ToE when plasmid X is alone in the bacterial population (model 1). In **(A,B)** plasmid X is the target of the effect ξX while in **(C,D)** plasmid Y is the target of effect ξY (models 5.1 and 5.2 respectively). The analysis is always in function of the focal plasmid X. The blackline represents *y* = *x*.

In Figure 5B, we show how times of extinction differ from when the focal plasmid is alone in the population (model 1). There is a significant difference (Kruskal–Wallis test, *df* = 2, *P* = 8.2 × 10^{–7}, followed by *post hoc* Dunn test with Bonferroni correction, Supplementary Table S9), such that ξ_{X} = 10 displays a higher median. On the other hand, a qualitative analysis revealed no significant differences (χ^{2} test with Yates continuity correction, *df* = 2, *P* = 0.05) such that the time of extinction always decreases, except in three cases (where ω_{X} = 0.9, γ_{X} = 10^{–11}, δ < 10^{–8} and ξ_{X} = 10).

In conclusion, an increasing intercellular conjugation effect can increase the time of extinction but with a very limited impact and still not allowing plasmid maintenance.

#### Plasmid Y as the Target of Interactions

The focal plasmid could not survive in any of the 1296 cases when the competitor plasmid was the target of intercellular conjugation effects. Figure 5C displays times of extinction in terms of the intercellular conjugation effect ξ_{Y}. There is a significant effect and the three ξ_{Y} parameter values produce significantly different results (Kruskal–Wallis test, *df* = 2, *P* < 2.2 × 10^{–16}, followed by Dunn test with Bonferroni correction, Supplementary Table S10). However, times of extinction of plasmid X tend to increase when ξ_{Y} < 1 and to decrease when ξ_{Y} > 1, which is the opposite outcome of when plasmid X is the target of interaction instead of Y. Therefore, preventing the acquisition of plasmid Y reduces the rate of formation of less fit XY cells, that would otherwise accelerate extinction of plasmid X as they are more easily outcompeted in the population. When analyzed qualitatively ξ_{Y} > 1 leads to a decrease in all cases while ξ_{Y} < 1 increases times of extinction in 393 out of 432 cases (χ^{2} test with Yates continuity correction, *df* = 4, *P* < 2.2 × 10^{–16}, followed by Fisher pairwise tests of independence, Supplementary Table S10).

Next, we compared times of extinction with those from model 1, which lacks the competitor plasmid (Figure 5D). Unlike when plasmid X was itself the target of the intercellular interaction, there are significant effects if plasmid Y is the target instead (Kruskal–Wallis test, *df* = 2, *P* < 2.2 × 10^{–16}, Supplementary Table S11). This shows that if X cells acquired plasmid Y at lower rates then plasmid X can be maintained for longer times in the population, while if acquired at higher rates plasmid Y will have a detrimental effect on plasmid X’s persistence. The qualitative analysis reveals the same conclusion (χ^{2} test with Yates continuity correction, *df* = 2, *P* = 3 × 10^{–12} followed by Fisher pairwise tests of independence, Supplementary Table S11). With ξ_{Y} = 0.01 there are 26 cases where ToE increased, all of which have γ_{Y} = 10^{–11}. For the other values of ξ_{Y} tested, times of extinction decreased in all cases.

Therefore, plasmid X can persist longer if it reduces the rate of acquiring plasmid Y.

### Model 6: Interactions Affecting Loss

In this model, the plasmid interactions lead to increased loss rates, such that δ_{XY} = δ⋅σ, where σ∈ {1; 10}. This data set consists on 864 total cases, 432 per value of σ. The plasmid cannot survive in any of these cases, so we only analyzed the times of extinction. When normalized relatively to model 2, the times of extinction for plasmid X can vary significantly (Wilcoxon test, *P* < 2.2 × 10^{–16}, Supplementary Table S12) but with little deviation (median = 1, standard deviation = 7.6 × 10^{–3}; Figure 6A). Times of extinction increased in 170 cases, decreased in 15 (all when ω_{y} = 0.975 and σ = 10) and did not change in 247, thus being significantly different from when the plasmids did not interact, i.e., σ = 1 (χ^{2} test with Yates continuity correction, *df* = 2, *P* = 2.2 × 10^{–16}). When compared relatively to the condition when plasmid X is alone in the population (model 1), both σ = 1 and σ = 10 decreased the times of extinction (Figure 6B). Nonetheless, the outcomes of σ = 1 and σ = 10 are significantly different (Wilcoxon test, *P* < 2.2 × 10^{–16}, Supplementary Table S13), but the effect is very small (σ = 1, median = 9.66 × 10^{–1}, standard deviation = 8.87 × 10^{–2}; σ = 10, median = 9.67 × 10^{–1}, standard deviation = 8.83 × 10^{–2}). In conclusion, this interaction has a very limited effect on the time of extinction.

**Figure 6.** Relative time of extinction (ToE) when plasmid interactions (δ) affect loss rates. **(A)** Relative time of extinction (ToE) in presence and absence of interactions affecting loss rates: the vertical axis represents ToE of plasmid X when there are interactions (model 6) and the horizontal axis represents the respective ToE when there are no interactions (σ = 1; model 2). **(B)** Relative ToE in presence and absence of a competitor plasmid Y: the vertical axis represents ToE of plasmid X when Y is present and plasmid interactions affect loss rates (model 6) while the horizontal axis represents the respective ToE when plasmid X is alone in the bacterial population (model 1). The analysis is always in function of the focal plasmid X. The blackline represents *y* = *x*.

### Multiple Interactions

Our ultimate goal is to understand if the presence of a second plasmid can be beneficial, and which type of interactions are determinant to prevent/delay extinction. To attain this goal, we collected the results obtained with the models where both plasmids are present, including only those cases where there were interactions. From these 5184 unique cases, the focal plasmid went to extinction in 5056. We normalized their times of extinction relative to the same focal plasmid when alone (model 1), such that the plasmid benefits from a second plasmid if the ratio is greater than one. There was a total of 371 cases where there was a benefit. Then we analyzed the proportion of outcomes (benefit vs. no benefit) among the different variables, which retrieved a significant result (χ^{2} test with Yates continuity correction, *df* = 5, *P* < 2.2 × 10^{–16}, followed by Fisher pairwise tests of independence, Supplementary Table S14) showing that some interactions have stronger impacts than others (Figure 7). Epistatic interactions have the strongest effect with 23.54% cases of benefit – in these cases, the focal persists longer in the presence of the second plasmid than when there is no other plasmid in the population. The intracellular effect on plasmid X follows with 5.32% cases of benefit, thus being the second most important plasmid-interacting allowing plasmid persistence to increase. Then follows the conjugation effects on plasmid Y (which do not differ significantly), respectively, intracellular and intercellular effects with 2.43% and 3%. The intercellular effect on plasmid X and the effect on the rate of loss (not differing from each other) follow, respectively, with 0.35% and 0%, thus having the least determinant role on plasmid persistence. The important impact of epistasis is further highlighted as this was the only interaction allowing the focal plasmid to persist in 128 cases (of model 3).

**Figure 7.** Outcome of the different plasmid interactions. The vertical axis represents the ratio of time of extinction (ToE) when the plasmids interact (models 3–6) by the respective ToE of when the focal plasmid is alone in the population (model 1). The horizontal axis represents the different values of the variable annotated on the top facet. Epistasis (ε) and the remaining parameters are split in different panels to facilitate visualization of the points. Panel **(A)** represents the effect of epistasis (ε) while panel **(B)** refers respectively to the following parameters: intracellular interactions affecting the conjugative efficiency of plasmid X (α_{X}) and of plasmid Y (α_{Y}), interactions affecting loss rates (δ) and, intercellular interactions affecting the conjugative efficiency of plasmid X (ξ_{X}) and of plasmid Y (ξ_{Y}). Blue dots represent plasmid combinations where the interactions increase ToE relatively to when the plasmid was alone in the population, while yellow dots represent no change or a decrease in ToE.

To further understand the contribution of the different interactions to plasmid maintenance we simulated cases where multiple interactions may occur, thus getting all possible combinations from no interactions at all to all interactions happening simultaneously. To avoid an overwhelming number of cases, we considered a single competitor plasmid Y with ω_{Y} = 0.9 and γ_{Y} = 10^{–12}. Note that in the previous sections none of the combinations involving this set of Y plasmids could survive nor could the fitness of XY cells be higher than that of plasmid-free cells (ω_{XY} ≥ 1). Among these 23328 simulations, plasmid X could survive in 291 (1.25%) existing in both X and XY cells. These cases share the following features: ε = 0.1, α_{X} = 10, γ_{X} = 10^{–11} and ω_{X} > 0.85. This again suggests that epistasis and the intracellular conjugation effect on plasmid X are the interactions most determinant for plasmid fate. Interestingly, plasmid Y can survive in the same 291 cases, where Y cells constitute at most a proportion of 9.47% of the population. Plasmid-free cells are absent in 144 of these 291 cases – only when δ = 10^{–8}.

We performed a logistic regression on the 23037 cases of extinction considering the following interactions (as continuous variables): ε, α_{X}, α_{Y}, ξ_{X}, ξ_{Y} and σ. The independent variable comprises two categories: “increase” with 3222 cases where the time of extinction increases relatively to when the plasmid is alone in the population and “decrease/null” for the remaining 19815 cases. The results are shown in Supplementary Table S15. All the six interactions have significant effects (*P* < 0.05, 34.72% ≤ pseudo *R*^{2} ≤ 47.05%, Supplementary Table S15). As suggested above, ε and α_{X} display more determinant effects (higher deviance from the null model and higher individual pseudo *R*^{2}), while σ has the lowest effect on delaying extinction. We repeated this analysis now including the 291 cases of survival as “increase” and the result is qualitatively similar (*P* < 0.05, 38.07% ≤ pseudo *R*^{2} ≤ 50.09%, Supplementary Table S16). Lastly, we analyzed a subset of cases consisting only on the conditions that allow plasmid survival, i.e., all cases where ε = 0.1, α_{X} = 10, γ_{X} = 10^{–11} and ω_{X} > 0.85. This subset has 291 cases categorized as “survival” and 195 as “extinction.” Only α_{Y}, ξ_{X}, ξ_{Y} and σ were used as predictors for the logistic regression since ε and α_{X} display no variance in this subset. The result shows that only ξ_{X} and ξ_{Y} are significant (*P* < 0.05, 6.95% ≤ pseudo *R*^{2} ≤ 9.23%, Supplementary Table S17). Therefore, we can conclude that these four interactions alone have a weaker impact on the outcome than epistasis and intracellular interactions (ε and α_{X}).

## Discussion

It is crucial to understand the conditions favoring plasmid maintenance in natural communities, which can vary in complexity with an increasing number of hosts and plasmids occupying the same ecological habitat. This work aimed to evaluate which plasmid properties and interactions mostly affect the time of persistence in the bacterial population. Our initial results show that plasmids cannot persist in bacterial populations unless they provide a fitness advantage (strictly mutualistic) or have high conjugation rates (strictly parasitic), which supports previous works (Stewart and Levin, 1977; Simonsen, 1991). Expanding these models by allowing bacterial populations to carry an additional non-interactive plasmid reveals, as expected, that plasmid persistence decreases due to competition between plasmids. Our main aim was to examine the role of different plasmid interactions, namely affecting host fitness, conjugative efficiency and the rate of loss due to missegregation, on plasmid persistence. We conclude that although all interactions evaluated can affect plasmid persistence, there is a hierarchy such that interactions affecting fitness (epistasis) have a stronger impact on plasmid maintenance than those affecting conjugation, and lastly plasmid loss.

If plasmids display epistasis and interact simultaneously at no other level, they can be maintained stably in the population if their combined effects provide a fitness benefit to their host cells. Interestingly, in some cases, plasmid-free cells may even be absent from the final population. This, yet, was never observed when populations started only with a single (beneficial) plasmid. However, fitness alone does not ensure plasmid persistence, such that beneficial plasmids or plasmid combinations can still get extinct. Indeed, our results show that plasmid combinations with positive epistasis could only persist if associated with low fitness costs, high conjugation rates and low loss rates. Nonetheless, this situation might not be entirely realistic because the expression of conjugative traits is subject to fitness costs due consumption of energy/metabolites – resulting in a tradeoff such that plasmids may enjoy high conjugative transfers or low cost, but not both (Turner et al., 1998; Porse et al., 2016). In addition, conjugation may render plasmid-carrying cells susceptible to infection by the so-called male-specific phages (phages that only infect cells expressing sex-pili), which decreases their persistence (Jalasvuori et al., 2011; Ojala et al., 2013).

When plasmid interactions could affect multiple parameters (fitness, conjugation, loss) simultaneously, the highest level of positive epistasis could save plasmids from extinction. But, it depended on a high conjugation rate and on plasmid interactions that further increased conjugative transfer. Both epistasis and the intracellular interaction acting on the conjugative transfer of the focal plasmid X were shown to be asymmetric. This means that positive effects enhancing fitness and transfer had more determinant impacts on plasmid fate than negative effects that impaired those traits with an opposite or greater strength.

Among the conjugative interactions, α_{X} was the most critical for the survival of plasmid X, because the transfer of this plasmid increased when the competitor plasmid was present in the same cell. In such a case, the interaction α_{y} was less important as it targeted the competitor plasmid. Notwithstanding, one should consider these interactions also in the last set of simulations where the two plasmids can interact in multiple ways. In such cases, α_{X} = 10 was necessary for the survival of the focal plasmid X but also of the competitor plasmid Y. Therefore, if analyzed from the perspective of plasmid Y, the importance of α_{X} and α_{y} would be reversed. It would be preferable to increase the transfer of its competitor than its own. The reason is that this would favor the emergence of XY cells since γ_{X} > γ_{y}. Nonetheless, the intracellular interaction on conjugation would still be more determinant than the intercellular.

From our models, the four conditions required for plasmid survival were: not too high fitness cost coupled with positive epistasis and, high conjugation transfer coupled with an intracellular effect further increasing it. Indeed, positive epistasis has been pointed out as an explanation for the existence of combinations of small and big (mobile) plasmids (San Millan et al., 2014a). On the other hand, intracellular interactions between conjugative plasmids decrease transfer rates more often than not. Therefore, plasmid combinations simultaneously enjoying intracellular effects on conjugation and fitness (epistasis) might be rare. Moreover, positive epistasis may result from inhibition (rather than facilitation) of conjugative transfer, such that host fitness increases by saving resources that would be spent on the expression of the conjugative machinery (Chao et al., 2000). Studies that can reveal a linkage between epistasis and intracellular interactions on conjugation, or the lack of it, are important to understand how interference between plasmids determine their stability in bacterial populations.

As proof of concept, our models show that interactions between plasmids can enhance their stability. Nevertheless, we acknowledge some limitations. To reduce the number of simulations we considered both plasmids X and Y to have equal loss probabilities, which is unrealistic. Additionally, our simulations encompass long times but consider static features, such that the entities are immutable (neither plasmids nor bacteria change by mutation or by recombination). Therefore, we neglect compensatory mutations that could mitigate the fitness cost of one (San Millan et al., 2015) or several plasmids (Loftie-Eaton et al., 2017; Hall et al., 2020; Jordt et al., 2020) and recombination events between plasmids allowing acquisition of addiction systems that alter loss rates (Loftie-Eaton et al., 2016; Stalder et al., 2017), pleiotropic effects between mutations (Jordt et al., 2020) or even the formation of co-integrates. The latter is especially important from a clinical point of view, in particular upon a fusion of plasmids encoding different resistance mechanisms thus creating multidrug-resistance plasmids (Garcia et al., 2007; Desmet et al., 2018), or even between resistance and virulence plasmids (Dong et al., 2018), that could explain the positive correlation between the diversity of resistance and virulence genes across metagenomes (Escudeiro et al., 2019). We also only considered interactions between two plasmids, but interactions between multiple plasmids (Gama et al., 2017c) can be considered in the future. Despite the limitations of our models, they demonstrate the contribution of interactions between plasmids, and their relative roles, for plasmid survival and should pave the way for further studies with more complex models.

## Data Availability Statement

The code, datasets and analyses presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://github.com/GAMMAtri/simulations_plasmid_interactions.

## Author Contributions

JG, RZ, and FD conceived the study. JG and FD designed the experiments. JG performed the experiments. JG and FD analyzed the data. JG wrote the first draft of the manuscript, with contributions of RZ and FD. All authors contributed to manuscript revision, read and approved the submitted version.

## Funding

Fundação para a Ciência e Tecnologia supported cE3c through UID/BIA/00329/2013. UiT The Arctic University of Norway through UiT publication fund.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020.02033/full#supplementary-material

## References

Anderson, E. S. (1968). The ecology of transferable drug resistance in the enterobacteria. *Annu. Rev. Microbiol.* 22, 131–180. doi: 10.1146/annurev.mi.22.100168.001023

Arnold, J. B. (2019). *ggthemes: Extra Themes, Scales and Geoms for “ggplot2”. R Package Version 4.2.0.* Available online at: https://CRAN.R-project.org/package=ggthemes (accessed May 13, 2019).

Bergstrom, C. T., Lipsitch, M., and Levin, B. R. (2000). Natural selection, infectious transfer and the existence conditions for bacterial plasmids. *Genetics* 155, 1505–1519.

Bouma, J. E., and Lenski, R. E. (1988). Evolution of a bacteria/plasmid association. *Nature* 335, 351–352. doi: 10.1038/335351a0

Carroll, A. C., and Wong, A. (2018). Plasmid persistence: costs, benefits, and the plasmid paradox. *Can. J. Microbiol.* 64, 293–304. doi: 10.1139/cjm-2017-0609

Chao, L., Hanley, K. A., Burch, C. L., Dahlberg, C., and Turner, P. E. (2000). Kin selection and parasite evolution: Higher and lower virulence with hard and soft selection. *Q. Rev. Biol.* 75, 261–275. doi: 10.1086/393499

Dahlberg, C., and Chao, L. (2003). Amelioration of the cost of conjugative plasmid carriage in Eschericha coli K12. *Genetics* 165, 1641–1649.

De Gelder, L., Ponciano, J. M., Joyce, P., and Top, E. M. (2007). Stability of a promiscuous plasmid in different hosts: no guarantee for a long-term relationship. *Microbiology* 153, 452–463. doi: 10.1099/mic.0.2006/001784-0

De Gelder, L., Vandecasteele, F. P. J., Brown, C. J., Forney, L. J., and Top, E. M. (2005). Plasmid donor affects host range of promiscuous incp-1β plasmid pb10 in an activated-sludge microbial community. *Appl. Environ. Microbiol.* 71, 5309–5317. doi: 10.1128/AEM.71.9.5309-5317.2005

De Gelder, L., Williams, J. J., Ponciano, J. M., Sota, M., and Top, E. M. (2008). Adaptive plasmid evolution results in host-range expansion of a broad-host-range plasmid. *Genetics* 178, 2179–2190. doi: 10.1534/genetics.107.084475

Desmet, S., Nepal, S., van Dijl, J. M., Van Ranst, M., Chlebowicz, M. A., Rossen, J. W., et al. (2018). Antibiotic resistance plasmids cointegrated into a megaplasmid harboring the bla(OXA-427) Carbapenemase Gene. *Antimicrob. Agents. Ch.* 62:e01448-17. doi: 10.1128/AAC.01448-17

Di Luca, M. C., Sørum, V., Starikova, I., Kloos, J., Hülter, N., Naseer, U., et al. (2017). Low biological cost of Carbapenemase-encoding plasmids following transfer from *Klebsiella pneumoniae* to *Escherichia coli*. *J. Antimicrob. Chemother.* 72, 85–89. doi: 10.1093/jac/dkw350

Dionisio, F. (2005). Plasmids survive despite their cost and male-specific phages due to heterogeneity of bacterial populations. *Evol. Ecol. Res.* 7, 1089–1107.

Dionisio, F., Conceição, I. C., Marques, A. C. R., Fernandes, L., and Gordo, I. (2005). The evolution of a conjugative plasmid and its ability to increase bacterial fitness. *Biol. Lett.* 1, 250–252. doi: 10.1098/rsbl.2004.0275

Dionisio, F., Matic, I., Radman, M., Rodrigues, O. R., and Taddei, F. (2002). Plasmids spread very fast in heterogeneous bacterial communities. *Genetics* 162, 1525–1532.

Dionisio, F., Zilhão, R., and Gama, J. A. (2019). Interactions between plasmids and other mobile genetic elements affect their transmission and persistence. *Plasmid* 102, 29–36. doi: 10.1016/j.plasmid.2019.01.003

Dong, N., Yang, X. M., Zhang, R., Chan, E. W. C., and Chen, S. (2018). Tracking microevolution events among ST11 Carbapenemase-producing hypervirulent *Klebsiella pneumoniae* outbreak strains. *Emerg. Microbes Infec.* 7:146. doi: 10.1038/s41426-018-0146-6

Dowle, M., and Srinivasan, A. (2019). *data.table: Extension of ‘data.frame‘. R Package Version 1.12.8.* Available online at: https://CRAN.R-project.org/package=data.table (accessed July 24, 2020).

Escudeiro, P., Pothier, J., Dionisio, F., and Nogueira, T. (2019). Antibiotic resistance gene diversity and virulence gene diversity are correlated in human gut and environmental microbiomes. *mSphere* 4:e00135-19. doi: 10.1128/mSphere.00135-19

Foster, T. J. (1983). Plasmid-determined resistance to antimicrobial drugs and toxic metal ions in bacteria. *Microbiol. Rev.* 47, 361–409. doi: 10.1128/mmbr.47.3.361-409.1983

Frost, L. S., and Koraimann, G. (2010). Regulation of bacterial conjugation: balancing opportunity with adversity. *Future Microbiol.* 5, 1057–1071. doi: 10.2217/Fmb.10.70

Gama, J. A., Zilhão, R., and Dionisio, F. (2017a). Conjugation efficiency depends on intra and intercellular interactions between distinct plasmids: plasmids promote the immigration of other plasmids but repress co-colonizing plasmids. *Plasmid* 93, 6–16. doi: 10.1016/j.plasmid.2017.08.003

Gama, J. A., Zilhão, R., and Dionisio, F. (2017b). Co-resident plasmids travel together. *Plasmid* 93, 24–29. doi: 10.1016/j.plasmid.2017.08.004

Gama, J. A., Zilhão, R., and Dionisio, F. (2017c). Multiple plasmid interference - Pledging allegiance to my enemy’s enemy. *Plasmid* 93, 17–23. doi: 10.1016/j.plasmid.2017.08.002

Garcia, A., Navarro, F., Miro, E., Villa, L., Mirelis, B., Coll, P., et al. (2007). Acquisition and diffusion of bla(CTX-M-9) gene by R478-IncHI2 derivative plasmids. *Fems Microbiol. Lett.* 271, 71–77. doi: 10.1111/j.1574-6968.2007.00695.x

Garcillan-Barcia, M. P., and de la Cruz, F. (2008). Why is entry exclusion an essential feature of conjugative plasmids? *Plasmid* 60, 1–18. doi: 10.1016/j.plasmid.2008.03.002

Gordon, D. M. (1992). Rate of plasmid transfer among *Escherichia coli* strains isolated from natural populations. *Microbiology* 138, 17–21. doi: 10.1099/00221287-138-1-17

Gyles, C., and Boerlin, P. (2014). Horizontally transferred genetic elements and their role in pathogenesis of bacterial disease. *Vet. Pathol.* 51, 328–340. doi: 10.1177/0300985813511131

Hadi, H. W. (2020). Significance of plasmids in bacillus cereus sensu lato group? a review. *Plant Arch.* 20, 881–890.

Haft, R. J. F., Mittler, J. E., and Traxler, B. (2009). Competition favours reduced cost of plasmids to host bacteria. *ISME J.* 3, 761–769. doi: 10.1038/ismej.2009.22

Hall, J. P. J., Brockhurst, M. A., Dytham, C., and Harrison, E. (2017). The evolution of plasmid stability: are infectious transmission and compensatory evolution competing evolutionary trajectories? *Plasmid* 91, 90–95. doi: 10.1016/j.plasmid.2017.04.003

Hall, J. P. J., Harrison, E., Lilley, A. K., Paterson, S., Spiers, A. J., and Brockhurst, M. A. (2015). Environmentally co-occurring mercury resistance plasmids are genetically and phenotypically diverse and confer variable context-dependent fitness effects. *Environ. Microbiol.* 17, 5008–5022. doi: 10.1111/1462-2920.12901

Hall, J. P. J., Wood, A. J., Harrison, E., and Brockhurst, M. A. (2016). Source–sink plasmid transfer dynamics maintain gene mobility in soil bacterial communities. *Proc. Natl. Acad. Sci. U.S.A.* 113, 8260–8265. doi: 10.1073/pnas.1600974113

Hall, J. P. J., Wright, R. C. T., Guymer, D., Harrison, E., and Brockhurst, M. A. (2020). Extremely fast amelioration of plasmid fitness costs by multiple functionally diverse pathways. *Microbiology* 166, 56–62. doi: 10.1099/mic.0.000862

Harr, B., and Schlötterer, C. (2006). Gene expression analysis indicates extensive genotype-specific crosstalk between the conjugative F-plasmid and the E. coli chromosome. *BMC Microbiol.* 6:80. doi: 10.1186/1471-2180-6-80

Harrison, E., Dytham, C., Hall, J. P. J., Guymer, D., Spiers, A. J., Paterson, S., et al. (2016). Rapid compensatory evolution promotes the survival of conjugative plasmids. *Mob. Genet. Elements* 6:e1179074. doi: 10.1080/2159256X.2016.1179074

Harrison, E., Guymer, D., Spiers, A. J., Paterson, S., and Brockhurst, M. A. (2015). Parallel compensatory evolution stabilizes plasmids across the parasitism-mutualism continuum. *Curr. Biol.* 25, 2034–2039. doi: 10.1016/j.cub.2015.06.024

Harrison, E., Hall, J. P. J., and Brockhurst, M. A. (2018). Migration promotes plasmid stability under spatially heterogeneous positive selection. *Proc. Biol. Sci.* 285:20180324. doi: 10.1098/rspb.2018.0324

Heuer, H., and Smalla, K. (2012). Plasmids foster diversification and adaptation of bacterial populations in soil. *FEMS Microbiol. Rev.* 36, 1083–1104. doi: 10.1111/j.1574-6976.2012.00337.x

Howard, S. L., Strong, P. C. R., and Wren, B. W. (2009). “The three bears and virulence-associated plasmids in the genus yersinia,” in *Microbial Megaplasmids Microbiology Monographs*, ed. E. Schwartz (Berlin: Springer), 171–186. doi: 10.1007/978-3-540-85467-8_8

Jackman, S. (2020). *{pscl}: Classes and Methods for {R} Developed in the Political Science Computational Laboratory. R Package Version 1.5.5.* Sydney, NSW: University of Sydney.

Jalasvuori, M., Friman, V.-P., Nieminen, A., Bamford, J. K. H., and Buckling, A. (2011). Bacteriophage selection against a plasmid-encoded sex apparatus leads to the loss of antibiotic-resistance plasmids. *Biol. Lett.* 7, 902–905. doi: 10.1098/rsbl.2011.0384

Jordt, H., Stalder, T., Kosterlitz, O., Ponciano, J. M., Top, E. M., and Kerr, B. (2020). Coevolution of host–plasmid pairs facilitates the emergence of novel multidrug resistance. *Nat. Ecol. Evol.* 4, 863–869. doi: 10.1038/s41559-020-1170-1

Kamruzzaman, M., Shoma, S., Thomas, C. M., Partridge, S. R., and Iredell, J. R. (2017). Plasmid interference for curing antibiotic resistance plasmids in vivo. *PLoS One* 12:e0172913. doi: 10.1371/journal.pone.0172913

Kaper, J. B., Nataro, J. P., and Mobley, H. L. T. (2004). Pathogenic *Escherichia coli*. *Nat. Rev. Microbiol.* 2, 123–140. doi: 10.1038/nrmicro818

Kassambara, A. (2020). *ggpubr: “ggplot2” Based Publication Ready Plots. R Package Version 0.4.0.* Available online at: https://CRAN.R-project.org/package=ggpubr (accessed June 27, 2020).

Koraimann, G., and Wagner, M. A. (2014). Social behavior and decision making in bacterial conjugation. *Front. Cell Infect. Mi* 4:54. doi: 10.3389/fcimb.2014.00054

Kottara, A., Hall, J. P. J., Harrison, E., and Brockhurst, M. A. (2016). Multi-host environments select for host-generalist conjugative plasmids. *BMC Evol. Biol.* 16:70. doi: 10.1186/s12862-016-0642-z

Kottara, A., Hall, J. P. J., Harrison, E., and Brockhurst, M. A. (2018). Variable plasmid fitness effects and mobile genetic element dynamics across *Pseudomonas* species. *FEMS Microbiol. Ecol.* 94:fix172. doi: 10.1093/femsec/fix172

Lau, B. T. C., Malkus, P., and Paulsson, J. (2013). New quantitative methods for measuring plasmid loss rates reveal unexpected stability. *Plasmid* 70, 353–361. doi: 10.1016/j.plasmid.2013.07.007

Levin, B. R., Stewart, F. M., and Rice, V. A. (1979). The kinetics of conjugative plasmid transmission: fit of a simple mass action model. *Plasmid* 2, 247–260. doi: 10.1016/0147-619X(79)90043-X

Loftie-Eaton, W., Bashford, K., Quinn, H., Dong, K., Millstein, J., Hunter, S., et al. (2017). Compensatory mutations improve general permissiveness to antibiotic resistance plasmids. *Nat. Ecol. Evol.* 1, 1354–1363. doi: 10.1038/s41559-017-0243-2

Loftie-Eaton, W., Yano, H., Burleigh, S., Simmons, R. S., Hughes, J. M., Rogers, L. M., et al. (2016). Evolutionary paths that expand plasmid host-range: implications for spread of antibiotic resistance. *Mol. Biol. Evol.* 33, 885–897. doi: 10.1093/molbev/msv339

Lüdecke, D., Makowski, D., Waggoner, P., and Patil, I. (2020). *performance: Assessment of Regression Models Performance. R Package Version 0.4.6.* Available online at: https://CRAN.R-project.org/package=performance (accessed July 27).

Lundquist, P. D., and Levin, B. R. (1986). Transitory derepression and the maintenance of conjugative plasmids. *Genetics* 113, 483–497.

Mangiafico, S. (2019). *rcompanion: Functions to Support Extension Education Program Evaluation. R Package Version 2.3.7.* Available online at: https://CRAN.R-project.org/package=rcompanion (accessed February 09, 2020).

Microsoft and Weston, S. (2019). *foreach: Provides Foreach Looping Construct. R Package Version 1.4.8.* Available online at: https://CRAN.R-project.org/package=foreach (accessed March 03, 2020).

Microsoft Corporation and Weston, S. (2019). *doParallel: Foreach Parallel Adaptor for the “parallel” Package. R Package Version 1.0.15.* Available online at: https://CRAN.R-project.org/package=doParallel (accessed August 02, 2019).

Nojiri, H., Shintani, M., and Omori, T. (2004). Divergence of mobile genetic elements involved in the distribution of xenobiotic-catabolic capacity. *App. Microbiol. Biotechnol.* 64, 154–174. doi: 10.1007/s00253-003-1509-y

Nordström, K., and Austin, S. J. (1989). Mechanisms that contribute to the stable segregation of plasmids. *Annu. Rev. Genet.* 23, 37–69. doi: 10.1146/annurev.ge.23.120189.000345

Ogle, D. H., Wheeler, P., and Dinno, A. (2019). *FSA**: Fisheries Stock Analysis. R Package Version 0.8.25.* Available online at: https://github.com/droglenc/FSA (accessed June 1, 2020).

Ojala, V., Laitalainen, J., and Jalasvuori, M. (2013). Fight evolution with evolution: plasmid-dependent phages with a wide host range prevent the spread of antibiotic resistance. *Evol. Appl.* 6, 925–932. doi: 10.1111/eva.12076

Pedersen, T. L. (2020). *patchwork: The Composer of Plots.* Available online at: https://patchwork.data-imaginist.com (accessed June 22, 2020).

Peña-Miller, R., Rodríguez-González, R., MacLean, R. C., and San Millan, A. (2015). Evaluating the effect of horizontal transmission on the stability of plasmids under different selection regimes. *Mob. Genet. Elements* 5, 1–5. doi: 10.1080/2159256X.2015.1045115

Porse, A., Schonning, K., Munck, C., and Sommer, M. O. A. (2016). Survival and Evolution of a Large Multidrug Resistance Plasmid in New Clinical Bacterial Hosts. *Mol. Biol. Evol.* 33, 2860–2873. doi: 10.1093/molbev/msw163

R Core Team (2018). *R: A Language and Environment for Statistical Computing.* Vienna: R Foundation for Statistical Computing.

Radnedge, L., Davis, M. A., Youngren, B., and Austin, S. J. (1997). Plasmid maintenance functions of the large virulence plasmid of *Shigella* flexneri. *J. Bacteriol.* 179, 3670–3675. doi: 10.1128/jb.179.11.3670-3675.1997

Sagai, H., Iyobe, S., and Mitsuhashi, S. (1977). Inhibition and facilitation of transfer among *Pseudomonas*-Aeruginosa R Plasmids. *J. Bacteriol.* 131, 765–769. doi: 10.1128/jb.131.3.765-769.1977

San Millan, A. (2018). Evolution of plasmid-mediated antibiotic resistance in the clinical context. *Trends Microbiol.* 26, 978–985. doi: 10.1016/j.tim.2018.06.007

San Millan, A., Heilbron, K., and MacLean, R. C. (2014a). Positive epistasis between co-infecting plasmids promotes plasmid survival in bacterial populations. *ISME J.* 8, 601–612. doi: 10.1038/ismej.2013.182

San Millan, A., Peña-Miller, R., Toll-Riera, M., Halbert, Z. V., McLean, A. R., Cooper, B. S., et al. (2014b). Positive selection and compensatory adaptation interact to stabilize non-transmissible plasmids. *Nat. Commun.* 5:5208. doi: 10.1038/ncomms6208

San Millan, A., Toll-Riera, M., Qi, Q., Betts, A., Hopkinson, R. J., McCullagh, J., et al. (2018). Integrative analysis of fitness and metabolic effects of plasmids in *Pseudomonas aeruginosa* PAO1. *ISME J.* 12, 3014–3024. doi: 10.1038/s41396-018-0224-8

San Millan, A., Toll-Riera, M., Qi, Q., and MacLean, R. C. (2015). Interactions between horizontally acquired genes create a fitness cost in *Pseudomonas aeruginosa*. *Nat. Commun.* 6:6845. doi: 10.1038/ncomms7845

Schauberger, P., and Walker, A. (2019). *openxlsx: Read, Write and Edit xlsx Files. R Package Version 4.1.3.* Available online at: https://CRAN.R-project.org/package=openxlsx (accessed May 06, 2020).

Sezonov, G., Joseleau-Petit, D., and D’Ari, R. (2007). *Escherichia coli* physiology in luria-bertani broth. *J. Bacteriol.* 189, 8746–8749. doi: 10.1128/JB.01368-07

Sheppard, R. J., Beddis, A. E., and Barraclough, T. G. (2020). The role of hosts, plasmids and environment in determining plasmid transfer rates: a meta-analysis. *Plasmid* 108:102489. doi: 10.1016/j.plasmid.2020.102489

Silva, R. F., Mendonca, S. C. M., Carvalho, L. M., Reis, A. M., Gordo, I., Trindade, S., et al. (2011). Pervasive sign epistasis between conjugative plasmids and drug-resistance chromosomal mutations. *PLoS Genet.* 7:e1002181. doi: 10.1371/journal.pgen.1002181

Simonsen, L. (1991). The existence conditions for bacterial plasmids: theory and reality. *Microb. Ecol.* 22, 187–205. doi: 10.1007/BF02540223

Simonsen, L., Gordon, D. M., Stewart, F. M., and Levin, B. R. (1990). Estimating the rate of plasmid transfer: an end-point method. *Microbiology* 136, 2319–2325. doi: 10.1099/00221287-136-11-2319

Soetaert, K. (2009). *rootSolve: Nonlinear Root Finding, Equilibrium and Steady-State Analysis of Ordinary Differential Equations. R Package 1.6.*

Soetaert, K., Petzoldt, T., and Setzer, R. W. (2010). Solving Differential Equations in {R}: package de{S}olve. *J. Sta. Softw.* 33, 1–25. doi: 10.18637/jss.v033.i09

Sota, M., Yano, H., Hughes, J., Daughdrill, G. W., Abdo, Z., Forney, L. J., et al. (2010). Shifts in host range of a promiscuous plasmid through parallel evolution of its replication initiation protein. *ISME J.* 4, 1568–1580. doi: 10.1038/ismej.2010.72

Stalder, T., Rogers, L. M., Renfrow, C., Yano, H., Smith, Z., and Top, E. M. (2017). Emerging patterns of plasmid-host coevolution that stabilize antibiotic resistance. *Sci. Rep.* 7:4853. doi: 10.1038/s41598-017-04662-0

Starikova, I., Al-Haroni, M., Werner, G., Roberts, A. P., Sørum, V., Nielsen, K. M., et al. (2013). Fitness costs of various mobile genetic elements in Enterococcus faecium and *Enterococcus faecalis*. *J. Antimicrob. Chemother.* 68, 2755–2765. doi: 10.1093/jac/dkt270

Stevenson, C., Hall, J. P. J., Brockhurst, M. A., and Harrison, E. (2018). Plasmid stability is enhanced by higher-frequency pulses of positive selection. *Proc. Biol. Sci.* 285, 20172497. doi: 10.1098/rspb.2017.2497

Stewart, F. M., and Levin, B. R. (1977). The population biology of bacterial plasmids: a PRIORI conditions for the existence of conjugationally transmitted factors. *Genetics* 87, 209–228.

Toussaint, A., and Merlin, C. (2002). Mobile elements as a combination of functional modules. *Plasmid* 47, 26–35. doi: 10.1006/plas.2001.1552

Turner, P. E., Cooper, V. S., and Lenski, R. E. (1998). Tradeoff between horizontal and vertical modes of transmission in bacterial plasmids. *Evolution* 52, 315–329. doi: 10.2307/2411070

Yano, H., Wegrzyn, K., Loftie-Eaton, W., Johnson, J., Deckert, G. E., Rogers, L. M., et al. (2016). Evolved plasmid-host interactions reduce plasmid interference cost. *Mol. Microbiol.* 101, 743–756. doi: 10.1111/mmi.13407

Keywords: bacteria, plasmids, simulations, epistasis, facilitation, conjugation

Citation: Gama JA, Zilhão R and Dionisio F (2020) Plasmid Interactions Can Improve Plasmid Persistence in Bacterial Populations. *Front. Microbiol.* 11:2033. doi: 10.3389/fmicb.2020.02033

Received: 01 June 2020; Accepted: 31 July 2020;

Published: 31 August 2020.

Edited by:

Eva M. Top, University of Idaho, United StatesReviewed by:

Alvaro San Millan, Ramón y Cajal Institute for Health Research, SpainDavid Baltrus, University of Arizona, United States

Copyright © 2020 Gama, Zilhão and Dionisio. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: João Alves Gama, joao.gama@uit.no; Francisco Dionisio, fadionisio@fc.ul.pt; francisco.dionisio@gmail.com