Effect of genetic background on the evolution of Vancomycin-Intermediate Staphylococcus aureus (VISA)

Strains

Strains used for the evolution experiments: NRS70 (N315), NRS123 (MW2), and NRS384 (USA300-0114), were obtained from BEI resources (https://www.beiresources.org/).

Experimental evolution of vancomycin resistance

The strategy for the evolution experiment is outlined in Fig. 1. Parent strains were streaked on Brain Heart Infusion (BHI) plates from frozen stocks and single colonies used to establish shaking overnight cultures of BHI broth. We then propagated independent lines in BHI broth shaking at 250 rpm at 37 °C until cultures were turbid (OD600 > 1.5) before diluting cultures at a 1:500 dilution into fresh media. Strains were grown initially on BHI broth containing 1 µg/mL vancomycin then transferred to BHI broth with 2 µg/mL antibiotic. Once turbid, the cultures were similarly diluted to 4 µg/mL, then grown for 3 days before transfer to 8 µg/mL, and grown for a final 3 days. After growth in broth at the final antibiotic concentration, cultures were plated on BHI plates containing 2, 4, and 8 µg/mL of vancomycin. A single colony was picked from the highest concentration at 48 h growth. The colony was passaged twice on selective media and then examined by Gram stain to confirm species.

Genomic analysis

DNA extraction and library prep were performed as manufacturer’s instructions (Wizard Genomic DNA Purification Kit; Promega, Madison, WI, USA; Nextera XT DNA Library Prep Kit; Illumina, San Diego, CA, USA). Genome sequencing was performed on Illumina HiSeq and MiSeq with paired-end reads. Raw read data were deposited in the NCBI Short Read Archive under project accession number PRJNA525705.

Raw read data was processed by BBduk (v37.66) (Bushnell, 2016) to filter out adapters associated with Illumina sequencing and trim reads based on quality. Mutations were called based on comparison to reference sequences using breseq (Deatherage & Barrick, 2014) CONSENSUS mode. Reference genomes used for NRS70, NRS123, and NRS384 VISA strains were NC_002745.2 (BioProject: PRJNA224116, BioSample: SAMD00061099, Assembly: GCF_000009645.1), NC_003923.1 (BioProject: PRJNA224116, BioSample: SAMD00061104, Assembly: GCF_000011265.1), and NZ_CP027476.1 (BioProject: PRJNA224116, BioSample: SAMN07411405, Assembly: GCF_002993865.1) respectively. To compare mutations across the three strains, the protein sequences of NRS123 and NRS384 were aligned using BLAST against that of NRS70 (Camacho et al., 2009). Genes were considered orthologous to NRS70 if the amino acid sequences had at least 95% identity, 95% coverage, or an evalue = 0. To determine if large scale deletions were due to phage loss, PhageWeb was used to identify prophages within the genomes (Sousa et al., 2018).

Measuring fitness based on growth

Strains were grown in duplicate in a 96-well plate beginning at an OD < 0.1 and grown for 24 h at 37 °C with constant shaking in a Biotek Eon Microplate Spectrophotometer, with OD measurements every 10 min. To assess growth curves, OD readings were imported into R (R Core Team, 2016) and maximal growth rate (r) calculated using the growthcurver package (Sprouffske & Wagner, 2016) using the average of two biological replicates. Fitness was calculated as the ratio between the average r of each evolved strain to the average r of the parent strain. Fitness distributions between NRS70 evolved VISA strains and NRS384 evolved VISA strains were compared using a two-sided two-sample Kolmogorov-Smirnov test.

Statistical analysis

To determine if the difference in the prevalence of mutations across different VISA genes (walK, prs, rpoB, rpoC, vraS, yvqF) was significant, a Pearson’s Chi-squared test was performed, with subsequent post hoc pairwise comparisons using Fisher’s exact test and a Bonferroni multiple test correction. To analyze if the prevalence of mutations was also affected by background, binomial generalized linear models (GLMs) were fitted for each gene, and pairwise Wald tests were performed to test whether the difference in the coefficients for each background was non-zero. The statistical significance of gene correlation was assessed using a Spearman correlation with a Holm multiple testing correction. Total SNP distributions were compared using a Kruskal-Wallis test followed by pairwise Mann-Whitney U tests. Synonymous and nonsynonymous SNP distributions were compared with a Pearson’s Chi-squared test followed by Bonferroni corrected Fisher’s exact tests. The vancomycin and daptomycin MIC distributions between NRS70 VISA strains and NRS384 VISA strains were compared using two-sided two-sample Kolmogorov–Smirnov tests. To test if the linear relationship between vancomycin MIC and daptomycin MIC was dependent on background (NRS70 and NRS384), an Analysis of Covariance (ANCOVA) was performed.

To assess the effect of mutations in individual genes on vancomycin MIC, linear regression models were fitted using the lm function in R. walK, prs, rpoB, rpoC, vraS, and yvqF as binary variables (mutated or not) were used as predictors for log transformed vancomycin MICs. All mutations found in these genes were classified as nonsynonymous SNPs. The final model was chosen by backwards selection with the goal of minimizing Akaike information criterion (AIC) as in Eyre et al. (2017).

To investigate differences in the effect sizes of SNPs on vancomycin MIC from background, standard linear models were fitted with interaction terms for SNP presence and background. In these models, strain background was included as a fixed-effects predictor because the mixed-effects version (with background as a vector random effect affecting the main and interaction effects of the SNPs) was not computationally stable. To study potential pairwise interactions between SNPs affecting the MIC of vancomycin, linear mixed-effects models were fit with presence or absence of SNPs as fixed-effects predictors (main and interaction effects) and with background (NRS70, NRS384) as a scalar random effect. Mixed-effects models with vector random effects (affecting intercept and slope of certain interactions) were also tested, but these models were not stable enough for convergence. Given the sparsity of the data matrix, the predictors were restricted to SNPs present in at least 10% of the strains. Predictors included were the six universally mutated genes (walK, prs, rpoB, rpoC, vraS, yvqF) and sdrC. All SNPs with the exception of SNPs in sdrC, which were all synonymous, were nonsynonymous. The response variable (vancomycin MIC) was log transformed. These models were also fitted for another response (log transformed growth rate as a proxy of fitness). Similar to the vancomycin models, pairwise SNP interactions were assessed with mixed-effects models and SNP-background interactions with standard linear models. A Bonferroni corrected p-value (ɑ = 0.00056) was used to assess test significance. R scripts used for these analyses are available at https://github.com/crsl4/staph-visa.

Antimicrobial susceptibility testing

Antimicrobial susceptibility testing for daptomycin, vancomycin, and methicillin were performed according to Clinical Laboratory Standards Institute (CLSI) standards (CLSI, 2013). In brief, cultures were streaked from a frozen stock onto BHI agar and restreaked the next day. From the second day culture, colonies were resuspended in normal saline to a 0.5 McFarland standard. For plate-based testing, Etest for daptomycin and disc diffusion with cefoxitin, cultures were struck to create a lawn on cation-adjusted Mueller-Hinton broth (CAMHB) agar plates before application. Plates were incubated at 35 °C and read at 18 h. For broth microdilution (BMD) to determine strain vancomycin minimum inhibitory concentration (MIC), 0.5 McFarland standards were diluted 1:20 in normal saline, and 10 µL was added to 90 µL of CAMHB with the appropriate concentration of vancomycin. Plates were incubated at 35 °C without shaking and read at 24 h.

Patterns of acquisition of SNPs at VISA endpoints is influenced by genetic background

Three strains, representing major clinical lineages in the US: N315/NRS70 (Kuroda et al., 2001) (ST5/CC5), USA300-0114/NRS384 (Diep et al., 2006) (ST8/CC8) and MW2/NRS123 (Baba et al., 2002) (ST1/CC1) were experimentally evolved under selection pressure to achieve a vancomycin MIC of 4-10 μg/mL from an initial MIC of 1 μg/mL (Fig. 1). Selection of 120 of these mutants, 40 independent lines per genetic background, was achieved by sequentially raising the vancomycin concentration through serial passages of strains in BHI broth. Mutants were sequenced and compared to their isogenic parent using breseq (Deatherage & Barrick, 2014) to ascertain mutations acquired during vancomycin selection. Thirteen evolved strains are excluded from further analyses due to either low sequence quality or culture contamination. The final number of evolved VISA strains per background with high quality DNA sequences was 35 NRS70 strains, 34 NRS123 strains, and 38 NRS384 strains.

The majority of the genetic changes that occurred during the selection for vancomycin-intermediate resistance were single nucleotide polymorphisms (SNPs) and will be the focus of the rest of the analyses. Most deletions, insertions, and substitutions only occurred once (Tables S1–S6). Of note, prophage loss (StauST398-4 in NRS70, Pvl108 in NRS123, 23MRA in NRS384) was observed in four strains from each of the genetic backgrounds. Phage loss has been previously observed in vancomycin and daptomycin resistant isolates and has been associated with the induction of the SOS response by antibiotics (Maiques et al., 2006; Boyle-Vavra et al., 2011; Machado et al., 2020).

We found between 5 and 25 SNPs in evolved VISA strains compared to their respective parent strains, including non-synonymous mutations in several genes previously reported to be associated with VISA genes. As the genome sizes for NRS70, NRS123 and NRS384 are 2.81 MB, 2.82 MB, and 2.87 MB, respectively, we did not expect differences in genome size to significantly influence the number of mutations observed during the course of selection.

SNPs fell within 151 coding genes, and no gene was universally mutated in all strains, highlighting the existence of multiple pathways to intermediate vancomycin resistance in S. aureus. However, of the 151 genes, 121 genes had mutations in only one of the three backgrounds; 24 genes were mutated in two backgrounds; and only six genes (walK, prs, rpoB, rpoC, vraS, yvqF) had non-synonymous mutations in all three backgrounds. The six “universally mutated” genes have been previously implicated in VISA from sequencing of clinical isolates, but interestingly, some genes that had been previously reported did not acquire mutations frequently or in all three backgrounds. Notably, we did not find any agrR mutants, and graRS mutants were found in only NRS384 derived strains and walR (Howden, Peleg & Stinear, 2013) only in NRS70 and NRS123 derived strains.

The six universally mutated genes varied in their frequency across strains. On opposite sides of the spectrum, mutations in prs occurred in fewer than 10% of the evolved strains while mutations in walK occurred in 80–95% of the evolved VISA strains (Table 1). A Pearson’s chi-squared test (p < 2.2 × 10⁻¹⁶) and post hoc Fisher’s exact tests (all p < 5.91 × 10⁻⁰³) determined that the proportion of strains differed significantly between each of these six genes. There was also evidence that some gene mutations occurred more frequently in some backgrounds than others. By fitting binomial generalized linear models (GLMs) to the individual gene mutation distributions with the genetic background as a predictor, proportions of strains with rpoB (NRS70 vs. NRS123, p = 0.03; NRS70 vs. NRS384, p = 0.015), rpoC (NRS70 vs. NRS384, p = 0.019), vraS (NRS70 vs. NRS384, p = 0.04), and yvqF (NRS70 vs. NRS123, p = 0.035; NRS70 vs. NRS384, p = 0.028) were found to be significantly influenced by genetic background. Between NRS70 and NRS123, five other genes were mutated in both backgrounds. Similarly, between NRS123 and NRS384, six other genes were mutated. Finally, between NRS70 and NRS384, 12 other genes were mutated (Table 2). The higher proportion of other shared genes between NRS70 and NRS384 may indicate that these two strains share more similar evolutionary paths toward vancomycin resistance, but mutations in these genes were rare, and differences were not statistically significant.

We examined the co-occurrence of mutations in the six universally mutated genes (Figs. 2A–2D). Mutations in vraS and yvqF were almost mutually exclusive (coefficient −0.8, p < 0.0001), with minor differences between genetic backgrounds (Figs. 2B–2D, NRS70 and NRS384: −0.9; NRS123: −0.6). NRS70 was the only background where we detected significant associations between other genes (Fig. 2B): prs/vraS (coefficient 0.5, p = 0.011) and prs/yvqF (coefficient 0.5, p = 0.035).

To examine the distribution of amino acid changes caused by SNPs, we focused on yvqF and walK due to the high frequency of mutation in our VISA strains (Table 1). For yvqF, there were 32 SNPs in 68 strains: 3 SNPs are found in all three backgrounds, 10 SNPs are found in two, and 19 SNPs in only one. Despite walK mutations (41 SNPs across 93 strains) being more common, there are no SNPs in common between the three backgrounds, and 31 SNPs are found in only one. Furthermore, the most common SNPs within each background for these two genes were found in different amino acid residues (yvqF: NRS70/A152V, NRS123/P126S, NRS384/P174A; walK: NRS70/E236D, NRS123/I29M, NRS384/Q216E).

In population genetics, there is a common a priori assumption that synonymous SNPs have little or no functional impact compared to non-synonymous. However, we noted some genes with high rates of synonymous mutations. Strikingly, sdrC was mutated in 7 out of 38 NRS384 VISA strains. The sdr locus encodes sdrC, sdrD, sdrE which are members of the repeat-rich microbial surface components recognizing adhesive matrix molecules (MSCRAMM) family (Josefsson et al., 1998). These genes are not always conserved together (Liu et al., 2015), and their presence has been associated with bone infection (Sabat et al., 2006), resistance to host immunity (Sitkiewicz, Babiak & Hryniewicz, 2011; Askarian et al., 2017), biofilm formation (Barbu et al., 2014), and host-cell adhesion (Corrigan, Miajlovic & Foster, 2009; Cheng et al., 2009; Barbu et al., 2010; Askarian et al., 2016). Belikova et al. (2020) showed that in USA300 sdr genes undergo frequent within-genome recombination during growth in lab conditions and during infection. The high number of synonymous mutations may be the result of repeats expanding or collapsing in evolved strains, causing base-calling artefacts when mapped onto the reference genomes. Ultimately, long-read sequencing approaches are needed to deconvolute these complex structures, but there is a strong possibility that base-calling artefacts inflated synonymous mutation numbers in some strains. We also noted that there were significant differences between backgrounds in the total number of mutations that accrued during the evolution experiments (p = 1.1 × 10⁻⁵). NRS123 VISA (median 3 mutations) strains had fewer mutations than NRS70 (median 4, p = 0.004) and NRS384 (median 5, p = 2.8 × 10⁻⁵) VISA strains (Fig. 3). Furthermore, the number of synonymous SNPs relative to nonsynonymous SNPs was not the same (Table 3, Pearson chi square p = 1.8 × 10⁻⁵, NRS70 vs. NRS123 p = 0.046, NRS70 vs. NRS384 p = 0.038, NRS123 vs. NRS384 p = 1.0 × 10⁻⁷). However, when excluding sdr mutations, NRS70 and NRS123 are no longer statistically significantly different, but NRS384 remained statistically different from NRS70 and NRS123 (Table 3, Pearson chi square p = 0.00021, NRS70 vs. NRS384 p = 0.014, NRS123 vs. NRS384 p = 0.0011). Furthermore, the number of synonymous SNPs relative to nonsynonymous SNPs was not the same (Table 3, Pearson chi square p = 1.8 × 10⁻⁵, NRS70 vs. NRS123 p = 0.046, NRS70 vs. NRS384 p = 0.038, NRS123 vs. NRS384 p = 1.0 × 10⁻⁷).

Resistance phenotypes of evolved VISA strains

All evolved VISA strains at the endpoint of the iterative selection were able to grow in 8 μg/mL vancomycin BHI broth. However, to accurately measure the level of resistance, we performed BMD on NRS70 and NRS384 strains according to CLSI standards (unfortunately, all NRS123 frozen cultures were accidentally discarded while in −80 °C storage before we could test them). All endpoint NRS70 and NRS384 strains had MICs within the range to be considered VISA (Fig. 4A), but some (20 NRS70 and 28 NRS384 strains) were found to be below 8 μg/mL. The NRS70 VISA strain MICs ranged from 4 to 10 μg/mL, with a median and mode of 6 μg/mL, while the NRS384 VISA strain MICs ranged from 4 to 8 μg/mL, with a median and mode of 6 μg/mL. The difference in the MIC distributions was not statistically significant between NRS70 and NRS384 strains (p = 0.323). Mutated genes and/or mutated gene patterns associated with isolates with the highest level of vancomycin resistance were not appreciably different than those found in isolates with a lower resistance, likely indicating that individual mutations played a large role in the resistance achieved versus genes overall.

The cross-resistance of vancomycin and daptomycin in VISA strains has been well-studied (Cui et al., 2006; Nam et al., 2018). NRS70 and NRS384 endpoint strains were thus tested for development of daptomycin resistance by Etest. Parental strains had daptomycin MICs of 0.38 μg/mL. Thirty of the 35 NRS70 strains became resistant, and 33 of the 38 NRS384 strains became resistant (Fig. 4B). As observed with the vancomycin resistance genes, there was no clear mutational pattern associated with higher daptomycin resistance. The daptomycin distributions between NRS70 and NRS384 were not statistically different (p = 1), however, only NRS384 strains had a daptomycin MIC of 6 μg/mL versus 4 μg/mL for NRS70 strains. Greater levels of vancomycin resistance correlated to greater levels of daptomycin resistance (Fig. 5, p = 0.006). However, the regressions for each background were not significantly different (p = 0.822).

Collateral sensitivity with beta-lactams is a less frequently reported phenomenon in VISA strains (Sieradzki & Tomasz, 1999; Mwangi et al., 2007; Barber et al., 2014). Within our evolved NRS70 and NRS384 strains, only two NRS70 strains became methicillin sensitive. Known mutations for beta-lactam sensitivity in a MRSA background such as loss of SCCmec or inactivation of mecA were not found, and there the mechanism for collateral sensitivity in these strains was unknown.

Fitness decreases during evolution of higher resistance but is modulated by genetic background

We examined the fitness of intermediate and endpoint strains of each NRS70 and NRS384 experimentally evolved lineage using maximum growth rate in batch culture log(r) as an indicator. Isolates from later stages of selection, able to tolerate greater concentration of vancomycin, generally equated to lower fitness (Figs. 6A–6B). Looking at the fitness distributions at each selection stage (4 μg/mL, 6 μg/mL, 8 μg/mL) between NRS70 and NRS384, there was a striking difference at 4 μg/mL (p = 7.08 × 10⁻⁶) and 6 μg/mL (p = 1.86 × 10⁻⁶), which disappeared at 8 μg/mL (p = 0.273) (Figs. 7A–7C). Thus, NRS70 strains appeared to be able to tolerate lower levels of vancomycin resistance with lower fitness cost compared to NRS384 strains. Linear mixed models with SNP presence/absence as fixed effects and background as a scalar random effect were used to determine interaction effects between individual SNPs and log(r). However, no results were statistically significant after Bonferroni correction. The strongest interaction was between sdrC and rpoC whereby mutations in sdrC had a negative effect on growth rate in the presence of rpoC mutation from a previous positive effect (uncorrected p = 0.023, Fig. S1), but as noted above, the caveat is that sdrC mutations may be surrogates for more complex recombination events.

Genetic predictors of vancomycin MIC

To investigate the interactions of mutation and S. aureus strain background to vancomycin resistance, we used two modelling approaches. In the first, we used multiple linear regression models to determine which mutated genes were most predictive of vancomycin resistance. To reduce the number of parameters and for ease of comparison, only the six universally mutated genes were used in the initial models (Table 4). To make the models more meaningful, parameters were evaluated for inclusion in the final models using AIC (Eyre et al., 2017). In the final model for NRS70, walK (coefficient estimate 0.298), rpoC (0.228), vraS (0.989), and yvqF (1.058) were statistically significantly associated with vancomycin resistance, defined as log(MIC) (R² = 0.583, p = 2.15 × 10⁻⁶). In comparison, only walK (coefficient estimate 0.633) by itself was statistically significant in the final model of NRS384 (R² = 0.313, p = 9.95 × 10⁻⁴), although vraS (0.340) and yvqF (0.378) were included for improved model fit.

In the second approach, we fit linear regression models to determine if the effect sizes of SNPs on log(MIC) differed by background. No results were statistically significant after Bonferroni correction, but there was a suggestive result that rpoC may contribute positively to vancomycin MIC in NRS70 strains but not NRS384 strains (uncorrected p = 0.028, Fig. S2), which is consistent with our previous linear regression models. Linear mixed models were used to determine interaction effects between SNPs. No results were statistically significant after Bonferroni correction, but there was a suggestive result that the effect of a sdrC mutation (either increasing or decreasing vancomycin MIC) depended on walK (uncorrected p = 0.0472, Fig. S3). Overall, the models indicated that while the six universally mutated core genes were broadly responsible for vancomycin resistance in NRS70 and NRS384, they were of differential importance within each background. In addition, we did not find strong evidence for pairwise interaction between SNPs influencing MIC.