Positive selection and compensatory adaptation interact to stabilize non-transmissible plasmids

Plasmids are important drivers of bacterial evolution, but it is challenging to understand how plasmids persist over the long term because plasmid carriage is costly. Classical models predict that horizontal transfer is necessary for plasmid persistence, but recent work shows that almost half of plasmids are non-transmissible. Here we use a combination of mathematical modelling and experimental evolution to investigate how a costly, non-transmissible plasmid, pNUK73, can be maintained in populations of Pseudomonas aeruginosa. Compensatory adaptation increases plasmid stability by eliminating the cost of plasmid carriage. However, positive selection for plasmid-encoded antibiotic resistance is required to maintain the plasmid by offsetting reductions in plasmid frequency due to segregational loss. Crucially, we show that compensatory adaptation and positive selection reinforce each other’s effects. Our study provides a new understanding of how plasmids persist in bacterial populations, and it helps to explain why resistance can be maintained after antibiotic use is stopped.


Data cloning convergence diagnostics plots.
Each row corresponds to a different bacterial type, from top to bottom: . In all cases we used a uniform prior and kept 5000 posterior MCMC draws after discarding burnin. Left and middle columns show parameter estimates as a function of the number of clones for ̅ ⁄ and respectively. Note how the mean values of the posterior distribution (red lines) are converging towards the maximum likelihood estimates and the standard errors (vertical lines) are getting smaller as the number of clones increases. Plots in the right column illustrate how the standardized maximum eigenvalue also converges to zero as the number of clones increases, indicating estimability of the parameters. Examples of clone nomenclature of clones in this work: 30S+_3_2; 30, clone isolated at the end of the experiment. S, it has been subjected to one step of selection pressure from day 8 to 9. +, the clone carries pNUK73. _3, population 3. _2, clone 2 within the population 3. 30-_1_1: 30, clone isolated at the end of the experiment. -, the clone is plasmid-free. _1, population 1. _1, clone 1 within population 1. 8S+_2_2: 8, clone isolated after antibiotic treatment at day 8 (technically it was isolated after day 9). +, the clone carries pNUK73. _2, population 2. _2, clone 2 within the population 2. Bacterial isolates used to determine growth kinetic parameters for each bacterial type defined in the evolutionary model. Parameter estimates were obtained by fitting growth curves of each isolate acquired by measuring optical densities at 600nm in a 96-well plate every 20 minutes using the same environmental conditions as the competition experiment (with 4 replicates, except for the parental strain that was replicated 3 times). Parameter estimates [95% confidence intervals] obtained using the MCMC algorithm with 2 × 10 7 iterations kept after burn-in. MCMC 1: priors used were uniform (0,1 ×10 -8 ) for ̅ and uniform(0,1 × 10 11 ) for ρ. MCMC 2: uses lognormal priors with mean 0 and standard deviation equal to 1 for both parameters. MCMC 3: uses a beta(1,1) prior for ̅ and a uniform(0,1 × 10 11 ) prior for ρ. MCMC 4: uses a gamma(0.001,0.001) prior for ̅ and uniform(0,1 × 10 11 ) prior for ρ.

Model parametrization
First, we will determine growth kinetic parameters for different bacterial strains isolated from the plasmid stability experiment (see Supplementary Table 2). Then, once we have obtained estimates for the parameters that characterize the growth of each strain in a single season, we will simulate a serial transfer experiment using the evolutionary model described in the manuscript. The goal of this approach is to predict the ecological dynamics of a competition experiment based on observations of how each strain grows independently. The remaining parameters of the model will be determined using values obtained from the literature (rate of point mutation) or using other experimental data (antibiotic susceptibility, drug degradation and rate of plasmid segregation).
Let us begin by denoting the bacterial density at time t with the variable and with the concentration of environmental resource. Then bacterial growth in a homogeneous environment with resource limitation can be modelled with the following equations: with initial conditions Here ρ denotes a resource conversion coefficient, ̅ represents the maximum growth rate and the cell's affinity for the resource. This model is indeed very simple, but it has been extensively used before to describe bacterial growth in batch reactor experimental systems [8,17] and will be the core of the evolutionary model used in the main text to study plasmid dynamics.
Despite the advantages of using mechanistic models, this approach also requires some caveats. For instance, it is known that ecological and biophysical models can present identifiability problems [9,13,5] and as a consequence parametrization algorithms may fail to converge or to provide reliable estimates for all parameters. In particular, microbial growth models that contain Michaelis-Menten type nonlinearities, like the one described by equations (1-2) can be nonidentifiable due to ̅ and being highly correlated [6,8]. It has been reported that the system can be structurally identifiable if we can measure the initial concentration of limiting resource [3], but not if we only measure bacterial optical density, which is the case in our experimental setup. To overcome this limitation, it has been proposed that instead of trying to estimate values for ̅ and K, the ̅ ratio should be used to assess competitive fitness between different strains [4,7], a quantity referred to in the literature as specific affinity [8].
It is important to highlight that the objective of the mathematical model presented in this paper is not to describe the internal biological processes of the cell or to assign biological interpretation to the estimated growth kinetic parameters, but to quantify the fitness cost associated with plasmid-bearing, as well as the fitness advantage of acquiring a compensatory mutation relative to the parental strain. Therefore we will use the simple Monod model described by Supplementary Equations (1-2) as an empirical metabolic model whereby the specific affinity ( ̅ ) and the resource conversion coefficient (ρ) will provide us with a measure of the relative fitness of different strains. Both parameters were jointly estimated by fitting the Monod model to growth curve data (optical density measurements) assuming normally distributed errors using a Metropolis-Hastings Markovchain Monte Carlo (MCMC) method implemented in R, a free open-source statistical package [14], with scripts available in a public repository [16].
The estimated parameter values for different prior distributions are summarized in Supplementary Table 3. In all cases, we used 2.5 × 10 7 iterations with a burn-in period of 5 × 10 6 and a thinning of 100 iterations. In order to adjust the acceptance probabilities for the proposed updates we used an adaptive proposal variance method as a first step of the MCMC algorithm. During this period no samples were stored in the chain until the acceptance ratio reached a target value, then we fixed the proposal variance and started storing the MCMC iterations. Convergence of the Markov chains was assessed by visual inspection.
The resulting posterior distributions, as well as the diagnostic plots of the MCMC algorithm for each one of the bacterial strains are presented in Figures S5-S8. Note how, although the MCMC chains appear to have converged and the estimated parameters provide a good fit to the data capturing the essential qualitative features, this simple model systematically underestimates the bacterial density during the first few hours of exponential phase. This could be a consequence of the rich media used in the experiments (LB broth) containing multiple limiting substrates [15]. While we could, of course, pose more complex models to capture these features [1,12], for the questions we address in this paper we consider the simpler evolutionary model to be more appropriate [18].
Furthermore, in order to assess the estimability of the parameters of our model and to obtain maximum likelihood estimates for each parameter, we implemented a data-cloning method as described in [10,11]. This method is based on a hypothetical situation whereby an individual performs simultaneously k independent experiments and by coincidence produces the same results. This is, of course, an unrealistic scenario, but we can simulate it by generating k clones of the original data and using a Metropolis-Hastings MCMC method modified to consider a new likelihood function given by the original function raised to the kth power. Lele and co-authors demonstrated in [10,11] that as the number of clones increases, the marginal posterior distribution converges to a multivariate normal distribution with mean equal to the maximum likelihood estimate with approximate variance corresponding to k times the posterior variance. Supplementary Fig. 9 illustrates the results obtained after applying the data cloning algorithm to our model and optical density data.
Finally, we used the ANOVA test proposed in [2] to confirm that the parameters are estimable and that the estimates obtained using the MCMC algorithm are reliable. The objective of this test is to reject the null hypothesis that there are significant differences in estimates when changing the number of clones or when considering different prior distributions. For instance, the plasmid-bearing parental strain with k = {1000,2000} clones, four different prior distributions (described in Supplementary Table 3) and 40000 posterior samples kept after burn-in produced a p-value for cloning effect of 0.7012 for ρ and 0.4025 for ̅ , and a p-value for prior effect of 0.8152 and 0.9106 respectively (analogous pvalues were also obtained for the other bacterial types). These values indicate that there are no significant cloning or prior effects and as a consequence we conclude that the MCMC chains have converged to the maximum likelihood estimates. We will use the