Demography and Intercontinental Spread of the USA300 Community-Acquired Methicillin-Resistant Staphylococcus aureus Lineage.

UNLABELLED: Community-acquired methicillin-resistant Staphylococcus aureus (CA-MRSA) was recognized worldwide during the 1990s; in less than a decade, several genetically distinct CA-MRSA lineages carrying Panton-Valentine leukocidin genes have emerged on every continent. Most notably, in the United States, the sequence type 18-IV (ST8-IV) clone known as USA300 has become highly prevalent, outcompeting methicillin-susceptible S. aureus (MSSA) and other MRSA strains in both community and hospital settings. CA-MRSA bacteria are much less prevalent in Europe, where the European ST80-IV European CA-MRSA clone, USA300 CA-MRSA strains, and other lineages, such as ST22-IV, coexist. The question that arises is whether the USA300 CA-MRSA present in Europe (i) was imported once or on very few occasions, followed by a broad geographic spread, anticipating an increased prevalence in the future, or (ii) derived from multiple importations with limited spreading success. In the present study, we applied whole-genome sequencing to a collection of French USA300 CA-MRSA strains responsible for sporadic cases and micro-outbreaks over the past decade and United States ST8 MSSA and MRSA isolates. Genome-wide phylogenetic analysis demonstrated that the population structure of the French isolates is the product of multiple introductions dating back to the onset of the USA300 CA-MRSA clone in North America. Coalescent-based demography of the USA300 lineage shows that a strong expansion occurred during the 1990s concomitant with the acquisition of the arginine catabolic mobile element and antibiotic resistance, followed by a sharp decline initiated around 2008, reminiscent of the rise-and-fall pattern previously observed in the ST80 lineage. A future expansion of the USA300 lineage in Europe is therefore very unlikely. IMPORTANCE: To trace the origin, evolution, and dissemination pattern of the USA300 CA-MRSA clone in France, we sequenced a collection of strains of this lineage from cases reported in France in the last decade and compared them with 431 ST8 strains from the United States. We determined that the French CA-MRSA USA300 sporadic and micro-outbreak isolates resulted from multiple independent introductions of the USA300 North American lineage. At a global level, in the transition from an MSSA lineage to a successful CA-MRSA clone, it first became resistant to multiple antibiotics and acquired the arginine catabolic mobile element and subsequently acquired resistance to fluoroquinolones, and these two steps were associated with a dramatic demographic expansion. This expansion was followed by the current stabilization and expected decline of this lineage. These findings highlight the significance of horizontal gene acquisitions and point mutations in the success of such disseminated clones and illustrate their cyclic and sporadic life cycle.

INTRODUCTION

Staphylococcus aureus remains one of the most challenging and costly sources of bacterial infection worldwide. It asymptomatically colonizes about one-third of the human population and may cause infections with outcomes that range from mild to severe and are occasionally life threatening (1). Notably, it is the most common causative agent of nosocomial infections and a leading cause of death in hospitalized patients (2). Until the mid-1990s, methicillin-resistant S. aureus (MRSA) infections were reported exclusively in hospital settings and most hospital-associated MRSA (HA-MRSA) diseases resulted from a limited number of successful clones (3). However, in the beginning of the 2000s, MRSA infections began to be reported in healthy individuals without risk factors or discernible connections to health care institutions (4, 5). These community-acquired MRSA (CA-MRSA) strains had genetic backgrounds distinct from those of traditional HA-MRSA strains. Moreover, international strains of CA-MRSA belonged to a series of different lineages and specific clones were predominant on different continents (6). The dominant CA-MRSA clone in the United States, referred to as the USA300 North American (USA300-NA) MRSA clone in this study, belongs to pandemic multilocus sequence type 8 (ST8) constituted by methicillin-susceptible S. aureus (MSSA) and MRSA strains with the same pulsed-field gel electrophoresis type, USA300. This epidemiologically successful clone carries the IVa subtype of the staphylococcal cassette chromosome mec (SCCmec) element, agr allele 1, arginine catabolic mobile element (ACME) type I, and a set of virulence genes including lukS-PV/lukF-PV, sek, and seq (7–10).

It is generally recognized that this USA300-NA MRSA clone descended from an ancestral USA500-like MSSA strain by the acquisition of various mobile genetic elements (MGEs) and shows very recent clonal expansion (2, 11). These MGEs, namely, SCCmec type IV carrying the mecA gene, S. aureus pathogenicity island 5 containing sek and seq, phage phiSA2 carrying the Panton-Valentine leukocidin genes, and ACME type I, are thought to contribute to the success and high virulence of the commonly known USA300-NA MRSA strains (12–16). ACME type I is found exclusively in the USA300-NA MRSA lineage (16–21). It has been shown that the functional modularity (arginine deiminase system [arc] and the speG gene, which encodes a polyamine resistance enzyme) plays a major role in the enhanced success of this clone during colonization and skin infections (22, 23). In recent years, multiple studies have shown that the USA300-NA MRSA clone has spread worldwide, with reports of outbreaks in South America, the Middle East, the western Pacific, and Europe (for a detailed review, see reference 24). Interestingly, in South America, the most prevalent CA-MRSA is a USA300 variant, the so-called Latin American variant or USA300-LV (i.e., ST8, spa type t008, SCCmec type IVc, lukSF-PV+, arcA) present in Colombia, Venezuela, and Ecuador (24–26). It has recently been proposed that the USA300 MRSA lineage has evolved since the 1980s into two parallel epidemics, one in North America with the acquisition of the ACME by the ancestral USA300 lineage, forming the USA300-NA MRSA epidemic, and one in South America, with the acquisition of a novel copper and mercury resistance (COMER) element by the ancestral lineage (27).

In Europe, CA-MRSA epidemiology is characterized by clonal heterogeneity, although the most common European strain is the European clone (ST80-IV, pvl+) (28). This clone originally emerged from a West African ancestor (29) and expanded in the Mediterranean area, the Middle East, and North Africa, as many of the first patients infected with this clone in Europe had histories of travel to these regions. Although the prevalence of this “European” clone is relatively high in some countries, like Algeria and Tunisia, with a widespread distribution in both hospital and community settings (28–30), in Europe, the prevalence of this clone, as well as that of CA-MRSA in general, remains low (31, 32) and even seems to be declining (33, 34). This decay is in keeping with the calculation of a recent slow decrease in its effective population size based on a Bayesian skyline model (29). However, the question arose as to whether the CA-MRSA clone (USA300-NA or USA300-LV) could replace ST80 CA-MRSA and expand in a worrisome scenario similar to the one observed in the United States in the past decade. Consistent with this hypothesis is the fact that the USA300-NA MRSA clone has been reported to be present in many European countries for several years (31, 32, 35–45) and is occasionally associated with infection clusters (46–48). A recent study showed that an increase in CA-MRSA USA300 prevalence in 2013 in Geneva, Switzerland, was due to distinct importations from the North and South American continents (49).

This study addressed the question of the USA300 lineage’s capacity to expand successfully in Europe by applying whole-genome sequencing to a collection of French USA300 CA-MRSA strains responsible for sporadic cases, as well as micro-outbreaks, over the past decade, followed by a comparison with available genomes of USA300 CA-MRSA, USA300 MSSA strains, and other ST8 isolates from the United States (50, 51). Genome-wide phylogenetic relationships and coalescent-based analyses showed that the population structure of the French isolates reflects multiple introductions of the USA300-NA MRSA clone. Furthermore, the coalescent-based demography of USA300-NA lineage confirmed its success in the late 1990s but also suggested a sharp decline initiated around 2008, reminiscent of the rise-and-fall pattern observed in the ST80 lineage.

RESULTS

USA300 phylogeography and resistance makeup.

Prior to phylogenetic reconstruction, we checked if the data set provided some evidence of homologous recombination. We first made a visual inspection of the concatenated single-nucleotide polymorphisms (SNPs), and no contiguous SNPs consisting of three or more SNPs, mirroring homologous recombination, were detected. In a second step, neighbor nets were inferred in order to detect putative recombination signatures that would result in networks rather than trees. No major splits were found, and the pairwise homoplasy index (PHI) test failed to detect recombination signatures (P = 0.475). We then performed whole-genome phylogenetic analysis of the 498 ST8 isolates (including MRSA, MSSA, and USA300-LV and USA300-NA MRSA isolates) obtained from France (n = 67) and the United States (n = 431) within a time span of 15 years (see Table S1 in the supplemental material). A total of 12,840 SNPs were used to generate a maximum-likelihood (ML) tree that shows, as previously reported, that strains from San Diego (51) and New York City (50) are interspersed (Fig. 1). The closest relative sister group of the USA300-NA lineage corresponded to seven strains, all isolated from the New York City area. These seven isolates carry the COMER element characteristic of the USA300-LV South American clade (27) and belong to this lineage. Adding French USA300 CA-MRSA strains to the phylogenetic analysis did not unravel any new branching pattern or hidden sublineages and revealed that all of the French isolates belong to the USA300-NA lineage. Among the U.S. isolates, a fluoroquinolone-resistant lineage with the same mutations in gyrA (Ser84Leu) and parC (Ser80Tyr) has emerged within the USA300-NA MRSA clone and disseminated globally (Fig. 1). We observed a similar distribution in France, with 18 fluoroquinolone-susceptible strains, including 10 from an outbreak in Le Puy-en-Velay (a town in central France) (47) and 46 resistant isolates, including the 28 isolates from an outbreak in a long-term care facility (Paris 16th District) showing the same two mutations in the gyrA and parC loci. Plasmid analysis showed that all except three French isolates carry plasmids related to previously described plasmid p18805-p03 (52). This plasmid encodes multiple antibiotic and heavy metal resistance genes [aminoglycosides, ant(6)-Ia and aph(3′)-III; beta-lactams, blaZ; macrolides mph(C); macrolides-lincosamides-streptogramin B, msr(A); cadmium, cadD and cadX]. These plasmids showed almost 100% nucleotide sequence identity indicative of a single origin but different regions of deletion. Similar plasmids are also present in the vast majority of the American isolates (50, 51) but absent from the strains of the USA300-LV lineage. It is therefore likely that an ancestor of this multiresistance plasmid was acquired by the common ancestor of the USA300-NA clone and contributed to its expansion. The 28 isolates from the outbreak in a long-term care facility (Paris 16th District) carried in their chromosome a 3.2-kb transposable element containing the dfrG gene, which encodes a dihydrofolate reductase conferring trimethoprim resistance.

FIG 1 

Phylogenetic reconstruction of ST8 and its derived USA300 lineage. ML tree based on 498 genomes and a total of 12,840 concatenated SNPs, rooted by using distantly related MSSA strain ERS092996. Branch colors indicate the geographic origins of the strains: red for New York State; yellow for California; green for Minnesota, and blue for France. Outer circles represent ACME and antibiotic resistance profiles, whereas the inner circle represents the years of isolation of the strains. The stars correspond to 12 predicted independent intercontinental USA300-NA CA-MRSA transfers from North America to France. The different assignments in the ST8 “sublineages” are also indicated by colored areas within the circular tree. Fq, fluoroquinolone.

Temporal Bayesian analysis reveals the complexity of USA300 demography.

In a second step, we used the Bayesian coalescent analysis implemented in BEAST (53) to generate a phylogenetic time tree and to derive estimated dates for the common ancestors of the strains (see Fig. S1 in the supplemental material). The tree was calibrated by using the sampling dates of the isolates ranging from 2000 to 2014. We tested the performance of various demographic models that favored a Bayesian skyline model with a relaxed molecular clock (data not shown). The estimated short-term mutation rate corresponded to 1.34 × 10−6 (95% confidence interval of 1.11 × 10−6 to 1.56 × 10−6) substitutions per nucleotide site per year, similar to that reported previously in USA300 MRSA (50, 54) and other S. aureus evolutionary lineages (55, 56). On the basis of this substitution rate, the estimated time of the most recent common ancestor (TMRCA) of the USA300-NA clone is 1996 (see Fig. S1 in the supplemental material), a result that slightly shifts the ancestor of the CA-MRSA clone toward the 21st century, compared to the estimates obtained by Uhlemann and colleagues (between 1970 and 1993) (50). In addition, our analyses revealed that the acquisition of ACME, methicillin resistance, and a multiresistance plasmid dates back to 1996 and that fluoroquinolone resistance was gained around 1998 (see Fig. S1 in the supplemental material). It is noteworthy that, because of the lack of intermediate strains harboring only one of those mobile elements, we could not determine whether the acquisition of the SCCmec element harboring the methicillin resistance gene mecA, the multiresistance plasmid, and the ACME type I element occurred simultaneously or one event occurred very shortly after the other. The divergence with the USA300-LV MRSA South American lineage dates back to 1983, an estimate slightly more recent than that proposed by Planet et al. (around 1975) (27).

We then generated a Bayesian skyline plot that estimated the pathogen’s demographic changes over time (Fig. 2) on the basis of the 498 ST8 genomes. The ST8 USA300 effective population size was relatively stable from the early 1960s to the mid-1990s. Then, a bottleneck followed, giving rise to the so-called USA300-NA clone. Interestingly, the USA300-NA clone went through a sharp two-phase expansion event that perfectly matches the acquisitions of the ACME and the antibiotic resistance elements, thus giving rise to the USA300-NA MRSA lineage, followed by the acquisition of fluoroquinolone resistance, as mentioned above. Strikingly, each of these major gains of function was accompanied by a 1-order-of-magnitude population size increase, resulting in a global 100-fold effective population size increase (Fig. 2). However, after 2008, the coalescent-based analysis detects a sharp decline of the USA300-NA MRSA clone. These results need further investigation to see if ongoing epidemiological surveys can confirm this model. We next searched for additional SNPs under positive selection that might explain the success of the USA300 lineage. To address this issue, we used the PCAdapt software (57, 58) based on a Bayesian factor model to identify candidate mutations (Fig. 3). In Table S2 in the supplemental material are shown the 20 SNPs with the highest scores (log10 Bayesian factors of >6 × 1050). Although such statistic-based predictions must be viewed cautiously and require experimental validation, the 20 mutations were all located in protein coding sequences (see Table S2 in the supplemental material). The six synonymous mutations might affect the expression of the encoded protein. Strikingly, several mutated genes might be involved in the interaction with the human host or the environment, as they encode surface proteins, a lipoprotein, an autolysin, a putative transporter involved in teichoic acid transport, and a protein involved in lipid biosynthesis.

FIG 2 

Bayesian skyline plot indicating population size changes in the USA300 lineage over time with a relaxed molecular clock. The shaded area represents the 95% confidence interval, and the arrows point to major phenotypic events driven by ACME and antibiotic resistance that might have contributed directly to the success of the USA300-NA MRSA lineage.

FIG 3 

SNP-based Bayesian factor model analysis for detecting genes involved in positive selection in the USA300-NA lineage. (A) Latent factors of the 12,840 SNPs and 498 strains with the first two factors. (B) Manhattan plot representing the selection scan and the outliers that are related to the different latent factors. The dotted line highlights the top 1% of the SNPs associated with the highest Bayes factor (BF) values.

French isolates are interspersed in United States isolate collections.

Strikingly, the 67 strains in the French panel clustered in 12 lineages of 1 to 38 isolates. Furthermore, these clusters are scattered within the North American MRSA isolates (USA300-NA) (Fig. 1). This indicates a common origin of USA300-NA isolates and European USA300 CA-MRSA isolates, likely resulting from several transfers from the United States to Europe. Interestingly, only one case could be related to travel to the United States, in 2003 (strain HT20030124 from Annecy; see Table S1 in the supplemental material).

We then calculated the genome-wide pairwise distances of the most closely related French-American strains and compared them with pairwise distances between strains collected from the same patient at multiple sites, multiple colonies from the same nasal sample, or strains collected from the same patient upon different hospital admissions (59, 60) (Fig. 4). The median number of SNP differences that separated the closest French-American epidemic link was 74 (range, 59 to 120) and significantly differed from the other two settings mentioned above (P > 0.0001). However, this number of differences remains smaller than the median number of SNP differences (n = 104) for USA300 CA-MRSA isolates collected from different households in a New York City community (50). In addition, bootstrap support for clades encompassing the French strains and their closest American homologs was >70% for >50% of the different subsamples (see Fig. S2 in the supplemental material). Therefore, the epidemic links we detected between American and French strains are likely to reflect direct, or at least very short, transmission paths.

FIG 4 

Box plots representing pairwise SNP comparisons of strains from the same patient (black), the closest strains from France and their American relative (blue), and strains from different hospital admissions of the same patient (grey). These data were gathered from the present study and those of Golubchik et al. (59) and Price et al. (60).

The largest French cluster contains 28 isolates from an outbreak in a Parisian long-term care facility but also 3 strains isolated in Aubervilliers nearby Paris, as well as 6 unlinked isolates from different parts of France (Fig. 5A and B). This cluster is very likely the product of a single introduction from the United States. Its TMRCA is 1998 (Fig. 5C; see Fig. S1 in the supplemental material), suggesting that this subclone had been circulating for at least a decade in France before being noticed.

FIG 5 

Evolutionary relationships of the strains from an outbreak at a long-term care facility (LTCF) in the Paris area. (A) Minimum spanning tree based on SNP data. Circle sizes are a function of the number of strains with the same genotype. Dashed lines correspond to links with >25 SNPs; thin lines correspond to links with <25 SNPs, and thick lines correspond to links with ≤10 SNPs. Dates of isolation are indicated. Colors represent the geographical or epidemiological origins according to the key. (B) ML tree based on the alignment of 488 SNPs of the dominant French lineage. Numbers of SNPs are mapped on the branches. Strains from the outbreak group in three clusters colored green, blue, and red. (C) Posterior distribution of the Paris LTCF TMRCA generated with the BEAST algorithm. The common ancestor (median) appeared in Europe in 1998 (95% highest posterior density, 1996 to 2003).

DISCUSSION

Here we show that, in most instances, the USA300 CA-MRSA cases identified in France in the last decade corresponded to sporadic and independent importations from the United States (USA300-NA MRSA clone) without further indication of spreading in French territory. This is consistent with another observation in Switzerland, where the local increased prevalence of USA300 CA-MRSA in 2013 in Geneva was shown to result from multiple importations from America (49). As the Institut de Veille Sanitaire (a disease surveillance agency of the Ministry of Health) strongly recommend referring strains associated with multiple cases of skin and soft tissue infections (SSTIs) (within a household or other communities) to the French National Reference Laboratory, it is thus likely that the lack of observed spreading truly is due to the absence or a very limited number of secondary cases and not to underreporting. The reasons for the apparent lack of success of these USA300-NA CA-MRSA isolates in diffusing within the French community are unknown. However, it is in keeping with the low prevalence of CA-MRSA observed in Europe and specifically in France (32), despite the isolation of various CA-MRSA lineages on the European continent as early as 1993 (38, 61). Moreover, these observations also argue against the hypothesis that the USA300-NA MRSA clone is intrinsically more successful than European CA-MRSA ST80 (3, 62, 63) and thus that the USA300-NA MRSA clone could be a potential threat in Europe.

In one instance of the present study, however, one imported lineage appeared to have a remarkable geographic spread, which was subsequently responsible for two outbreaks, one in a long-term care facility (Paris 16th District) and another limited outbreak in a facility for the disabled also in the Paris area (Aubervilliers) in 2013 and 2014 (Fig. 5A). The same lineage was also responsible for sporadic cases scattered in French territory, i.e., in Paris (75019, June 2011), northeastern France (Strasbourg, January 2013), the French Alps (Grenoble, September 2011), central France (Orleans, January 2011), and western France (St-Nazaire, June 2012). Our Bayesian analysis predicts a TMRCA for this clone of 1998 (Fig. 5C). This observation suggests that introduction of the USA300-NA MRSA clone can be successful but to a limited extent. We did not detect genes specific to this lineage or determine which genetic mechanisms might be behind the successful dissemination of certain isolates.

A recent report by Planet et al. (27) suggested that the history of the USA300 CA-MRSA lineage followed a parallel evolution since the mid-1970s on the American continents; one epidemic sublineage disseminated in South America (so-called USA300-LV) and another disseminated in North America (referred to here as USA300-NA MRSA). The separation of these two sublineages is concomitant with the acquisition of the ACME by the North American USA300 MRSA clone, while USA300-LV acquired another mobile element—the COMER element. They also acquired, likely independently, two different subtypes of the SCCmec elements, type IVc in the LV clade and type IVa in the NA clade. The demography of the USA300-NA lineage based on Bayesian analysis is in agreement with this model, showing a considerable demographic expansion around 1996, which coincides with the acquisition of the ACME and the SCCmec element by this USA300-NA lineage, followed 2 years later by a second expansion peak that likely corresponds to the acquisition of fluoroquinolone resistance. This period of expansion assessed by Bayesian analysis is in line with the observed global spread of the USA300-NA CA-MRSA clone in North America and worldwide, although with a more limited impact in Europe (9, 32, 35, 37, 38, 44, 45, 64–67). Strikingly, the effective population size of the lineage is predicted to be declining during the first decade of this century, a prediction that should be interpreted cautiously since there is no surveillance of the prevalence of USA300-NA CA-MRSA at the worldwide level. However, a recent meta-analysis retrieving 1,004 publications covering the three continents showed that USA300 CA-MRSA is far from expanding on the various continents (68) and is instead declining in many countries (33, 34). Observations on the European continent suggest that the USA300-NA lineage, which has been described for decades as having a rather low prevalence and, according to our study, was imported into France on multiple instances many years ago, is apparently not expanding (32, 34, 44, 69, 70). Overall, it appears that the USA300-NA clone is behaving like many other MRSA clones, in particular, HA-MRSA, for which a model of cyclic periods of expansion, equilibrium, and decline has been observed (56, 71–76). It is worth mentioning that in only very few cases could a plausible account be given to explain the final population expansions, such as acquisition of resistance to antibiotics or accumulation of mutations and increased fitness (71–73). In the cases of USA300 CA-MRSA and the worldwide emergence of CA-MRSA in general, there is no satisfactory explanation for the concomitant emergence and expansion of the various CA-MRSA clones in different countries or on different continents at the end of the 20th century (5).

The most plausible scenario explaining CA-MRSA clone dynamics seems to be a complex combination of acquisition and loss of traits under selection (resistance to antibiotics, resistance to human host defenses), stochastic random processes, and eventually, replacement of bacterial communities by others competing for the same niche. The latter factor might explain the population decline observed in our Bayesian analyses; however, we have to keep in mind that the natural habitat of S. aureus is the nose and the skin. Therefore, focusing on population genomics of only pathogenic lineages might blur our global understanding of species communities and population dynamics.

MATERIALS AND METHODS

Bacterial strains.

The bacterial strains used in this study were from 24 sporadic cases that were geographically and temporally dispersed and from three limited outbreaks reported to the French National Reference Center for Staphylococci (CNR Staphylococci). All strains were isolated from cases (infection and/or carriage) that occurred in continental French territory in the past 11 years (see Table S1 in the supplemental material). The first outbreak occurred in Le Puy en Velay, a small town located in a rural area in Auvergne, France, and involved 12 individuals with community-acquired SSTIs, i.e., 6 children attending the same day care facility and 6 adult relatives in six households, between June 2011 and May 2012 (47). Two outbreaks occurred in long-term health care units. One was in a nursing home for mentally disabled patients in Aubervilliers in the Paris area; four patients were detected as carriers or infected with USA300 CA-MRSA between July 2013 and February 2014. The other outbreak occurred in a long-term care facility in the Paris 16th District between June 2012 and December 2013 and involved 10 infected patients and 12 carriers.

Four strains from the United States corresponded to sporadic cases of SSTIs that occurred before April 2003 and were provided by Barry Kreiswirth (77); three other strains originated from patients in Minnesota before 2001 and were provided by Timothy Naimi (78). The genomic backgrounds of these strains were assessed by diagnostic DNA microarray analysis (Identibac S. aureus Genotyping; Alere, Jena, Germany) as previously described (79). This microarray covers 332 different target sequences corresponding to 185 distinct genes and their allelic variants. The affiliation of isolates with clonal complexes was determined by a comparison of their hybridization profiles with those of reference strains previously characterized by multilocus sequence typing (79).

Published genomic data used for comparative analysis.

Genomic data from 36 isolates from the first large-scale genomic studies on USA300 CA-MRSA were retrieved (51). These strains had been isolated between 2003 and 2007 in California. In addition, the genomic data from 387 ST8 isolates from an extensive analysis of a New York City community performed between 2009 and 2011 were included. Metadata from these two studies were graciously provided by the authors (50, 51).

DNA sequencing and SNP detection.

S. aureus genomes were sequenced by using the Illumina HiSeq 2000 (101 nucleotide reads) or MiSeq (150 nucleotide reads) sequencer, with a coverage of >75×. Libraries were constructed with the Illumina TruSeq kit. Sequence reads from previously described USA300 strains were downloaded from the European Nucleotide Archive website (http://www.ebi.ac.uk/ena) or obtained directly from the authors. Sequence reads were aligned with the first completely sequenced and annotated US300-NA genome, FPR3757 (GenBank accession no. CP000255.1) (10) by using the Burrows-Wheeler Alignment tool (BWA mem 0.7.5a) (80). SNP calling was done with the Genome Analysis Toolkit (GATK 2.7-2) (81) Unified Genotyper by following Broad Institute best practices. Candidate SNPs were further filtered by requiring coverage of greater than half of the genome mean coverage and 95% read agreement to validate the call. SNPs selected with this stringent filter were regarded as true positives and were searched for in the original set for every strain in which they had been filtered out. SNPs, short indels, and coverage were visualized with SynTView (82). Specific analysis of clones from the main French lineage was done by manually checking read alignments around SNPs with Tablet (83). Four regions corresponding to MGEs not present or showing a high density of SNPs in other USA300 isolates and in closely related ST8 strains were removed from the analysis, leaving a total of 12,840 polymorphic sites. These regions are the ACME and the SCCmec cassette (41 to 125 kb), pathogenicity island 3 (encoding enterotoxins K and Q; 881 to 896 kb) (84), and integrated phages phiSA2usa (encoding the Panton-Valentine toxin; 1,546 to 1,644 kb) and phiSA3usa (2,084 to 2,127 kb). These regions represent a total of 240 kb or 8.3% of the strain USA300_FPR3757 chromosome.

Accessory genome analysis.

Plasmids and chromosomal regions specific to the French isolates compared to the FPR3757 reference genome were analyzed following de novo assembly of unmapped reads. Assembly was performed with Velvet by using an optimized k value (85). Plasmids and phages were identified by BLASTn searches with plasmid and phage sequences from S. aureus as the queries. Antibiotic resistance genes were identified with the Center for Genomic Epidemiology web tool (http://www.genomicepidemiology.org/).

Recombination detection.

Most of the analyses developed in our analytical framework (phylogenetics and Bayesian inference) are based on the assumptions that S. aureus evolution is mostly clonal and that recombination can be neglected. Therefore, in a preliminary step, we tested for the presence of mosaic genomes with the algorithm SplitsTree v4.13.1 (86). Putative recombination signatures were inferred with Neighbor-Net (87), and each data set was analyzed for the presence of recombinant sequences with the PHI test in SplitsTree with an alpha value of 0.001.

Phylogenetic analyses.

Phylogenic reconstructions were performed by considering the 12,840 polymorphic sites retained in the core genome of the 498 isolates. Phylogenetic relationships were reconstructed by the ML approach implemented in PhyML 3.0 (88). The robustness of the ML tree topology was assessed with bootstrapping analyses of 1,000 pseudoreplicated data sets. A transversion substitution model was selected on the basis of Akaike’s information criterion with jModelTest 2.1.3 (89). Phylogenies were rooted with MSSA strain ERS092996.

Coalescent-based analyses.

Evolutionary rates and tree topologies were analyzed with the generalized time reversible (GTR) and Hasegawa-Kishino-Yano (90) (HKY) substitution models with gamma distributed among-site rate variation with four rate categories (Γ4). We tested both a strict molecular clock (which assumes the same evolutionary rate for all of the branches of the tree) and a relaxed clock that allows different rates among the branches. Constant-size, logistic, exponentially growing coalescent models were used. We also considered the Bayesian skyline plot model (91), based on a general, nonparametric prior that enforces no particular demographic history. We used a piecewise linear skyline model with 20 groups and then compared the marginal likelihood of each model with Bayes factors estimated in Tracer 1.5. Bayes factors represent the ratio of the marginal likelihood of the models being compared. Approximate marginal likelihoods for each coalescent model were calculated via importance sampling (1,000 bootstrap replications) with the harmonic mean of the sampled likelihoods. A ratio between 3 and 10 indicates moderate support of the idea that one model fits the data better than another, whereas values of >10 indicate strong support. For each analysis, two independent runs of 100 million steps were performed and the chain was sampled every 10,000th generation. Examination of the Markov chain Monte Carlo (MCMC) samples with Tracer 1.5 indicated convergence and adequate mixing of the Markov chains, with effective sample sizes for each parameter in the hundreds or thousands. The first 10% of each chain was discarded as burn-in. We found the maximum clade credibility topology with TreeAnnotator 1.7.5 (52), and we reconstructed the Bayesian skyline plot with Tracer 1.5. The relaxed clock models provided a better fit to the data (Bayes factor, >12) and under the different models tested, the Bayesian skyline model provided a (marginally) better fit, overall.

Analysis of genes under positive selection.

To capture SNPs under positive selection, we used the PCAdapt software to perform a genome scan based on a Bayesian factor model (57). We chose K = 2 factors because the third and the fourth factors did not correspond to population structure and distinguished individuals within the same clades. The factor analysis was performed on the centered genotype matrix that was not scaled. The MCMC algorithm was initialized by singular value decomposition, and the total number of steps was 400 with a burn-in of 200 steps.

Time-calibrated Bayesian evolutionary tree of the ST8 and USA300 lineages produced in BEAST under the HKY substitution and Bayesian skyline models. The vertical bars correspond to phenotypic and resistance trait acquisitions. The strains highlighted in red correspond to the fluoroquinolone-resistant strains of the USA300-NA sublineage. Download

Figure S1, PDF file, 0.9 MB

Bootstrap consensus tree based on 498 S. aureus strains by using 12,840 variable positions. French isolates are colored blue, and node support after 100 bootstrap replications is shown on the branches. We used the subtree-pruning-regrafting maximum-parsimony search method implemented in MEGA 5.10 (K. Tamura, D. Peterson, N. Peterson, G. Stecher, M. Nei, and S. Kumar, Mol Biol Evol, 28:2731–2739, 2011) to generate phylogenetic trees. In order to test the phylogeny’s robustness, a bootstrap method with 100 replications was used. We then selected the bootstrap consensus tree. As a cladogram, branch lengths are not correlated with the number of SNP differences. Finally, the software FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/) was used to generate the graphic representation. Download

Figure S2, PDF file, 0.1 MB

Strains used in this study, their geographical origins, and the presence of the ACME and of nonsynonymous mutations associated with fluoroquinolone resistance.

Table S1, XLS file, 1.7 MB

List of SNPs acquired by the USA300-North American last common ancestor and under putative positive selection.

Table S2, DOCX file, 0.1 MB