A genomic portrait of the emergence, evolution, and global spread of a methicillin-resistant Staphylococcus aureus pandemic.

The widespread use of antibiotics in association with high-density clinical care has driven the emergence of drug-resistant bacteria that are adapted to thrive in hospitalized patients. Of particular concern are globally disseminated methicillin-resistant Staphylococcus aureus (MRSA) clones that cause outbreaks and epidemics associated with health care. The most rapidly spreading and tenacious health-care-associated clone in Europe currently is EMRSA-15, which was first detected in the UK in the early 1990s and subsequently spread throughout Europe and beyond. Using phylogenomic methods to analyze the genome sequences for 193 S. aureus isolates, we were able to show that the current pandemic population of EMRSA-15 descends from a health-care-associated MRSA epidemic that spread throughout England in the 1980s, which had itself previously emerged from a primarily community-associated methicillin-sensitive population. The emergence of fluoroquinolone resistance in this EMRSA-15 subclone in the English Midlands during the mid-1980s appears to have played a key role in triggering pandemic spread, and occurred shortly after the first clinical trials of this drug. Genome-based coalescence analysis estimated that the population of this subclone over the last 20 yr has grown four times faster than its progenitor. Using comparative genomic analysis we identified the molecular genetic basis of 99.8% of the antimicrobial resistance phenotypes of the isolates, highlighting the potential of pathogen genome sequencing as a diagnostic tool. We document the genetic changes associated with adaptation to the hospital environment and with increasing drug resistance over time, and how MRSA evolution likely has been influenced by country-specific drug use regimens.

Methicillin-resistant Staphylococcus aureus (MRSA) remains a major global cause of health-care-associated infections since it was first described five decades ago following the introduction of penicillinase-stable β-lactam antibiotics into clinical practice (Klein et al. 2007). During this period, MRSA strains have emerged several times by acquiring variants of staphylococcal cassette chromosome mec (SCCmec) elements that carry the mecA methicillin resistance determinant (Robinson and Enright 2003; Nübel et al. 2008). However, the vast majority of MRSA isolated worldwide belong to a limited number of clones, some of which are associated with global epidemics.

Epidemic MRSA-15 (EMRSA-15) has proven to be particularly successful over the last two decades, transmitting rapidly within and between hospitals, as well as to different countries. EMRSA-15 carries a type IV SCCmec element and belongs to multilocus sequence type ST22. Initially isolated in the southeast of England in 1991 (Richardson and Reith 1993), the spread of EMRSA-15 in the UK health-care setting was rapid; by 2000, EMRSA-15 accounted for over 60% of MRSA nosocomial bacteremias in England (Johnson et al. 2001). Over this period the proportion of MRSA among S. aureus bacteremia in the UK increased from <2% to 40% (Johnson et al. 2001).

Shortly after EMRSA-15 was reported in England, MRSA indistinguishable from EMRSA-15 were reported to cause outbreaks in other geographic regions, including Germany, the Czech Republic, Portugal, New Zealand, Australia, and Singapore (Pearman et al. 2001; Witte et al. 2001; Hsu et al. 2005; Melter et al. 2006; Amorim et al. 2007; Grundmann et al. 2010). These belong to multilocus sequence type ST22 or its close relatives. In a recent structured survey, ST22 was the most frequently found MRSA across Europe (Grundmann et al. 2010), and it was among the three most frequent MRSA clones in seven out of 15 countries for which quantitative data were available (Grundmann et al. 2010). Although reported in many countries, ST22-MRSA accounted for only 0.2% of MRSA in a surveillance study recently performed in the USA (Limbago et al. 2009) and has not, as yet, been reported from South America (Sola et al. 2008) and from large parts of Asia (Arakere et al. 2005; Ko et al. 2005). ST22-MRSA has repeatedly demonstrated its ability to supplant and replace other formerly established MRSA strains (Hsu et al. 2005; Amorim et al. 2007; Aires-de-Sousa et al. 2008; Witte et al. 2008). Of further concern, ST22-MRSA containing the lukS/F genes (encoding the Panton-Valentine leucocidin) has been reported to cause community-associated infections (Linde et al. 2005; Witte et al. 2007).

The mechanisms that lead to the emergence and spread of MRSA strains are of great interest. Despite intense research in this area, the spatiotemporal dynamics of the population structure and epidemic spread remain poorly understood, partially due to the limited resolution provided by contemporary genotyping tools. In this study we have investigated the origins and evolution of EMRSA-15 by whole-genome sequencing 193 ST22 isolates from 15 countries, collected between 1990 and 2009. Integrating evolutionary and spatial information, we have reconstructed the spatial and temporal dynamics underpinning the expansion of this clone and ascertained the genetic changes correlating with its enhanced spreading success.

Results

Phylogenomic reconstruction of the ST22 lineage

We determined full-genome sequences from 193 ST22 isolates (Supplemental Table S1) including isolates both from hospital and community settings, and also 12 isolates from the earliest description of EMRSA-15 (Richardson and Reith 1993). Once mobile genetic elements (MGEs) were excluded, a total of 8095 single nucleotide polymorphisms (SNPs) in the “core” genome were used to reconstruct the evolutionary history of S. aureus ST22 (Fig. 1). These data revealed that the core genome has accumulated variation at a rate of 1.3 × 10−6 substitutions per nucleotide site per year (95% Bayesian credible intervals, 1.2 × 10−6 to 1.4 × 10−6; Supplemental Table S2). Substitution rates varied slightly among different phylogenetic clades within ST22, but these differences were not significant (Supplemental Table S2). Our rate estimate was similar to those previously reported for other Staphylococcus aureus lineages (Lowder et al. 2009; Harris et al. 2010; Nübel et al. 2010). This estimate then provided a temporal calibration with which to date the emergence and geographic spread of EMRSA-15, and evolutionary events that have contributed to its success.

Figure 1.

Phylogeny of ST22 and the emergence of MRSA clones. (A) Maximum likelihood phylogenetic tree of ST22 isolates. The tree was rooted by using the distantly related S. aureus isolate MSSA476 as an outgroup. Colors indicate the isolates' countries of origin. Roman numerals indicate acquisitions of structurally different SCCmec elements, which cause methicillin resistance. (B) Maximum clade credibility tree of the ST22-A clade based on BEAST analysis using a variable clock rate (uncorrelated lognormal) model. Tips of the tree are constrained by isolation dates, the time scale is shown at the bottom. Gains and losses (Δ) of genetic determinants for resistance to methicillin (SCCmec IVh), fluoroquinolones (point mutations in grlA and gyrA), erythromycin (plasmid-encoded ermC), and clindamycin (mutations in ermC leader peptide region, c-ermC) have been mapped on the tree by applying the parsimony criterion.

The level of homoplasy in the SNP data set was extremely low (homoplasy index, 0.009), reflecting the clonal propagation of the core genome. The majority of isolates (162 of 193, 84%), collected from multiple countries, clustered in a single clade (ST22-A) (Fig. 1A), which encompasses the epidemic strain EMRSA-15. All but four of the ST22-A isolates carried a single nucleotide deletion within the ureC gene, resulting in a frame-shift inactivation of the urease alpha subunit and the consequent nonproduction of urease (Supplemental Table S1), which is a phenotypic characteristic of EMRSA-15 (Pearman et al. 2001). The four isolates lacking this deletion are the most basal isolates of the ST22-A clade (for a list of the genetic variation that distinguishes the ST22-A clade from the rest of the ST22 population, see Supplemental Table S3).

The origin of 173 of the isolates was recorded. Whereas the vast majority (97%) of ST22-A isolates were collected in hospitals, this proportion was significantly smaller in non-ST22 isolates (P < 0.001; Supplemental Table S1). This suggests that members of ST22-A have adapted to transmission and survival in health-care settings. In contrast, other clades within ST22 encompass methicillin-susceptible S. aureus isolates, MRSA with different types and subtypes of SCCmec elements, and isolates with genes for the Panton-Valentine leucocidin toxin. The phylogeny indicates that these isolates have distinct evolutionary origins from the health-care-associated clone ST22-A isolates (Supplemental Table S1) and represent the ST22 founding population from which ST22-A/EMRSA-15 recently emerged.

Included in the sequenced collection are 12 isolates from the original description of EMRSA-15 (Richardson and Reith 1993). At the time, Richardson and Reith noted that the majority of EMRSA-15 in the UK were ciprofloxacin sensitive; however, they identified a group of variant isolates from the English West Midlands (Birmingham region) that were ciprofloxacin resistant. All eight ciprofloxacin-sensitive EMRSA-15 isolates from the Richardson and Reith study are found in the ST22-A1 clade (Fig. 1B; Supplemental Table S1), whereas the four ciprofloxacin-resistant isolates form a basal group in clade ST22-A2. It therefore appears that ST22-A1, while spreading among hospitals in the British Isles, formed a progenitor population from which a fluoroquinolone-resistant pandemic clone (i.e., ST22-A2) emerged.

EMRSA-15 genotype

A combination of gene acquisition and core genome mutation shapes the genome of ST22. Therefore, in addition to SNP discovery based on mapping sequencing reads to a reference genome, we also determined the full genetic repertoire from each of the isolates, including MGEs, by de novo assemblies of sequencing reads. Using the phylogenetic reconstruction, we have been able to examine the evolution of ST22-A, and also pinpoint the genetic events that accompanied the emergence of EMRSA-15.

A significant component of the S. aureus genome is derived from MGEs that contribute to the accessory genome and includes SCCmec, prophages, S. aureus Pathogenicity Islands (SaPI), IS elements, transposons, and plasmids (Lindsay and Holden 2004). Comparative analysis reveals the composition and fluidity of the ST22 accessory genome (Fig. 2), and this accounts for between 4% and 14% of CDS within the ST22 genomes. Comparison of the ST22-A and non-ST22-A strains reveals that, typically, the accessory genome of ST22-A strains is larger and less variable than that of non-ST22-A strains (Fig. 2A).

Figure 2.

Accessory genome diversity in ST22. CDSs in the accessory genome of each isolate were clustered and assigned to homology groups (HGs). (A) Conservation of accessory genome HGs across the ST22 population. The isolates are ordered according to the phylogeny displayed along the left and top of the figure. The conserved HGs in the pairwise comparison are displayed as a heat map matrix and colored according to the number of conserved HGs (range 64 to 235 HGs) per pair; see figure for scale. (B) Size and composition of each isolate's accessory genome based on the number of HG CDSs. The HGs have been subdivided into different mobile genetic element types based on matches to an annotated MGE reference set included in the clustering (see legend to figure).

Analysis of the origins and location of CDSs in the ST22 accessory genomes show MGEs belonging to a range of classes (Fig. 2B), with prophages being the largest component. The accessory genome variation within the ST22-A isolates results from both the loss of ancestrally acquired elements such as prophages and SCCmec (Fig. 2) and the gain on multiple occasions of plasmids and prophages (Supplemental Table S1).

The ST22 reference genome, HO 5096 0412, is a representative of EMRSA-15 and is found in the ST22-A2 clade (Supplemental Table S1). Among the MGEs it contains are the following: an SCCmec IVh element (Milheirico et al. 2007), which carries the mecA gene encoding methicillin resistance; a Tn554 like-transposon, encoding a β-lactamase; two prophages (φSa2HO 5096 0412, which encodes a BcgI-like restriction enzyme, and φSa3HO 5096 0412, which encodes staphylokinase [Sako and Tsuchida 1983]) and the immune evasion proteins CHIPS and SCIN (van Wamel et al. 2006); and a 2.4-kb ermC carrying plasmid (Köser et al. 2012). Comparative genomic analysis of the ST22 isolates showed that these elements were conserved in the ST22-A population, with few exceptions (Supplemental Table S1).

Using the phylogeny we showed that the SCCmec IVh element was acquired once, just prior to the emergence of ST22-A. Similarly, φSa2 appears to have been acquired prior to the emergence of ST22-A (Supplemental Table S3). Prophage φSa3 seems to have been acquired before the emergence of the ST22-A lineage because it is present in members of the non-ST22-A population. A SaPI element encoding enterotoxin C (Novick 2003) is absent from the reference genome, but present in ∼68% of the ST22-A isolates (Supplemental Table S1). This element is also found in non-ST22-A isolates at the same site (3′ region of GMP synthase, SAEMRSA1503430), suggesting that like φSa3, its association with the ST22 lineage predates the emergence of ST22-A. In addition, other MGEs are variably present in the pandemic clone, several of which are associated with the transfer of antibiotic resistance genes; but, unlike the elements above, none of these have been stably maintained since ST22-A2 emerged.

In addition to the MGE-induced changes in the genomic architecture of ST22-A, all ST22-A genomes possess five indels in the core genome, including a 2268-bp deletion in the fibronectin-binding protein (FnBP) locus (Supplemental Table S3). This FnBP deletion was caused by homologous recombination between the C-terminal region of fnbA and the N-terminal region of fnbB, resulting in a gene fusion, fnbA–fnbB. FnBPA and FnBPB bind the host extracellular-matrix proteins fibronectin, fibrinogen, and elastin (Greene et al. 1995; Roche et al. 2004), allowing the S. aureus to connect to cellular integrins and trigger internalization (Fowler et al. 2000). In the ST22 isolates containing fnbA and fnbB (non-ST22-A), these two homologs are ∼71% identical (∼67% identical at the amino acid level) and 2871 bp and 3045 bp in length. The predicted recombination site occurs ∼1850 bp into fnbA, thereby generating a gene fusion that retains most of the sequence of fnbA. The fusion protein is the same length as FnBPA and is 98.5% identical, whereas it is only 68.6% identical to FnBPB. It is therefore likely that the encoded fusion protein retains most of the binding properties of FnBPA (Roche et al. 2004).

From the phylogeny we defined the genetic changes that differentiate the ST22-A2 and ST22-A1 populations (Supplemental Table S4). That differentiation was not marked by any unique acquisitions of MGEs. Instead, the variation in the core genome that marks this split includes a single insertion and 30 SNPs, of which 18 result in nonsynonymous substitutions, and one nonsense mutation. The nonsense mutation introduces a premature stop codon in the ECM-binding protein homolog Ebh. This very large surface protein binds fibronectin and contains a C-terminal transmembrane domain that anchors the intact protein in the cell envelope (Clarke et al. 2002). The truncated variant caused by the nonsense mutation lacks this transmembrane domain, and ST22-A2 isolates likely do not express this adhesin on their surface.

Demography and phylogeography

The high proportion of ST22-A isolates in the sample and their very low sequence diversity, points to ST22-A being a highly fit, recently emerged MRSA strain. This inference is further supported by coalescence analyses estimating that the population size of ST22-A has grown exponentially by 24% each year (Bayesian credibility intervals, 20%–29%) since the 1990s, which is four times faster than the average growth of other ST22 strains (non-ST22-A) (Fig. 3). The rapid population expansion is accompanied by the star-shaped branching pattern of the ST22-A clade (Fig. 1A).

Figure 3.

Exponential growth rates (λ) of ST22-A and other ST22 isolates (posterior probability density functions from Bayesian analysis). ST22-A is estimated to have grown significantly faster than non-ST22-A. The population size N at time t was modeled by N(t) = N(0) × e(λ × , where λ is the exponential growth rate and t is measured in years. Accordingly, the population size change per year is eλ.

A Bayesian phylogeographic reconstruction of ancestral nodes in the phylogenetic tree supports a root of the ST22-A clade in the UK. Within a few years of its emergence there, multiple exports spread it from the UK to Australia as well as to various countries in Europe, Asia, and Africa (Fig. 1B; Supplemental Fig. S1). Early ST22-A isolates (clade ST22-A1 in Fig. 1B) were susceptible to fluoroquinolones, and apparently did not spread beyond the British Isles. Bayesian modeling of the spread of fluoroquinolone-resistant ST22-A2 in the UK suggests that ST22-A2 originated within the English Midlands (Fig. 4), and was restricted to this geographical region until the late 1980s. The same analysis further suggests that isolates from the Midlands seeded infection to the north, reaching Scotland by around 1993, but also to the south, reaching London around 1998 (Fig. 4; Supplemental Movie S1). By 2000, the epidemic had spread to most locations within the UK via multiple routes and sources (e.g., Bristol, from the Southern England and Midlands; Manchester, from the Midlands and the North; Aberdeen from elsewhere in Scotland and North England). The UK epidemic then spread globally through multiple international transmission events (Fig. 1).

Figure 4.

Bayesian reconstruction of the spread of ST22-A2 in the UK. A continuous spatial diffusion model was used to reconstruct the finer-scale geographical dispersal of ST22-A2 within the UK and to predict the origin of fluoroquinolone resistance. Lines indicate the inferred routes of spread with confidence displayed as green ovals representing 80% of the highest posterior density (HPD) for latitude and longitude. The timing of transmission events are represented by red (old) to black (recent) lines and light- to dark-green oval shading (for the animation of the reconstruction of the spread, see Supplemental File S2). The maps are based on satellite pictures made available at Google Earth (http://earth.google.com).

Estimates for the timings of first appearances of ST22-A in countries beyond the UK (Supplemental Fig. S1) are concordant with available epidemiological evidence for EMRSA-15, suggesting epidemic expansion in those countries soon after its introduction (Richardson and Reith 1993; Pearman et al. 2001; Witte et al. 2001; Hsu et al. 2005; Melter et al. 2006). Our data reveal multiple introductions from the UK into Germany, Denmark, and Australia, possibly via the migration of colonized health-care workers (Pearman et al. 2001). Country-specific clades in the phylogenetic tree point to subsequent diversification within many recipient countries, and the data support additional spreading events, e.g., from Portugal to Germany (Fig. 1B). Because basal isolates outside the ST22-A clade are also geographically widespread, ST22 was almost certainly present in many countries prior to the emergence of ST22-A (Fig. 1A; Supplemental Table S1). Indeed, methicillin-susceptible ST22 is known to form a common component of the nasal carriage population in many places (Feil et al. 2003; Holtfreter et al. 2007; Ruimy et al. 2009).

The development of drug resistance in ST22

The broad range of resistance traits within ST22-A contrasts strikingly with the more susceptible and diverse population from which it emerged (Supplemental Table S1). For ST22-A, the number of resistance traits per isolate increased with time (P < 0.001) (Fig. 5), whereas the trend in the rest of the population (non-ST22-A) was not statistically significant (P > 0.1).

Figure 5.

Increase in antibiotic resistance traits within ST22-A population over time. Number of genotypic resistance determinants per ST22-A isolate (red, no. 162) and non-ST22-A isolates (blue, no. 31) from all countries, 1990–2008, including regression lines with 95% confidence interval for each group. Size of circle corresponds to number of isolates.

In order to investigate the evolutionary events responsible for the increase in resistance, we analyzed the whole-genome sequences to elucidate the antibiotic resistance genotype of ST22. Molecular determinants within the genomes explained 99.8% of the measured phenotypic resistance traits (847 resistance traits in total). Only two phenotypes did not correlate with the predicted genotypes. HO 7258 0475 05 was erythromycin resistant, but lacked ermC and all other erythromycin resistance mechanisms that have been documented in the literature. Strain 09-00678 contained thyA and dfrA, but was not phenotypically cotrimoxazole resistant. For both isolates it is possible that the drug-resistance determinants were lost during lab culturing prior to DNA preparation or phenotypic testing.

Our genomic analyses showed that evolution of resistance in ST22 has occurred in a stepwise fashion encompassing both the core and accessory genomes. Acquisition of a transposon-encoded β-lactamase BlaZ confers penicillin resistance, and the widespread distribution of this transposon in ST22 is consistent with a single, primordial development of β-lactam resistance, and probably reflects usage of penicillin since the 1940s. Only six ST22 isolates in our collection are sensitive to penicillin, four of which have lost the Tn552-like transposon carrying the blaZ gene, and one of which carries a blaZ gene inactivated through an IS element insertion in its promoter region (Supplemental Table S1). The penicillin-resistant ST22 population then gave rise to several methicillin-resistant clones, which emerged through imports of SCCmec elements into their chromosomes. The phylogenomic distribution of different SCCmec variants in the phylogenomic context provides evidence of at least seven independent imports into ST22, one replacement of SCCmec type IV by SCCmec type V, and two secondary losses of SCCmec through precise excision by homologous recombination between flanking direct repeat sequences (Fig. 1; Supplemental Table S1). As described above, ST22-A emerged after the acquisition of SCCmec IVh, which is predicted to have occurred over a decade before the first report of EMRSA-15 (Fig. 1; Richardson and Reith 1993).

Around 1986 (95% Bayesian credible intervals 1984–1988) ST22-A developed fluoroquinolone resistance through two point mutations generating amino acid substitutions Ser80Phe in topoisomerase IV (GrlA) and Ser84Leu in gyrase A (GyrA). This genotype is found in the ST22-A2 clade (Fig. 1B) that spread throughout the UK and globally. We note that our dating estimate for acquisition of fluoroquinolone resistance coincides closely with the introduction of fluoroquinolone drugs into clinical medicine in the UK, ciprofloxacin being first licensed in February 1987. Further, phylogeographic modeling suggested that fluoroquinolone resistance emerged in the English Midlands in the early 1990s (Fig. 4). Notably, clinical studies with ciprofloxacin occurred in the Midlands prior to the licensing of ciprofloxacin in 1987 (Crump et al. 1983; Finch et al. 1986; Silverman et al. 1986).

In the subsequent decades the pandemic has continued along a trajectory of increasing antibiotic resistance, developing resistance to a range of antibiotics of different classes (tetracyclines, aminoglycosides, cotrimoxazole, fusidic acid, mupirocin, and oxazolidinones) on multiple occasions (Supplemental Fig. S2; Supplemental Table S1). The development of these resistances has been rather sporadic and also less stably maintained than those against β-lactams and fluoroquinolones, but nevertheless demonstrates clear evidence of the strong selective pressure exerted by antibiotic usage and clinical practice. For example, our data set documents multiple independent events of erythromycin-resistance acquisition based on the uptake of plasmids encoding the ErmC erythromycin ribosomal methylase (Fig. 1B). Phylogenetic analysis also indicates that the ermC plasmid was lost multiple times, and in several cases these losses affect entire clades, suggesting that they represent genuine events in the natural population rather than laboratory artefacts (Fig. 1B). ErmC additionally confers resistance to clindamycin when mutations in the leader peptide region of the ermC gene cause its constitutive expression (Levin et al. 2005). All clindamycin-resistant isolates possess genetic rearrangements that disrupt the promoter region upstream of ermC, including deletions of part or all of the ermC leader peptide region or insertion of an IS element in ermC leader peptide region. At least eight independent, structurally different, mutations have generated clindamycin resistance in the population examined (Supplemental Table S1) and five of these were found in ST22-A (denoted “c-ermC” in Fig. 1B). The data also suggest that country-specific variation in antibiotic usage may be shaping the geographic distribution of resistance. The prescription level of clindamycin is much higher in Germany than in the UK (Coenen et al. 2006; Meyer et al. 2006; de With et al. 2011), causing selective pressure for a high prevalence of clindamycin resistance. Among the 34 erythromycin-resistant isolates from Germany, 32 (94%) are also clindamycin resistant, whereas none of the 45 UK isolates are resistant to clindamycin (significant difference, P < 0.001, Fisher's exact test).

Discussion

EMRSA-15 represents a “measurably evolving population” (Drummond et al. 2003), accumulating genetic variation over epidemiological timescales at an average rate similar to previous estimates for other S. aureus strains (Lowder et al. 2009; Harris et al. 2010; Nübel et al. 2010). One consequence of this fast nucleotide substitution rate is that evolutionary events could be dated, including the stepwise acquisition, refinement, and loss of multiple drug-resistance traits. Furthermore, past temporal changes in population size and geographic expansion could be analyzed retrospectively based on the molecular traces that these processes left in the genomes. These conclusions were possible because our genome sequences span a wide range of time and space, which is a prerequisite for such coalescence-based analyses. Our data revealed much more detail about the pathogen's population history than would have been possible with a snapshot survey from a single time point.

In principle, changes in the frequency of a genetic variant may be due to stochastic events or genetic adaptation and selection (Keim and Wagner 2009). For the highly successful MRSA strains described here, it seems plausible that they are more specialized to the restricted ecological niches in the hospital environment, where exposure to antibiotics is common and survival depends on efficient transmission between patients in close contact, than is the ancestral population from which they emerged. This scenario is supported by the observation that hospital-associated MRSA rarely establish themselves in the community. Aside from the development of antibiotic resistance, little is known about other adaptive genetic changes that promote the success of S. aureus (Gomes et al. 2005; Chambers and Deleo 2009). A recent microarray-based investigation did not reveal any concordant differences between gene contents of highly successful (“epidemic”) versus sporadic S. aureus strains (Kuhn et al. 2010). Recently, a phage-encoded virulence determinant (which is not present in ST22) has been shown to be associated with enhanced spreading of MRSA in Asia (Li et al. 2012). Phenotypic changes may be caused by genetic changes that are subtler than the gain of entire genes through horizontal gene transfer, or their loss (Kennedy et al. 2008; Olsen et al. 2010). However, our analysis of the genetic changes that accompanied the emergence of EMRSA-15 failed to pinpoint a single determinant, but, in contrast, identified a series of gene acquisition events and mutations that potentially promote the success of this clone in the hospital setting (Supplemental Tables S3, S4). Prominent among these are those associated with antibiotic resistance.

Evidence of the selective pressure exerted by the widespread use of broad-spectrum antibiotics in hospitals and beyond is manifest in the ST22 isolates' genomes. In order to unravel the evolutionary events surrounding the development of resistance, we examined the genotypic basis of the antibiotic resistance phenotypes. These analyses permit an evaluation of the utility of genome sequencing as a diagnostic tool for the prediction of antibiotic resistance. Our in silico analysis of genome sequences identified the molecular correlates for >99% of the phenotypic antimicrobial resistance traits. This strong correlation documents the impressive success of past research into resistance against established antibiotics, leaving few undiscovered resistance mechanisms in MRSA. Of note, several mutations in the ermC promoter region were structurally related, but not identical to previously reported mutations (Levin et al. 2005), and oxazilidinone resistance in our isolate collection was due to a mechanism that had been recently described from in vitro experiments (Locke et al. 2009), but had not previously been reported from clinical isolates (Jones et al. 2009). These observations highlight the potential of large-scale sequence analyses for discovering novel mechanisms. While various molecular methods were used in the past to detect resistance determinants in clinical isolates (Bergeron and Ouellette 1998), whole-genome sequences enabled a screen for all known genes and mutations in a single experiment. For the 12 antibiotics, we screened for 116 different alleles that were known to confer antibiotic resistance by searching antibiotic resistance genes in the de novo assemblies or SNPs in the core genome (Supplemental Table S5). Accordingly, we were able to demonstrate that it should be possible to predict with high accuracy the full resistance profile from the sequence data alone. Such predictions may be more challenging for other organisms whose mechanisms of resistances are poorly studied, or more complex, such as those associated with efflux and permeability.

The reduction in costs and time to generate genomic sequences has recently led to genome sequencing being developed as a clinical diagnostic tool that can be used to investigate and stop outbreaks in hospitals (Eyre et al. 2012; Harris et al. 2012; Köser et al. 2012; Snitkin et al. 2012). The observations from this study herald an important development in the progress of transitioning whole-genome sequencing from a research tool to a diagnostic application. Whole-genome sequencing can provide both high discriminatory power for molecular epidemiology and a prediction of antibiotic resistance, which is valuable for “near real-time” infection control and clinical management, and also for surveillance, horizon scanning of the clinical population for the evolution of drug resistance, and the identification of emerging threats.

With the benefit of hindsight it is clear that the emergence of the ST22-A2 variant of the EMRSA-15 population in the UK in the 1980s marked one such threat. Originating in the English Midlands, it spread rapidly during the 1990s, becoming endemic in UK hospitals, and was associated with an upsurge in the number of MRSA infections (Johnson et al. 2001).

Of the genetic events that define the emergence of ST22-A2, the two nonsynonymous SNPs associated with fluoroquinolone resistance would appear to be good candidates to explain the success of the clone. Fluoroquinolones are readily excreted in sweat and can suppress the growth of the normal microbiota, which includes S. aureus that colonizes skin, nose, and throat (Hawkey 1997). Fluoroquinolone resistant ST22-A2, therefore, may have a competitive edge, promoting colonization and survival in the hospital environment where this class of antibiotic is frequently used. Indeed, widespread fluoroquinolone usage has previously been reported to promote the spread of fluoroquinolone-resistant MRSA, even though this class of antibiotics is not normally used to treat S. aureus infections (Weber et al. 2003; Salangsang et al. 2010). A recent study examining MRSA rates in a tertiary hospital over a 10-yr period identified a significant increase in MRSA infection rates after an increase in fluoroquinolone prescription levels following a period of restricted usage (Parienti et al. 2011). Conversely, significant falls in MRSA rates have been linked to a decrease in the usage of fluoroquinolones over a 4-yr period (Lafaurie et al. 2012). In a prevalence study looking at elderly patients admitted to hospitals in the UK, a strong association was found between MRSA colonization and prior ciprofloxacin exposure (Hori et al. 2002), with EMRSA-15 being the dominant lineage detected in this study (75.9% of MRSA isolates).

The biological and ecological advantage of fluoroquinolone resistance would therefore appear to be clear. However, to understand the remarkable success of ST22-A2 requires insights into the historic clinical context in which it emerged. Our analysis suggests that development of fluoroquinolone resistance coincided closely with the introduction of fluoroquinolone drugs into routine clinical medicine in 1987. Initially, ST22-A2 isolates were restricted to the English Midlands until 1990, after which point the clone spread rapidly. During this time the use of fluoroquinolones in the hospital and community settings increased rapidly, doubling from 1990 to 1993 (Livermore et al. 2002). Throughout the rest of the 1990s, usage in the community stabilized, but the usage in hospitals continued to grow steadily both in absolute terms and also as a proportion of all antibiotic usage, constituting 31.5% of total usage in 1999 compared with 18.9% in 1992 (Livermore et al. 2002). The emergence and rapid spread of ST22-A2 in the UK in the 1990s, therefore, perhaps reflects changes in clinical prescription regimens. In this regard it is worth noting that we also uncovered evidence of the recent influence that country-specific antibiotic regimens (clindamycin usage in Germany vs. the UK) have had in shaping the MRSA population within Europe, illustrating the link between pathogen population structure and national antibiotic prescription policymaking.

The pandemic success of EMRSA-15 cannot, however, be simply characterized as a case of the right mutation at the right time. Fluoroquinolone resistance is not unique to ST22-A2, but is also found among non-ST22-A isolates in the collection as the result of homoplasic substitutions (Supplemental Table S1), as well as in other epidemic MRSA clones (Fluit et al. 2001; Woodford and Livermore 2009; Tenover et al. 2011). Hence, other facets of the ST22-A2 must have contributed to its outstanding success. Fluoroquinolone resistant ST22-A2 emerged from a hospital-adapted genetic background (ST22-A1) (Fig. 1) that was already successful, and spread through English hospitals by the start of the 1990s (Richardson and Reith 1993).

An important class of S. aureus proteins involved in host colonization are surface-expressed binding proteins (Edwards et al. 2012). Analysis of the genetic events that define the emergence of ST22-A identified mutations in two surface-expressed binding proteins, FnBP and Ehb, both of which bind fibronectin, a glycoprotein component of the host's extracellular matrix. Ebh is one of the largest bacterial proteins (∼1.1 MDa) (Clarke et al. 2002). The pandemic clone of EMRSA-15 has a nonsense mutation that ablates expression of Ehb at the cell envelope. Previous studies demonstrated that mutation of ebh had no effect on the growth rate or pathogenicity in a murine skin abscess model of infection (Clarke et al. 2002). Notably, this protein is intact in some successful MRSA lineages (Gill et al. 2005; Diep et al. 2006; Holden et al. 2010) and mutated in others (Kuroda et al. 2001; Holden et al. 2004). Therefore, it is unlikely that mutation of Ehb is the sole determinant of increased transmissibility of ST22-A. Its effect may be more pronounced due to other mutations, for example, the deletion in the fibronectin-binding protein locus.

Recombination between the C-terminal regions of fnbA and fnbB results in an fnbA–fnbB gene fusion and an effective reduction in fnb gene copy number at this locus due to the deletion of fnbB. Studies with clinical isolates suggested that strains associated with invasive disease are significantly more likely to have two fnb genes (Peacock et al. 2000). Molecular studies have also shown that both FnBPA and FnBPB are required for the establishment of septic infection (Shinji et al. 2011). However, in S. aureus that have both homologs, deletion of either gene does not abolish fibronectin-binding capacity of the cell (Greene et al. 1995). The functional consequences of the deletion and resulting gene fusion on the host–cell interaction and transmissibility are unclear, especially when considering that there is also a mutation in ebh.

Interestingly, although most other representatives of the main MRSA lineages contain two copies of the fnb genes (Kuroda et al. 2001; Gill et al. 2005; Diep et al. 2006; Holden et al. 2010), EMRSA-16, an unrelated epidemic MRSA lineage that originated in the UK around the same time as EMRSA-15, also contains a deletion-mediated fibronectin-binding protein gene fusion (Holden et al. 2004). Like EMRSA-15, EMRSA-16 also has evolved fluoroquinolone resistance via gyrA and grlA mutations just prior to its emergence (McAdam et al. 2012). Both of these epidemic clones spread throughout UK hospitals since the 1990s (Johnson et al. 2001), albeit EMRSA-15 has been the more successful of the two clones (Ellington et al. 2010). It is therefore intriguing that they share the fluoroquinolone resistance and fnbA–fnbB genotypes that affect resistance and host cell interaction. However, pinpointing the genetic basis of their success is more complex. In the EMRSA-16 background, mutations in the accessory gene regulator C and α-hemolysin (DeLeo et al. 2011) and squalene desaturase (McAdam et al. 2012) have been identified and hypothesized to contribute to its epidemic success. However, these are not found in EMRSA-15. In order to fully elucidate the genetic basis for the relative success of EMRSA-15, extensive experiments investigating the colonization and transmission properties of the hospital-adapted clone are required.

Methods

Bacterial isolates

Sources and properties of 193 S. aureus isolates are listed in Supplemental Table S1. In assembling the collection in this study, we included isolates that encompass wide geographical and temporal ranges, but also included isolates that reflect the documented history and prevalence of EMRSA-15 and ST22.

Susceptibility to antibiotics was tested by using the standardized broth microdilution method according to DIN58940 instructions (Deutsches Institut fuer Normung 2007), and phenotypic resistance was defined by applying minimum inhibitory concentration breakpoints from the European Committee on Antimicrobial Susceptibility Testing (EUCAST; http://www.eucast.org/). In total, 18 antibiotics were tested (penicillin, PEN; oxacillin, OXA; gentamycin, GEN; linezolid, LNZ; erythromycin, ERY; clindamycin, CLI; ciprofloxacin, CIP; fusidic acid, FUS; mupirocin, MUP; moxifloxacin, MFL; co-trimoxazole, SXT; tetracyline, TET; vancomycin, VAN; teicoplanin,TPL; rifampicin, RAM; phosphomycin, PHO; tigecycline, TGC; daptomycin, DAP), of which 12 had resistance in at least one isolate in the collection (Supplemental Table S1). Bacterial typing was performed as described previously (Strommenger et al. 2008).

Genomic library creation and multiplex sequencing

Unique index-tagged libraries for each sample were created, and up to 12 separate libraries were sequenced in each of eight channels in Illumina Genome Analyzer GAII cells with 54-base paired-end reads of Illumina HiSeq with 75-base paired-end reads. The index-tag sequence information was used for downstream processing to assign reads to the individual samples.

Detection of SNPs in the core genome

The paired-end reads were mapped against the chromosome of S. aureus HO 5096 0412 (accession no. HE681097) (Köser et al. 2012), a representative of EMRSA-15 that was isolated from a fatal neonatal infection in Suffolk, UK in February 2005. SNPs were identified as described in Croucher et al. (2011). Indels were identified using Dindel (Albers et al. 2011). Unmapped reads and sequences that were not present in all genomes were not considered as part of the core genome, and therefore SNPs from these regions were not included in the analysis. SNPs falling within MGEs regions were also excluded from the core genome, as well as those falling in high-density SNP regions, which could have arisen by recombination. The core genome was curated manually to ensure a high-quality data set for subsequent phylogenetic analysis and comprised of 2,643,131 bp.

Phylogenetic analysis

A phylogeny was drawn for ST22 using RAxML v0.7.4 (Stamatakis et al. 2005) to estimate the trees for all SNPs called from the core genome. The general time-reversible model with gamma correction was used for among-site rate variation for 10 initial trees.

Comparative analysis

Raw Illumina data were split to generate paired-end reads and assembled using a de novo genome assembly program, Velvet v0.7.03 (Zerbino and Birney 2008), to generate a multicontig draft genome for each of 193 S. aureus isolates. Gene prediction was performed using Prodigal (Hyatt et al. 2010). Predicted CDSs were compared using BLASTP (Ausubel et al. 1995), and the proteins were clustered into homology groups using TribeMCL (Center for Mathematics and Computer Science and EMBL-EBI) (Enright et al. 2002) with a cut-off of 1e−50.

Included in the clustering data set were the complete genome sequences of S. aureus available in the public sequence databases. MGEs in these sequences were identified using comparative genomic analysis and manual curation. Homology groups (HG) were categorized according to their origins based on matches to CDSs from the MGE reference set that clustered within the group. For the purpose of the preliminary accessory genome analysis, HG origins were subdivided into Prophage, Conjugative transposon, Plasmid, S. aureus Pathogenicity Island (SaPI), SCCmec element, Transposon, Hypothetical proteins, and others. Conjugative transposon/Plasmid reference sequences included Tn916-like (Clewell et al. 1985) and Tn5801-like elements (Kuroda et al. 2001) as well as integrated conjugative elements such as ICE6013 (Smyth and Robinson 2009), whereas Transposon included Tn554 (Phillips and Novick 1979), Tn552 (Rowland and Dyke 1989), and IS elements. Hypothetical protein HGs did not have any matches with either core or MGE CDSs from the reference set. HGs defined as “Other” indicate that there were no matches with MGE CDSs, but there were matches with CDSs from the core regions of the reference genomes. In some cases, there were examples of clustered proteins from more than one origin being present in the single HG; for example, Prophage and SaPI and Conjugative transposons and Plasmid. This reflects the shared evolutionary origins and mechanisms of transfer MGEs. In these cases, HGs were assigned arbitrarily.

Detailed comparisons of genome sequences were conducted on the de novo assemblies using BLASTN (Ausubel et al. 1995), and were facilitated by using the Artemis Comparison Tool (ACT) (Carver et al. 2005).

Determination of antibiotic resistance genotype

The literature was extensively mined for the genetic basis of antibiotic resistance in S. aureus. A database was compiled that contained horizontally acquired genes and core gene mutations that conferred antibiotic resistance. Antibiotic resistance conferred by SNPs in components of the core chromosome were identified by manual curation; substitutions were compared with antibiotic resistance mutations in S. aureus proteins detailed in the literature. The presence of antibiotic resistance genes of MGEs was investigated using BLASTN (Altschul et al. 1990) to compare de novo assemblies against a database of S. aureus resistance genes compiled from literature and public sequence database searches, and also by mapping sequence reads to a pseudo-molecule consisting of concatenated antibiotic resistance genes as detailed by Köser et al. (2012).

The genotype of each isolate was searched for the resistance genes and/or mutations for the antibiotics tested; for the 12 antibiotics for which resistance was detected, 847 incidences in total, 116 different genes or alleles were searched for. The presence of the following known resistance genes and mutations/deletions of determinants conferring resistance were identified: mecA (Utsui and Yokota 1985), ermC (Catchpole et al. 1988), blaZ (Voladri and Kernodle 1998), tetK (Khan and Novick 1983), fusC (O'Neill et al. 2007), fusA (mutation) (Chen et al. 2010), dfrA (Rouch et al. 1989), thyA (deletion) (Leelaporn et al. 1994), ileS-2 (Hodgson et al. 1994), aacA-aphD (Rouch et al. 1987), grlA (mutation) (Griggs et al. 2003), gyrA (mutation) (Griggs et al. 2003), rplC (mutation) (Locke et al. 2009), and ermC plasmid (deletion) (Weisblum 1995) (see Supplemental Table S1).

Assessing the association of the number of genetic determinants of resistance with time

The total number of resistance genes and mutations/deletions conferring resistance were summed for each strain. Linear regression was used to assess the association of the number of genetic determinants of resistance per strain with time. The dependent variable was the number of genetic determinants of resistance per strain, and the independent variable was time (year of isolation). Significance was assessed at P < 0.05. This analysis was applied to all 162 ST22-A isolates, from all countries, and separately to the 77 ST22-A isolates from the UK. This latter analysis was performed to evaluate the association of time with the number of genetic resistance determinants in a single country. Linear regression with the non-ST22-A strains demonstrated that the number of genetic determinants of resistance per strain was not significantly associated with time (P > 0.10).

Assessing the association of ST22-A strains with hospital origin

Out of the 162 global ST22-A strains, there were 147 where it was known whether or not they were hospital associated: 142 (96%) were hospital associated, and five were not. Of the 31 non-ST22-A strains, there were 26 where it was known whether or not they were hospital associated: 13 (50%) were hospital associated, and 13 were not. These proportions were significantly different according to Fisher's exact test (P < 0.05).

Bayesian analyses

We used the Bayesian software package, BEAST (v1.7.1) (Drummond et al. 2012) to investigate the temporal, spatial, and demographic evolution of ST22. The core genome was retained for the analyses after discarding all ambiguous bases and strains that showed admixture, resulting in a final alignment of 193 taxa and 2,154,246 bp. The different clades of ST22 (A and non-A) have different demographic histories, as inferred from the shape of their coalescent trees. Non-ST22-A individuals coalesce further back in time and at a slower rate than ST22-A individuals, and the latter have a more star-like genealogy. For this reason, we analyzed these two clades separately (25 taxa for non-ST22-A and 162 taxa for ST22-A). We used an uncorrelated lognormal distribution (variable clock rate) to model the rate of evolution and constrained the tips of the phylogeny to their dates of isolation to calibrate that rate, as well as giving it a default uniform prior of between 0 and 1. We used the HKY85 model of nucleotide substitution (Hasegawa et al. 1985) with four discrete gamma-distributed rate categories, and to model the relative node ages, ran the two data sets under both the constant size and exponential size coalescent tree prior. All other default parameters in the program BEAUti v1.6.2 were used. To model the global geographical spread, we calculated the probabilities of transfer between each discrete location (i.e., country) via a continuous-time Markov chain with a nonreversible infinitesimal rate matrix, using this to find the location at ancestral nodes (Lemey et al. 2009; Edwards et al. 2011). To reduce the number of parameters to estimate, we utilized Bayesian stochastic search variable selection (BSSVS), where at any point in the Monte Carlo Markov chain (MCMC) the transition probabilities of some states can be set to zero, with the prior designed to minimize the number of location changes as described in Lemey et al. (2009). To infer the geographic spread within the UK and to predict the root state of the fluoroquinolone-resistant isolates, we modeled diffusion over a continuous two-dimensional space assuming a Brownian motion (or Weiner process), with the rate of diffusion assumed to be constant over the entire tree and estimated from the data (Lemey et al. 2010). For this analysis, we used ST22-A strains and discarded strains that were fluoroquinolone-sensitive and those that did not have exact latitude and longitude information, leaving a total of 60 taxa. The uncertainty in root location is represented as 80% of the highest posterior density of the continuous traits: latitude and longitude. We calculated a posterior distribution of all parameters using two separate MCMCs and terminated the runs after checking for efficient mixing, convergence, and checking that 90% of the distribution was post “burn-in.” To estimate the timing of the first introduction to each country (Supplemental Fig. S1), we took the date of the youngest node for any individual location in the posterior distribution of trees.

Data access

The data have been deposited in the European Nucleotide Archive (ENA) (http://www.ebi.ac.uk/ena/) under study numbers ERP000150 and ERP000633.