Mutators can drive the evolution of multi-resistance to antibiotics

Antibiotic combination therapies are an approach used to counter the evolution of resistance; their purported benefit is they can stop the successive emergence of independent resistance mutations in the same genome. Here, we show that bacterial populations with ‘mutators’, organisms with defects in DNA repair, readily evolve resistance to combination antibiotic treatment when there is a delay in reaching inhibitory concentrations of antibiotic—under conditions where purely wild-type populations cannot. In populations of Escherichia coli subjected to combination treatment, we detected a diverse array of acquired mutations, including multiple alleles in the canonical targets of resistance for the two drugs, as well as mutations in multi-drug efflux pumps and genes involved in DNA replication and repair. Unexpectedly, mutators not only allowed multi-resistance to evolve under combination treatment where it was favoured, but also under single-drug treatments. Using simulations, we show that the increase in mutation rate of the two canonical resistance targets is sufficient to permit multi-resistance evolution in both single-drug and combination treatments. Under both conditions, the mutator allele swept to fixation through hitch-hiking with single-drug resistance, enabling subsequent resistance mutations to emerge. Ultimately, our results suggest that mutators may hinder the utility of combination therapy when mutators are present. Additionally, by raising the rates of genetic mutation, selection for multi-resistance may have the unwanted side-effect of increasing the potential to evolve resistance to future antibiotic treatments.

written more compactly using vector notation and taking the dot product, ( , ) = ì · ì , where ì and ì each have length + 1. The reader may be familiar with binary logistic regression with = 2 outcomes, usually with = 1 defined as 'success' and = 0 as 'failure'. The function ( , ) is linked to the probability of observing outcome by taking the log of the odds-ratio of the two outcomes, i.e. the logit function.
where ì is the vector of values taken by the explanatory variables ì for the observation .
For > 2 outcomes, the multinomial logit can be thought of as computing − 1 independent logistic regression models with respect to a consistent reference level [1]. If is chosen as the reference level, 0, is defined as the 'intercept' and all other elements of ì are equal to zero. This results in The fact that =1 ( = ) = 1 allows calculating P ( = ) = 1/ 1 + −1 =1 ì · ì , which can then be used to solve the other probabilities. For any outcome , the general form of P ( = ) is thus given as which can then be used to estimate the coefficients ì through various methods. The particular method we used, Bayesian categorical regression, is described in the next section.

Fitting the statistical model and hypothesis testing using Bayesian categorical regression
In our particular analysis, represents the type of resistance observed, with five possible outcomes: 'no resistance', 'rifampicin resistance', 'nalidixic acid resistance', 'mixed resistance' and 'double resistance'. For 'mixed resistance', populations grew on selective plates containing either rifampicin or nalidixic acid, but not on plates containing both rifampicin and nalidixic acid, whereas 'double resistant' populations grew on plates containing both antibiotics. We note that these outcomes refer to detection of resistance, rather than fixation, i.e. the frequency of resistant individuals within each population has a value between > 0 and ≤ 1.
The independent variables, ì , were the initial mutator frequency, the antibiotic treatment applied, which microtitre plate the population inhabited, and position within each microtitre plate. Initial mutator frequency was treated as a categorical predictor (with levels 'zero', 'low', 'medium' or 'high', using 'zero' as the reference level), as we have no a priori expectation of a linear relationship between the proportion of mutators and the estimated coefficient. Antibiotic treatment was a categorical predictor (with levels 'no antibiotic', 'rifampicin', 'nalidixic acid', or 'combination', using 'no antibiotic' as the reference level). Initial mutator frequency and antibiotic were each treated as a fixed effect,as variance estimates from random effects variables with fewer than five levels tends to be imprecise [2]. Plate number and position within each microtitre plate were treated as random effects.
Incorporating both fixed and random effects requires fitting a 'mixed-effects model' to the data. However, mixed-effects models for categorical data are not straightforward to fit using standard frequentist inference methods. While it is possible to fit such models in principle [3], there is not, to our knowledge, a readily-available software implementation. However, mixed-effects categorical models can be fit using recently-developed tools that use Bayesian inference methods [4,5]. We provide an example below. In addition to the software availability, there are additional benefits of using a Bayesian inference approach [6].
To estimate the coefficients, ì , we fit a mixed-effects categorical model using brm() from the brms package [4,5] in R [7]. To fit a categorical model, 'family' was set to 'categorical' (with the default link 'logit'). To ensure convergence was achieved, we set max_treedepth = 15 and adapt_delta = 0.99. We ran four chains of 2000 iterations each, with 1000 burn-in iterations. To permit hypothesis testing on point estimates, samples from specified priors were drawn by setting sample_prior="yes". The choice of priors was based on preliminary data, and is described in detail in a later section. Default values were used for other settings.
To evaluate whether the interaction between initial mutator frequency and antibiotic treatment was important, we compared the full model (i.e. main effects and interaction) to a model with a main-effects only using 'Pareto-smoothed importance sampling leave-one-out cross-validation' (PSIS-LOO [8]). Incorporating interactions effects into the model did not significantly improve fit (PSIS-LOO difference in fit: −2.6 ± 5.4 S.E.), hence we present estimates from the main effects model.
An example of the function call is shown below (full R scripts are also available, see 'Data Availability' statement in the main text).

Establishing priors for the model
Priors were established through a preliminary experiment in which pure populations of either wild-type or mutator bacteria were subjected to increasing concentrations of single and combination antibiotic treatments (see Methods in the main text). For this preliminary experiment, we assayed the proportion of populations by measuring OD at 600 nm relative to a 'blank' well with no bacteria present, using a BMG POLARstar OPTIMA plate reader (BMG Labtech, Ortenberg, Germany). A population was considered to be 'alive' with OD > 0.1, otherwise 'extinct'. The data are presented in Fig A. We observed that, in the absence antibiotic treatment, all populations survived. For single antibiotic treatments, all mutator populations survived, but some wild-type populations went extinct. For the combination antibiotic treatment, 0/54 wild-type and 39/54 mutator populations survived. We use this information to first establish priors for the intercepts of the model ( Fig A).
As we observed no alive wild-type populations in the combination treatment, we therefore have some confidence that the proportion of multi-resistant outcomes is then double < 1/54. Therefore, the intercept for multi-resistance should be less than logit(1/54) ≈ −3.97. We therefore used a -distribution with mean = −5 and broad and heavy tails with standard deviation = 2.5, and degrees of freedom = 7, to reflect uncertainty. This distribution covers both double ≈ 0 in the left tail, and double = 39/54, i.e. the number of alive populations in pure mutator populations, in the right tail. Given that the probability of being single-drug resistant in the absence of antibiotic is likely greater than double resistant, intercepts for other resistance outcomes are likely covered by this distribution, and so we use the same prior for the single-resistance intercepts.
Next, we use this information to establish priors for the other coefficients. We observed that the presence of mutators increases survival in the presence of both single and combination antibiotics, hence the effect of mutators on total resistance (i.e. the sum of all resistance states observed) is likely to be positive. However, it is possible that different resistance states may have a negative relationship with the proportion of mutators, if for example, elevated mutation rates pushes double resistance to arise in the background of single-drug resistance, single-drug resistance may have a negative relationship with initial mutator frequency. Hence for the effects of initial proportion of mutators, we use a -distribution centred at = 0 with the same broad and heavy tails = 2.5, and = 7 to reflect uncertainty. This covers the observed proportion of alive populations in a purely mutator population = 39/54; we should however expect to observe fewer resistance events when the proportion of mutators is less than 1, as was the case in the selection experiment described in the main text. It also allows for the extreme possibilities of (nearly) zero or (nearly) all resistance, albeit with less weight given.
From this experiment, we have limited direct evidence for how the presence of antibiotics should affect the probability of observing resistance. On one hand, the presence of antibiotics decreased the number of alive populations (and an extinct population cannot be resistant). On the other, our growth measurements of resistant strains in the presence of antibiotics suggest a positive effect of being resistant in the presence of antibiotics, which would allow them to spread to high frequency and thus escape loss due to genetic drift. The effect of antibiotics is likely to be in the same range as for mutators (which includes 'no effect'), hence we use the same prior distribution (i.e. a -distribution with = 0, = 2.5, and = 7).

Estimated model coefficients
Estimated model coefficients for the main-effects only model are shown in Table A. These are reported as treatment contrasts (i.e. relative to the treatment of no antibiotics and no mutators). We used 95% credible intervals (95% C.I.s) for all hypotheses tested. Coefficients were estimated on the logit scale (i.e. log-odds, which can assume any value between −∞ and ∞, corresponding to proportions of outcomes of resistance state = 0 and 1 respectively). A coefficient is estimated for each possible combination of response outcome (resistance state) and the predictors (initial mutator frequency and antibiotic), which indicates how a specific combination of predictor levels (e.g. 'low' mutator frequency and 'rifampicin' treatment) influences the probability that a given resistance state is observed (e.g. 'rifampicin resistance').
Mutators generally have a positive effect on resistance, as indicated by positive treatment contrasts (Table A). However, as previously noted, although we expect an overall positive association between mutator frequency and resistance, we do not necessarily expect a straightforward relationship with each resistance state. This is because the 'mixed resistance' and 'double resistance' states follow from a single-drug resistance state. Indeed, this is the case for the 'low' mutator treatment and nalidixic acid resistance state, which is overtaken by 'mixed resistance'. Recall that the mixed resistance state contains both single-drug resistance types. Thus, to assess the effect of mutators on the presence of nalidixic acid resistance types, we can combine these treatment contrasts. We can perform Bayesian hypothesis tests on combinations of treatment contrasts using hypothesis(). The combined treatment contrast for the 'low' mutator frequency for the combination of 'nalidixic acid resistance' and 'mixed resistance' is 3.36 [95% C.I. of (1.74, 5.00)], which is positive as predicted. The same combination of treatment contrasts could be performed for rifampicin resistance or for the other initial mutator frequencies, but as these already exclude zero, the combinations of their contrasts will also exclude zero.

Checking robustness against choice of priors
A potential consequence of using informative priors in Bayesian inference is that they may influence the posterior distribution for estimated coefficients unduly if they are not chosen appropriately. However, the use of non-informative priors does not necessarily mitigate these problems [for an in-depth discussion, see 9]. To check the robustness of our model against the originally chosen priors, we refit the model using different priors. We set Student-priors on the intercept with arbitrarily chosen of -10, -20, -30, -40, with (as before) = 2.5 and = 7. The estimated coefficients resulting from using these priors are quantitatively similar to our original model.
As a second approach to evaluating robustness, we also set strongly-informative priors on each coefficient of the model. We estimated means for each coefficient using a fixed-effects categorical model (i.e. without the random effects) by maximum likelihood using the multinom() function from the nnet package [10], which refers to this type of model as 'multinomial logistic regression'. Note that assigning priors in this fashion is used here only as a diagnostic technique, and is not recommended as a basis for assigning priors more generally. For the model with these strong priors, the majority of coefficients were again similar to the model with weaker priors. The exception was for coefficients associated with the 'double resistance' outcome. Here, because there were zero double resistance events associated with the reference levels of the main effects of the model (i.e. no mutators, no antibiotics), the maximum likelihood estimated the intercept associated with this outcome to be very small (log-odds of −38.34, equivalent to an odds ratio of approximately 4.5 × 10 −16 ). The posterior distribution for the intercept was dominated by this strong prior. Consequently, the estimated coefficients for the coefficients associated with antibiotic treatment and the presence of mutators where there were double resistant outcomes observed were much larger than those estimated when using the original weakly-informative priors. However, we note that in all cases, the qualitative outcomes with respect to antibiotic treatment and the presence of mutators remain unchanged.

Growth of strains derived from fluctuation tests Defining the statistical model
Here we determined under which conditions resistant strains achieved a growth advantage (see Fig A in S3 Appendix). Growth was characterised using area under the curve (AUC) of growth curves from OD measured at 600 nm. We compare the wild-type (E. coli K-12 BW25113) with single-and double-drug resistant mutants selected in the BW25113 genetic background through fluctuation tests (see Methods in the main text). The model fitted to the data Fig 2 in the main text is a Bayesian two-way factorial model, with 'strain' and 'antibiotic' as predictor variables. We treated AUC measured in each antibiotic treatment as a multivariate response. This was on the basis that each strain was measured at several different concentrations, hence may be non-independent. We used Student's priors, as opposed to a Gaussian priors, because doing so is more robust against extreme values, i.e. values that appear to be 'outliers', but where there is no evidence of error in data collection [11].

Establishing priors for the statistical model
The intercept of this model is the mean AUC from the growth curves of OD, which is always greater than zero. To calculate empirical AUC, we used SummariseGrowth() from the growthcurver, which uses the trapezoidal rule to approximate the integral under the curve.
OD values on this BMG FLUOstar OPTIMA plate reader are typically less than 1.

Defining the statistical model
Here we determine whether initial mutator frequency had an effect on the growth of double resistant strains that evolved during selection. Decreased growth with a higher mutator frequency may occur if deleterious variation accumulated by the elevated mutation rate hitch-hikes to high frequency along with resistance. Alternatively, increased growth may be realised if the elevated mutation rate allowed the accumulation of beneficial mutations, or increased the clonal interference among resistance mutations. We fit a Bayesian multivariate Student's mixed-effects model. We treated growth (measured by AUC) in the presence (20 mg/l) and absence (0 mg/l) as a bivariate response variable, and initial mutator frequency as the sole population-level factor ('low', 'intermediate', 'high', with the growth of double resistant strains derived in the BW25113 wild-type background as the reference level). As before, we used a bivariate model to account for correlations arising from measuring the AUC of each strain multiple times in different environments and a Student's model to be robust to 'outliers' in the data. AUC at each concentration was measured over several 'replicate' experiments, which was treated as a varying factor common to both response variables.

Establishing priors for the bivariate linear model
The same experimental procedure was used to measure growth of the selection experiment strains, hence we use the same priors as for the fluctuation test-derived strains (Fig B).

Fitting the statistical model and hypothesis testing
As previously, statistical model fitting was performed using brm(). To use a Student's model, 'family' was set to 'student' (with the default link 'identity'). To permit hypothesis testing on point estimates, samples from specified priors were drawn by setting sample_prior="yes". To ensure convergence, we set max_treedepth = 15 and adapt_delta = 0.99. Default values were used for other settings. As before, we used 95% C.I.s for hypotheses. Estimated coefficients are given in Table B