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

Vancomycin-intermediate Staphylococcus aureus (VISA) typically arises through accumulation of chromosomal mutations that alter cell-wall thickness and global regulatory pathways. Genome-based prediction of VISA requires understanding whether strain background influences patterns of mutation that lead to resistance. We used an iterative method to experimentally evolve three important methicillin-resistant S. aureus (MRSA) strain backgrounds—(CC1, CC5 and CC8 (USA300)) to generate a library of 120 laboratory selected VISA isolates. At the endpoint, isolates had vancomycin MICs ranging from 4 to 10 μg/mL. We detected mutations in more than 150 genes, but only six genes (already known to be associated with VISA from prior studies) were mutated in all three background strains (walK, prs, rpoB, rpoC, vraS, yvqF). We found evidence of interactions between loci (e.g., vraS and yvqF mutants were significantly negatively correlated) and rpoB, rpoC, vraS and yvqF were more frequently mutated in one of the backgrounds. Increasing vancomycin resistance was correlated with lower maximal growth rates (a proxy for fitness) regardless of background. However, CC5 VISA isolates had higher MICs with fewer rounds of selection and had lower fitness costs than the CC8 VISA isolates. Using multivariable regression, we found that genes differed in their contribution to overall MIC depending on the background. Overall, these results demonstrated that VISA evolved through mutations in a similar set of loci in all backgrounds, but the effect of mutation in common genes differed with regard to fitness and contribution to resistance in different strains.

It is vital for any future gene-based test for VISA to know if there are epistatic interactions between S. aureus strain background and drug resistance mutations. Recent publications have reported important epistatic interactions that drive patterns of antibiotic-resistance (Schubert et al., 2018;Ma et al., 2020). There are also examples of strain-specific effects in S. aureus: clinical CC30 S. aureus strains were found to have elevated (average 100 fold higher) persister formation compared to CC5, CC8, CC30, and CC45 strains (Liu et al., 2020), and under experimental evolution, CC398 strains were found to frequently evolve ciprofloxacin resistance due to a lineage-specific IS element that allowed amplification of the norA efflux pump (Papkou et al., 2020). In this study, we aimed to evaluate the relationship between genetic determinants of intermediate vancomycin resistance and strain background to determine if epistasis can direct the evolutionary trajectory of this phenotype. We evolved VISA strains from three common genetic backgrounds encountered clinically and analyzed the differences in mutation, fitness, and resistance levels.

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.

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 For each of the three S. aureus genetic backgrounds, up to 40 independent lineages were evolved from single colonies isolated from frozen culture streaks. The lineages were cultured in BHI broth containing increasing concentrations of vancomycin. At each stage, single colonies were isolated and frozen and tested for antibiotic MIC using Etest strips (vancomycin and daptomycin) and fitness using a growth curve based assays. Full details in the materials and methods section.

Statistical analysis
To determine if the difference in the prevalence of mutations across different VISA genes (walK, prs, rpoB, rpoC, vraS, yvqF)  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. evolved under selection pressure to achieve a vancomycin MIC of 4-10 mg/mL from an initial MIC of 1 mg/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.

Patterns of acquisition of SNPs at VISA endpoints is influenced by genetic background
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 −16 ) and post hoc Fisher's exact tests (all p < 5.91 × 10 −03 ) 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 −5 ). NRS123 VISA (median 3 mutations) strains had fewer mutations than NRS70 (median 4, p = 0.004) and NRS384 (median 5, p = 2.8 × 10 −5 ) VISA strains (Fig. 3). Furthermore, the number of synonymous SNPs relative to nonsynonymous SNPs was not the same (

Resistance phenotypes of evolved VISA strains
All evolved VISA strains at the endpoint of the iterative selection were able to grow in 8 mg/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 mg/mL. The NRS70 VISA strain MICs ranged from 4 to 10 mg/mL, with a median and mode of 6 mg/mL, while the NRS384 VISA strain MICs ranged from 4 to 8 mg/mL, with a median and mode of 6 mg/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 mg/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 mg/mL versus 4 mg/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 mg/mL, 6 mg/mL, 8 mg/mL) between NRS70 and NRS384, there was a striking difference at 4 mg/mL (p = 7.08 × 10 −6 ) and 6 mg/mL (p = 1.86 × 10 −6 ), which disappeared at 8 mg/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 2 = 0.583, p = 2.15 × 10 −6 ).
In comparison, only walK (coefficient estimate 0.633) by itself was statistically significant in the final model of NRS384 (R 2 = 0.313, p = 9.95 × 10 −4 ), although vraS (0.340) and yvqF (0.378) were included for improved model fit. Figure 6 Fitness of NRS70 and NRS384 evolved strains. 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 and maximal growth rate (r) calculated using the growthcurver package. Fitness was calculated as the ratio between the average r of each evolved strain to the average r of the parent strain. (a) Fitness (r) of NRS70 evolved strains plotted against the concentration of vancomycin that the strains were evolved. NRS70-7 at vancomycin selection stage 6 and 8 and NRS70-14 at vancomycin selection stage 6 were excluded due to lack of growth or too large a variance between replicates. MICs were determined by BMD. 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 and maximal growth rate (r) calculated using the growthcurver package. 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.  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.

DISCUSSION
In this project we looked for evidence of strain specificity in patterns of genetic change responsible for transition to vancomycin intermediate resistance in S. aureus. We found evidence that mutations in the same gene had differential effects on strain fitness and MIC depending on the genetic background. There were also statistically significant correlations in the patterns of co-occurrence of mutations linked to VISA. We did not probe the mechanisms of the putative epistatic interactions in this study. It is possible that, since many VISA-linked genes affect global cellular regulatory pathways, there is cross-talk with fixed mutations in the strain background that set global state changes (Priest et al., 2012). The strains also differ in their accessory gene content (Lindsay et al., 2006;McCarthy & Lindsay, 2010;Aanensen et al., 2016), which may also affect gene expression and protein interactions responsible for the VISA phenotype. This work adds to a growing body of literature on the importance of strain background in this pathogen. S. aureus lineages have been shown to be heterogeneous in several significant clinically relevant phenotypes such as toxin production, biofilm formation, and host immune resistance (King et al., 2016;Su et al., 2020). Epistatic interactions between antibiotic resistance and toxicity have also been explored (Yokoyama et al., 2018).
VISA was found to be linked to mutations (primarily SNPs) in a limited number of genes (e.g., walKR, rpoB/C, vraTSR). Six previously implicated VISA-associated genes (walK, prs, rpoB, rpoC, vraS, yvqF) (Howden et al., 2010) were mutated in our three genetic backgrounds. walK followed by yvqF were the most commonly seen mutations. For rpoB, rpoC, vraS, and yvqF, NRS70 and NRS384 mutation incidence were most different, and NRS123 mutation incidence was intermediate between the others. These results may be partially explained by the closer genetic relationship between CC1 and CC8 than CC5 and CC8 (Petit & Read, 2018). Other genes were found to have been mutated in two of the three genetic backgrounds and may have a role in vancomycin resistance or act as compensatory mutations. Several are already validated as being associated with vancomycin resistance (walR (Howden et al., 2011), vraR (Kato et al., 2010Baek et al., 2017;Asadpour & Ghazanfari, 2019), mprF (Ruzin et al., 2003;Chen et al., 2018)). NRS70 and NRS384 backgrounds had more of these genes in common than with NRS123, indicating that NRS70 and NRS384 strains may experience more similar selective pressures in the larger mutational landscape outside of the previously discussed six genes.
The overall vancomycin resistance levels achieved after experimental evolution in NRS70 and NRS384 VISA strains were similar, suggesting both backgrounds are equally capable of achieving lower and higher levels of vancomycin intermediate resistance.
Linear regression demonstrated that genetic predictors of vancomycin resistance differed between backgrounds. walK, rpoC, vraS, and yvqF were all significantly associated with vancomycin MICs in NRS70, with vraS and yvqF contributing most significantly, but only walK was significant in NRS384. As NRS384 strains carried more mutations than the other two backgrounds, this may indicate that while a subset of the six universally mutated genes (walK, prs, rpoB, rpoC, vraS, yvqF) are necessary for vancomycin resistance, other genes may play a significant role in determining the level of vancomycin resistance achieved in NRS384 VISA strains but not NRS70 strains.
We found that vancomycin-intermediate resistance imposed a fitness cost (measured by maximal growth rate compared to parent strain) to S. aureus that linearly scaled with the level of vancomycin resistance in most cases. Fitness distributions between these two genetic backgrounds were markedly different at stages 4 mg/mL and 6 mg/mL, with vancomycin resistance imposing less of a fitness cost on NRS70 VISA strains than NRS384 VISA strains. This may have clinical implications as it suggests that NRS70 S. aureus strains may evolve vancomycin resistance quickly, not be rapidly selected against due to growth defects, and may reach fixation if vancomycin concentrations are sustained and lead to the extinction of the ancestral vancomycin susceptible strain.
The laboratory evolution approach used in this study gave us the power to evolve VISA multiple times independently and achieve statistical significance using general linear models. With these data, we were able to investigate the effect of different SNPs with the strain background as a fixed effect. However, we cannot draw strong conclusions about clinical relevance as clinical VISA strains are the result of within-host evolution, and immune pressure is not present within our experimental design. Other factors that are present in human infection that were omitted in our experiment include but are not limited to nutrient limiting conditions and microenvironments. These undoubtedly reduce the mutational space S. aureus has available for adaption. In addition, the mutations in our studies were grouped by gene for simplicity, but likely not all mutations in VISA genes contribute equally to vancomycin resistance, and some may have no effect. The sheer number of mutations makes it impossible to recapitulate these mutations individually. We did not run control evolution experiments without vancomycin, so we cannot say whether some of the mutations in the most common six gene were due to random genetic drift, although it is well established that VISA strains do not evolve from VSSA in the laboratory in the absence of antibiotic selection. Another study design limitation was that fitness determinations by measuring maximal growth rate do not always correlate with competitive indices of evolved strains to ancestor strains and have been shown to be sensitive to environmental conditions. Thus, our measured fitness values may not accurately reflect within-host fitness, especially in environments with competing microbes.
The results of this study offer some encouragement that development of VISA strains may be predicted with reasonable accuracy directly from the genome sequence, an important future technique in diagnostic clinical microbiology (Su, Satola & Read, 2019). This is because the most common VISA mutations occurred in all three backgrounds, suggesting a "dictionary" of common mutations could be compiled. However, the finding that individual gene contribution to MIC may be background dependent (Table 4) suggests that prediction of the level of resistance in strains from different clonal complexes may be a much harder problem. In order to obtain the data to make exhaustive catalogs of VISA mutations and their epistatic dependencies or to build an accurate classifier, much more extensive versions of the experiment approach outlined here need to be performed, with a greater number of representative genetic backgrounds.

CONCLUSIONS
We found that there was a complex relationship between genetic background and VISA. There was clear evidence for parallel evolution of VISA through mutations in a common set of genes across strains (especially: walK, prs, rpoB, rpoC, vraS, yvqF). Some of these genes (particularly yvqF and vraS) showed negative epistasis in all strains. However, we found evidence the effects of mutations on vancomycin MIC and the relationship between fitness and MIC were dependent on strain background. These results indicate that prediction of the levels of vancomycin resistance based on genome sequence may require extensive databases of mutants from different genetic backgrounds.

Author Contributions
Michelle Su conceived and designed the experiments, performed the experiments, analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft. Michelle H. Davis performed the experiments, analyzed the data, authored or reviewed drafts of the paper, and approved the final draft. Jessica Peterson performed the experiments, analyzed the data, authored or reviewed drafts of the paper, and approved the final draft. Claudia Solis-Lemus analyzed the data, prepared figures and/or tables, authored or reviewed drafts of the paper, and approved the final draft. Sarah W. Satola analyzed the data, authored or reviewed drafts of the paper, and approved the final draft. Timothy D. Read conceived and designed the experiments, analyzed the data, authored or reviewed drafts of the paper, and approved the final draft.

DNA Deposition
The following information was supplied regarding the deposition of DNA sequences: Raw read data are available in the NCBI Short Read Archive: PRJNA525705.

Data Availability
The following information was supplied regarding data availability: The raw data are available in the Supplemental Files. Antimicrobial susceptibility testing results for VISA strains are given for vancomycin, daptomycin, and cefoxitin. These data were used for analyses regarding MIC distribution and fitness comparisons. The raw data used for fitness/growth rate (r) analyses is given: OD measurements every 10 min for 24 h.

Supplemental Information
Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.11764#supplemental-information.