Historic methicillin-resistant Staphylococcus aureus: expanding current knowledge using molecular epidemiological characterization of a Swiss legacy collection

Abstract

Background

Few methicillin-resistant Staphylococcus aureus (MRSA) from the early years of its global emergence have been sequenced. Knowledge about evolutionary factors promoting the success of specific MRSA multi-locus sequence types (MLSTs) remains scarce. We aimed to characterize a legacy MRSA collection isolated from 1965 to 1987 and compare it against publicly available international and local genomes.

Methods

We accessed 451 historic (1965–1987) MRSA isolates stored in the Culture Collection of Switzerland, mostly collected from the Zurich region. We determined phenotypic antimicrobial resistance (AMR) and performed whole genome sequencing (WGS) using Illumina short-read sequencing on all isolates and long-read sequencing on a selection with Oxford Nanopore Technology. For context, we included 103 publicly available international assemblies from 1960 to 1992 and sequenced 1207 modern Swiss MRSA isolates from 2007 to 2022. We analyzed the core genome (cg)MLST and predicted SCCmec cassette types, AMR, and virulence genes.

Results

Among the 451 historic Swiss MRSA isolates, we found 17 sequence types (STs) of which 11 have been previously described. Two STs were novel combinations of known loci and six isolates carried previously unsubmitted MLST alleles, representing five new STs (ST7843, ST7844, ST7837, ST7839, and ST7842). Most isolates (83% 376/451) represented ST247-MRSA-I isolated in the 1960s, followed by ST7844 (6% 25/451), a novel single locus variant (SLV) of ST239. Analysis by cgMLST indicated that isolates belonging to ST7844-MRSA-III cluster within the diversity of ST239-MRSA-III. Early MRSA were predominantly from clonal complex (CC)8. From 1980 to the end of the twentieth century, we observed that CC22 and CC5 as well as CC8 were present, both locally and internationally.

Conclusions

The combined analysis of 1761 historic and contemporary MRSA isolates across more than 50 years uncovered novel STs and allowed us a glimpse into the lineage flux between Swiss-German and international MRSA across time.

Background

The introduction of penicillin during the 1930s and early 1940s was quickly followed by the rise of penicillin-resistant Staphylococcus aureus due to penicillinase blaZ [1,2,3]. Methicillin, a penicillinase-resistant antibiotic, was introduced in 1961, and within a year, the first methicillin resistant S. aureus (MRSA) isolate carrying the mecA gene emerged and established itself as a major nosocomial pathogen [4]. The mecA gene is found on a mobile genetic element, the Staphylococcal cassette chromosome mec element (SCCmec), of which diverse types have evolved, and which spreads via horizontal gene transfer [5]. The acquisition of mecA by S. aureus has been dated to the mid-1940s [6], and it is hypothesized that the SCCmec cassette originated from the mecA gene of Staphylococcus fleurettii and its surrounding chromosomal region [7]. While one study advanced the hypothesis that the acquisition of the SCCmec cassette was a single evolutionary event [8], more recent studies supported the theory that MRSA emerged multiple times through independent events [9,10,11] and that cassette substitution can also occur, although far less frequently than cassette acquisition [12]. Molecular typing of MRSA isolates is an especially important part of any local, national, and international epidemiology management strategy [13]. Historic MRSA studies in European countries with low MRSA prevalence used various sampling strategies and typing methods [14, 15]. Since MRSA was largely a nosocomial issue for most of the latter twentieth century, most publications include samples isolated in hospitals either with no defined inclusion criteria or focused on specific types of infections such as bacteremia [16,17,18,19]. These publications employed genotyping methods such as phage typing, spa typing, and pulse field gel electrophoresis (PFGE). These techniques have been superseded by WGS and MLST, and direct comparisons between the older techniques and the molecular ones are often unclear. In more recent times, the falling costs of sequencing has enabled the high-resolution investigation of MRSA transmission in a variety of settings [20]. Due to resolution differences between technologies, re-sequencing of legacy collections is especially important to complement the historic picture. Current unified MRSA clone nomenclature incorporates information on the genomic ancestry as multi-locus sequence typing (MLST), methicillin resistance status, and the SCCmec cassette type (e.g., ST250-MRSA-I or ST8-MSSA) [9]. The emergence of antimicrobial resistant S. aureus throughout the twentieth century can be classified into four waves. First, the advent of penicillin resistance provided the first resistant S. aureus (1940-1960). Second came methicillin resistance, mostly linked to European hospitals (1960–1980), which spawned the term hospital-associated MRSA (HA-MRSA). In the third wave, methicillin resistance spread worldwide (1980–1990). The fourth wave defines the rise of community-associated and livestock-associated MRSA (CA-MRSA, LA-MRSA) (1990–2000) [21], coinciding with the development of resistance to ciprofloxacin and vancomycin [22, 23]. During this time, major clonal populations with specific STs emerged and disappeared through displacement by new, more successful epidemic clones [24]. The first epidemic MRSA clone was identified as ST250 in the 1960s, within clonal complex (CC)8. It is hypothesized that it arose as a single locus variant (SLV) of ST8-MSSA, which then acquired the SCCmec type I [9] and spread throughout Europe. ST250-MRSA-I was subsequently replaced by its SLV ST247-MRSA-I, first in Denmark around 1964 [17] and then globally, gaining further antimicrobial resistance [9]. ST247-MRSA-I was later displaced in Europe by other successful lineages in the late 1990s. ST239-MRSA-III, which also evolved from ST8-MRSA-III [9], was first identified in Australian hospitals in 1979 [25]. Subsequently, it was found in Brazilian hospitals in 1992 [26] and competed with ST247-MRSA-I in Portuguese hospitals throughout the 1990s [27, 28]. Later in the century, ST8-MSSA also acquired resistance to methicillin: one of the most successful contemporary CA-MRSA, USA300, belongs to ST8 [29, 30]. US300 has become endemic in North America and still finds niches to colonize around the globe, with repeated introductions of USA300 into Europe reported [29, 31,32,33,34]. ST22-MRSA-IV, another global lineage, has been reported in many countries since the 1990s. ST22-MRSA-IV became a dominant HA-MRSA in England (UK-EMRSA-15) [35,36,37] and by the year 2000 was responsible for 65% of British MRSA bacteremia episodes [38]. ST5 is a pandemic lineage reported in Europe, Asia, and North America [36]. Previous comparative WGS studies have highlighted a recent increase in MRSA diversity, combined with success in ecological niches [35, 39]. This success is bolstered by new genetic elements associated with either increased invasiveness or lower fitness cost [29, 40]. The epidemiological situation of MRSA within Switzerland was first assessed extensively at the national level in 1997 [41], where low rates of MRSA prevalence (< 5%) were reported, reflecting the situation in Germany, the Netherlands, Denmark, or Sweden [16]. A notable exception to this was Geneva University Hospital, which reported 23% MRSA prevalence, a figure comparable to hospitals of countries South of Switzerland [16]. This stark regional contrast was confirmed again in a study which explored the temporal trends of Swiss MRSA from 2004 to 2014 [42]. This study also showed how MRSA prevalence decreased from 26 to 12% in the French-speaking region and from 20 to 12% in the Italian-speaking region from 2004 to 2014, during which time the German region found a slight increase in prevalence (from 4 to 5%). The study lists some factors which could explain these regional differences, such as the geographical, cultural, and economic ties between the Swiss linguistic region and their neighboring countries or the fact that these countries adopted quite different control measures against MRSA dissemination. Following the 1997 report, Switzerland adopted a local approach to surveillance rather than a unified national approach, which may also have led to this regional heterogeneity [43]. Currently, surveillance of MRSA is conducted nationwide by the Swiss Center for antibiotic resistance ANRESIS [44], and its data show similar rates if MRSA in French-, Italian-, and German-speaking Switzerland [45]. Swiss MRSA epidemiology started being documented in the late 1960s by Prof. Fritz Kayser and his team at the Institute of Medical Microbiology at the University of Zurich [46,47,48]. The bacterial stocks generated during this research are now part of the Culture Collection of Switzerland (CCoS), a national repository for microorganisms [49]. These S. aureus clinical isolates were collected from Swiss hospitals in the greater Zurich area between 1965 and 1987 with inclusion based on scientific interest and availability from the diagnostic laboratory of the university hospital Zurich. At that time, no systematic sampling and biobanking was done. This is the time-period when MRSA was establishing itself as a major nosocomial pathogen, from when the number of publicly available whole genomes is quite limited, with just 103 from isolates between 1960 and 1992. Therefore, we aimed to perform an in-depth characterization of this Swiss legacy collection using whole genome sequencing and to position the results in a global context through comparison with publicly available genomes and modern Swiss MRSA.

Materials and methods

Whole genome sequencing and phenotypic AMR profiling

Bacterial genome assemblies from the CCoS and from the University Hospital Basel (USB) were generated at the Division of Clinical Bacteriology (USB) according to their ISO 17025 accredited standard procedures. DNA was extracted with a Qiagen BioRobot EZ1 using the QIAamp DNA Mini Kit (QIAGEN, Hilden, Germany) and according to the manufacturer’s guidelines. Library preparation was performed using an Illumina Nextera DNA Flex Library Prep Kit and multiplexing at 96-plex on a NextSeq 500 System using the Mid-Output Kit (Illumina, San Diego, USA). Resulting fastq files underwent a quality check, where the sequences with a phred score lower than 30 were discarded (only sequences with a base calling accuracy above 99.9% are kept). Sequencing data was quality controlled using FastQC (v 0.11.5) [50] and MetaPhlAn (v 2.7.7) [51]. Adaptors were trimmed using Trimmomatic (v 0.38) using default parameters (ILLUMINACLIP:2:30:10 SLIDINGWINDOW:4:15 MINLEN:125) [52]. Genome assemblies were created de novo using Unicycler (v 0.3.0b) [53] and checked with QUAST (v 5.0.2) [54] Further details can be found in [55]. For the isolates from CCoS, a phenotypic antibiotic resistance analysis was conducted as follows: bacteria were grown on CHROMID/MRSA Agar (bioMérieux, Marcy-l’Étoile, France) at 35 ± 2 °C for 48 h. Antibiotic susceptibilities were determined using a microdilution assay (Vitek2, bioMérieux, Marcy-l’Étoile, France) and interpreted according to breakpoints published by the European Committee on Antimicrobial Susceptibility Testing (EUCAST), version 9.0 [56].

Nanopore sequencing

Ten isolates from the CCoS were chosen for nanopore sequencing. DNA was extracted from the bacterial pellet using the DNA mini (Qiagen) kit on the QIAcube robot (Gram + Enzymatic lysis protocol) with 150 μl elution volume. Long read libraries were prepared using the ligation sequencing kit (SQK-LSK 109). The libraries were sequenced using GridION with R9.4 flow cells with a default 72 h run time. Basecalling and de-multiplexing were carried out within the inbuilt MinKnow (21.05.12) software. Basecalling was done using Guppy (5.0.12) high-accuracy model. Quality check was carried out using Nanoplot (1.35.4) [57] and MetaPhlAn (3.0.13) [51]. Trimming was done with Porechop (0.2.3) [58]. Short reads were removed with Filtlong (0.2.0) [59]. Assembly was performed with flye (2.8.1) [60] and polished with racon (1.4.7) [61], medaka (1.4.4) [62], and polypolish (v0.4.3) [63]. Raw data can be found at the European Nucleotide Archive under accession number PRJEB59014.

Collection of public genomes

The timeframe covered by the 451 S. aureus isolates from the CCoS is from 1965 to 1987. The parameters for the search of public genomes were set at a slightly wider timeframe (1960–1992) to allow for a wider selection of isolates. One hundred eight assembled genomes were downloaded from NCBI Pathogen Detection [64]. The 108 genomes, and their respective metadata, were obtained by applying the following criteria: species: Staphylococcus aureus, collection date: 1960 to 1992. One hundred fifty genomes were downloaded from Pathogenwatch [65, 66]. The assemblies were found by searching for Staphylococcus aureus genomes collected between January 1960 and December 1992.

Data filtering

Filtering criteria were applied to both the public and in-house sequenced isolates. Only isolates which satisfy the following conditions were kept: genome length within ± 10% of reference S. aureus genome NCTC 8325; where recorded, read depth greater than 20x; if antibiogram was present, resistance to oxacillin; otherwise, presence of mecA in the genomes was considered as denoting an MRSA. Forty-six isolates with an oxacillin-sensitive phenotype were removed to focus the analysis on MRSA; as oxacillin belongs to the same drug class, it shares its mode of action with methicillin and has replaced methicillin in clinical use. Twenty-three further isolates which possessed no phenotypic data and were predicted as MSSA during analyses were discarded. Additionally, 131 duplicate genomes were removed. Two CCoS isolates were removed for low read depth (< 20x). Finally, 1207 MRSA genomes isolated in multiple Swiss hospitals and sequenced at the University Hospital of Basel between 2007 and 2022 were selected to represent a modern collection. After this process, 451 MRSA isolates from the CCoS, 103 from public repositories (18 from Pathogenwatch and 85 from NCBI Pathogen Detection), and 1207 from the USB were used in further analysis. Metadata for all the isolates are available in Additional file 2.

Genotyping

This work used two genotypic typing methods: multi-locus sequence typing (MLST) and core genome MLST (cgMLST [67]). The genomes were typed in Ridom SeqSphere+ (v8.3.0), generating Minimum Spanning Trees (MSTs) of MLST or cgMLST data, and a world map displaying the origin of the isolates (Map data © Google, INEGI). For all visualizations, colorblind-friendly color palettes were created with the “coolors” website [68]. cgMLST clusters were generated with a minimum cluster size of 15 and a maximum cluster distance of 24 different alleles [67].

Phylogenies

To improve visualization comprehensibility for this step, the 1207 genomes from the USB were dereplicated with a 0.005 threshold in Assembly Dereplicator (v0.1.0) [69] resulting in 223 genomes. Assembly dereplicator clusters the genomes with dissimilarity lower than the threshold and keeps the assembly with the largest N50 for each cluster. An S. epidermis reference genome (NZ_CP035288.1) was added as an outgroup for rooting purposes. Prokka (v1.14.5) [70] was used to annotate the genomes. From the annotated genomes, the core genome alignment was generated with roary (v3.13.0) [70].This core genome alignment was processed with IQ-TREE 2 (v 2.2.1) [71]. Whole genome alignments of cgMLST clusters were calculated with SKA (v 1.0) [72], within-house hybrid assemblies used as a reference where available. For ST22, reference sequence NZ_CP053101.1 was used. From these alignments, phylogenies without recombination were generated with gubbins (v 3.2.1) [52] and rooted with BactDating (v1.1.1) [73]. Clusters were investigated with fastBAPS (v 1.0.8) [74]. The phylogenies were visualized in RStudio (v2022.07.1+554).

AMR and virulence prediction

In order to predict AMR and virulence, parts of the -finder software suite developed by the Danish center for genomic epidemiology were used. The Resfinder Software [75] was applied to all genomes. Leveraging the antibiograms of the isolates, they were compared with the software prediction to gauge software performance. Virulencefinder [76] was also run on all genomes. The detected genes which meet thresholds of 100% identity and same length of match and reference were arranged in a presence absence matrix. Subsequently, the virulence genes were classified by using the virulence factors database VFDB [77] and Uniprot [78] and by consulting relevant publications.

SCCmec cassettes

SCCmecfinder was also implemented for the genome set after dereplication using the SCCmecfinder web server [79].

R packages

The following R packages were used: pacman (v0.5.1) [80], tidyverse(v1.3.4) [81], rjson (v0.2.21) [82], ggtree (v3.2.1) [83], treeio (v1.18.1) [84], tidytree (v0.4.1) [85], ape (v5.6.2) [86], flextable (v0.8.2) [87], BactDating (v1.1.1) [73], fastbaps (1.0.8) [74], svglite (v2.1.0) [88], knitr (v1.4.0) [89].

Results

Epidemiological context of the CCoS

We collected and analyzed a legacy MRSA collection, mainly from Zurich, Switzerland (CCoS; n = 451, collected between 1965 and 1987), all publicly available genomes from a similar time-period (public repositories; n = 103, mainly from UK, Denmark, and Asia, collected between 1950 and 1992), and modern Swiss MRSA isolates sequenced at the University Hospital Basel (USB; n = 1207, deriving from different Swiss institutions, collected between 2009 and 2022). We characterized the genomes using MLST, cgMLST, and whole genome phylogeny. The complete dataset of 1761 genomes comprise 81 STs, of which 75 have been previously described (Fig. 1, Additional File 1A). Prior to 2000, 554 genomes were analyzed, belonging to 17 STs. Of these, six STs were previously undescribed and include 32 isolates. These six undescribed STs comprise two novel combinations of known loci (ST7844 and ST7842) and four exhibiting previously undescribed alleles (ST7638, ST7837, ST7839, and ST7843) (Fig. 1). In both historic datasets (Fig. 1A, C), most of the samples were collected between 1960 and 1972. The CCoS collection was dominated by ST247 isolates from the mid-1960s to the early 1970s, coinciding with low levels of ST250 and ST239 (Fig. 1C). A decade later, these STs were still present but we also registered ST254 and ST7844. The STs from public repositories of the same time period are similar (Fig. 1A). From 1980 onwards, the STs represented in the international collection are ST239, ST8, ST5, and ST22. Of note, all the previously mentioned STs, except for ST5 and ST22, belong to CC8, as seen in the MLST MST (Fig. 2). Among modern Swiss MRSAs, we observed many more STs with the most common being ST22, ST5, and ST8 (Fig. 1B). Other STs common in the CCoS are detected in small numbers among modern samples (ST239 and ST247). Among these contemporary isolates, ST228 MRSA is represented by a single sample isolated in Basel in 2021.

Fig. 1

Epidemic distribution of MLST sequence types by collection year and dataset (A Public repositories, B Culture Collection of Switzerland, C Contemporary Swiss MRSA, D Legend). Samples with no collection year are not represented, n = 6. Contemporary STs with less than 100 isolates are only shown as “other” (gray). Be aware that each subfigure has a different y-axis

Fig. 2

MLST minimum spanning tree of all MRSAs (historic Swiss, historic international, and modern Swiss) generated by Ridom SeqSphere+ (n = 1761). Nodes colored by ST with CC superimposed. CCs are defined as all STs which match their central genotype (ST) in four or more loci. Branch labels show number of allelic differences

Within CC8, ST250, ST254, and ST239 are single locus variants (SLVs) of ST8, while ST247 is a closely related SLV of ST250. The novel STs 7844 and 7842 are SLVs of ST239 and ST247, respectively. The results were visualized on a map, an interactive version can be found online under https://github.com/svannib/historic_MRSA/tree/main.

Identification of successful ancient lineages through cgMLST clustering

The average S. aureus genome contains 2872 coding sequences, of which 1861 are present in the S. aureus core genome MLST scheme. The cgMLST analysis of all ancient genomes (n = 554) yielded a cgMLST MST whose nodes have been colored by ST (Fig. 3), isolation country (Additional File 1, Figure B), and collection decade (Additional File 1, Figure C). Isolates of the same ST can be seen to generally cluster together. However, we observed that some STs exhibit higher diversity between isolates than others. For example, distances among ST247 isolates rarely exceed 30 alleles, while ST5 samples often lie more than 200 alleles apart. The isolates of these STs were sampled over very different time periods (ST247 mainly from 1963 to 1987 with three sample after 2000, while ST5 is more evenly distributed from 1980 to 2022, both with a gap in from mid-1990s to the mid-2000s), so this wide sampling could explain, at least in part, this higher diversity. Generally, STs with higher diversity contain samples isolated over a longer period. Another example is ST239, which shows a high diversity by cgMLST (allelic divergence between 50 and 120), with isolates from geographically distant countries (Australia, Singapore, United States, and Switzerland) and collected between 1972 and 2021. The genomes from the earliest ST239 (n = 5) are central to this group (Additional File 1, Figure B, highlighted in red). ST5 isolates are distant from other STs (1382 allele differences), displaying allelic distances up to 207 within this ST. ST22 samples are also distant from other groups (1599 allele differences), were isolated from different regions of Britain, and show allele differences of up to 37. The biggest allelic differences are seen between CCs, exceeding 1000 allele differences in all cases, while inside CCs distances never exceed 500 allelic differences. We identified seven putative transmission clusters, using a cluster cutoff of 24 allele differences (Fig. 3). Cluster 1 contains ST247 samples collected between 1960 and 1999, suggesting that isolates within this cluster remained successfully circulating throughout the decades. Cluster 2 also encompasses ST247 samples isolated in the 1960s from Switzerland and Denmark, where the first ST247 was identified. Cluster 3 includes both Swiss ( = 20) and British (n = 1) ST239 isolates, along with ST7844, the two STs separated by 11 allelic differences. Clusters 4 and 5 contain Swiss isolates mainly from the 1960s and 1970s, respectively. Cluster 6 contains British ST250 samples from the 1960s, representing the first described MRSA. Cluster 7 covers German ST254 whose closest relative is ST250 (203 allelic differences).

Fig. 3

cgMLST MST of CCoS and public repository genomes of MRSA (1960–1992, n = 554); nodes colored by sequence type and CCs are shaded. Clusters are shown with a maximal cluster distance of 24 allele differences and a minimal cluster size of 15. The earliest ST239 (1972) isolates are circled in red

Genotypic antimicrobial resistance (AMR) prediction and core genome phylogeny

To analyze the genomic data on a SNP rather than allele basis, core genome SNP phylogenies were calculated, calculated first on genomes from the CCoS and from public repositories sampled between 1960 and 1992 (n = 554; 2034 core genes present in over 99% of 554 isolates) (Additional File 1, Figure F) and subsequently with all MRSA covering a period from 1960 to 2022 (n = 1761; 1225 core genes across 777 genomes) (Additional File 1, Figure G). Both phylogenies show isolates from the same ST clustering together, as seen with cgMLST. We also observe the previously mentioned increased diversity among modern MRSA, with samples isolated after 2009 having a broader variety of STs. Prediction of AMR from the genome was limited to antibiotics that are important: (i) for classification of MRSAs (oxacillin); (ii) for the history of MRSA (penicillin); or (iii) current clinically relevant antimicrobials (tetracycline, erythromycin, clindamycin, gentamicin, and ciprofloxacin). The sensitivity of the methods, the ability to correctly predict AMR, and specificity, the ability to correctly predict susceptibility, were calculated by comparing genotypic predictions to available phenotypic data (Table 1). Both phenotypic and predicted antibiotic susceptibility data are displayed adjacent to the core genome phylogeny in Fig. 4 with a binary heatmap showing concordance between the two to facilitate interpretation.

Table 1 Sensitivity and specificity of AMR prediction for selected antibiotics and their respective antibiotic resistance encoding genes (ARG). The number of resistant and sensitive phenotypes given are based on the antibiograms (n = 451, EUCAST). “^*” indicates sample size of 1

Fig. 4

Predicted resistance/sensitivity, phenotypical resistance/sensitivity, and concordance between the two mapped to a maximum likelihood core genome SNP tree of ancient Swiss and international MRSA (1960–1992, n = 554). Leaves colored by sequence type. Outgroup (S. epidermis) line shortened through X-axis break for visualization purposes

Virulence genes

Across all datasets, 35 unique virulence encoding genes were identified, the most common coding for hemolysins and proteases, present in more than 90% of all isolates. Also common were staphylokinase (sak), staphylococcal complement inhibitor (scn), and toxin encoding genes such as lukD, lukE, sek, seq, and seb. Other staphylococcal enterotoxins (se_) were far less common, while the arginine catabolic mobile element (ACME) and toxic shock syndrome toxin (tst) were present in under 1% of samples. A virulence gene presence/absence heatmap was mapped to a phylogenetic core genome tree of the isolates (Fig. 5). The genomes show frequent co-occurrence of virulence factors sak and scn. Furthermore, international ST5 and ST22 isolates exhibit a general lack of serine-like proteases (splA/B/E) and toxins (lukD, lukE, seb, sek, seq) in favor of another group of staphylococcal enterotoxins (seg, sei, sem, sen, seo, seu, sec, sel). ST239 and ST7844 broadly lack seb while generally displaying sea presence. However, ST7844 and its most closely related Swiss ST239 group also lack sek and seq.

Fig. 5

Virulence gene presence/absence heatmap mapped to a maximum likelihood core genes phylogenetic tree of historic Swiss and international MRSA (n = 554). Leaves colored by sequence type. Outgroup (S. epidermis) line shortened through X-axis break for visualization purposes

SCCmec types

Among 777 dereplicated genomes, six unique SCCmec cassette types were detected, as well as several ambiguous predictions. For isolates which were predicted as methicillin-sensitive despite having a resistant phenotype (n = 5), no cassette was found in their genomes. SCCmec type I is by far the most prevalent (462/777, 59%), followed by type IV (144/777, 19%) and type III (61/777, 8%). SCCmec cassette types mapped to the core genome phylogeny (Additional File 1, Figure J) illustrate strong co-occurrence between cassette types and ST lineages.

International lineages shaped modern Swiss MRSA epidemiology

Based on cgMLST and core SNP trees, three groups were chosen for whole genome SNP analysis: ST239, ST7844, and closest STs (group 1); ST250 together with ST247 (group 2); and ST22 (group 3) . SNP phylogenies were generated using hybrid assembled Illumina/Nanopore reference genomes from within the cluster to best capture its diversity. The resulting phylogenic tree of ST239 (Additional File 1, Figure L) indicates that group 1 strains evolved from an ST239 ancestor, with SNP in MLST target genes resulting in SLV within the group. One example is a mutation in glpF which gave rise to ST7844 in a lineage which appears to previously have been successful, but which is not represented in modern Swiss MRSA. One cluster within this group contains Swiss, German, and British ST239 and ST7844, which were collected before 1990 (Bayesian cluster 3). It has no modern samples, hinting at its disappearance from Switzerland. Another cluster contains the earliest ST239 from Switzerland and international isolates from the 1980s and 1990s (cluster 4). Interestingly, some modern ST239 isolates (clusters 1–2) are closer to cluster 4 than to other modern samples. Cluster 5 covers both historic ST239 from Singapore and modern Swiss ST368 and ST241. The SNP phylogeny of ST247 and ST250 (group 2) shows that six Bayesian clusters of ST247 are present in the CCoS (Additional File 1, Figure M). Bayesian cluster 6 covers both Swiss and Danish samples collected in the 1960s. Modern ST247 play a small role in the modern Swiss epidemiology of MRSA, but these modern isolates cluster with historic Swiss, British, and Belgian ST250s. The ST22 phylogeny (group 3, Additional File 1, Figure N) again hints at the international origin of lineages circulating today in Switzerland, as British isolates from the 1980s are closely related to the ancestors of 85 modern Swiss samples (clusters 5 and 4).

Discussion

To put the CCoS in contemporary context, the high prevalence of ST247 among the CCoS sample set mirrors the epidemiological situation of other European countries at the time [9, 17, 90]. In addition, the Swiss ST239 from 1972 (n = 5) demonstrate that ST239 isolates were already present in Europe in the early 1970s and provides the earliest sequenced ST239 MRSA available [25]. Furthermore, ST228 MRSA, a sequence type with German origins endemic to Western Swiss hospitals [91], has been isolated only once in Basel, meaning this ST has not established itself in northwestern Switzerland.

Analyzing all MRSA in our collection, from 1960 to 2022 (Additional File 1, Figure E), we observe that CC8, dominant in 1960s–1970s, is not the main CC represented across modern samples. The only ST within CC8 which plays a main role after the turn of the century is ST8 (n = 141/1207, 12% of modern Swiss samples). This might be due to the international predominance of ST8 CA-MRSA in the last decades [92]. Particularly interesting is the strong modern presence of CC5 (n = 306/1207, 25% of modern Swiss samples) and CC22 (n = 210/1207, 17% of modern Swiss samples) isolates, which were present at the international (but not Swiss) level in the 1980s and 1990s. This suggests that international lineages arrived in Switzerland around the turn of the century. Furthermore, other complexes are present in modern Swiss hospitals such as CC1 (n = 85/1207, 7%), CC30 (n = 86/1207, 7%), and CC45 (n = 53/1207, 4%). In general, clones associated with HA-MRSA infections are the most prevalent both in ancient Swiss and public assemblies, possibly due to sampling bias, as most isolates were collected in hospitals. Modern global epidemiology shows less homogeneity, with lineages associated with community infections being more prevalent [93]. These changes in epidemiological landscape mirror the recent growth in MRSA biodiversity reported internationally [35, 39] and a more heterogenous sampling approach performed in the age of molecular microbiology.

The detection of antibiotic resistance genes (ARGs) of S. aureus to predict phenotypic resistance can reach specificity and sensitivity similar to routine susceptibility testing [94, 95]. In our case, sensitivity to clindamycin and erythromycin seem to be particularly challenging to predict (Table 1). As expected in an MRSA-only dataset, all isolates are resistant to meropenem, which was consequently removed from the visualization. However, four samples from the CCoS display an oxacillin-resistant phenotype despite no mecA or other mec variants being identified within the genomes. This could be due to beta-lactam resistance caused by mechanisms other than mecA, such as expression of penicillin-binding-proteins with low antibiotic affinity [96, 97] or overexpression beta-lactamases [98] or loss of the cassette after the phenotype was determined. Penicillin resistance (blaZ) is widespread in this dataset (99.8%, 450/451 penicillin-resistant phenotypes, 94.9%, 526/554 ancient genomes possess blaZ), being absent only from the majority of ST5 and some ST254 isolates (Fig. 4). MRSA possessing blaZ has been previously documented: since the late 1960s to this day, isolates are often resistant to penicillin [99, 100]. Also, among modern Swiss MRSA, blaZ is still present in the majority of isolates (89.5%, 1080/1207 genomes), a finding supported by other studies on European S. aureus [101]. Resistance to tetracycline is widespread (94.0%, 521/554) except for within ST5, ST22, and some ST250 genomes. Clindamycin resistance and erythromycin resistance are common among ST247 samples, except for ancient Danish genomes, present in some ST5 and ST22, and rare in ST250. Gentamicin resistance is a common feature among ST239, ST254, and ST7844, while being sparsely present throughout the rest of the dataset; ST22 and ST250 are uniformly gentamicin sensitive. Ciprofloxacin resistance is rare across the historic dataset and prediction of ciprofloxacin resistance fails to identify the single ciprofloxacin-resistant isolate (Table 1). The more recent introduction of the antibiotic on the market (1987) [102] and the multifactorial nature of ciprofloxacin resistance [95, 103] might explain the limited numbers of resistant isolates in ST239, ST5, and ST22. Higher incidence of ciprofloxacin resistance is registered in modern Swiss MRSA (37%, 452/1207). Where discrepancies were observed, the analysis was repeated on the sequenced reads, to exclude the possibility of assembly errors that could lead to an incorrect interpretation of the resistance. Of 44 AMR prediction discrepancies, two incorrect penicillin predictions and two incorrect gentamicin predictions were rectified. Furthermore, assuming that mecA also confers resistance to penicillin would correct eight further incorrectly classified isolates [104]. In addition, two isolates incorrectly predicted as clindamycin and erythromycin sensitive were found to possess the relevant genes which encode for the resistance (ermA, ermC) at lower coverage thresholds. Remaining discrepancies could be due to genetic instability (ermC) [105], errors in phenotypic testing and recording, or gaps in the knowledge base of resistance mechanisms. The data is listed in Additional File 3. This dataset suggests that ancient ST247 gained genes leading to a broader resistance profile (clindamycin, erythromycin, and tetracycline) compared to some early ST250. This might be one of the factors which contributed to the success of pandemic multidrug-resistant ST247 and to the decline of ST250, a hypothesis which was suggested by other publications investigating early MRSA epidemiology and historic genomes [6, 9, 106]. Furthermore, studies have shown that STs associated with broader resistance profiles tend to be more present in healthcare settings, where antimicrobial pressure is high and the fitness cost of resistance yields greater returns [107,108,109]. This balance between antibiotic resistance and fitness costs is another potential driver of success or failure of MRSA lineages. The early Danish ST247 samples in this dataset are erythromycin and clindamycin sensitive, but they may present a biased portrayal of isolates of the time, since widespread erythromycin and clindamycin resistance in ST247 isolates from the same period have been reported [17]. The phylogeny of the dereplicated dataset (Additional File 1, Figure H) suggests that modern Swiss MRSA have a higher variability of resistance patterns than ancient MRSA, even among closely related isolates. This could be due to the isolates coming from geographically separated Swiss hospitals or may be an artifact of the analytical dereplication performed on the modern Swiss MRSAs, as groups of very similar genomes are collapsed into one datapoint.

When analyzing virulence genes, strong co-occurrence of either serine-like proteases (splA/B/E) and toxins (lukD, lukE, seb, sek, seq) or other staphylococcal enterotoxins (seg, sei, sem, sen, seo, seu, sec, sel) may suggest the presence a genomic island. Previously described genomic islands which could play a role in the distribution and dissemination of the virulence genes seen in the dataset are vSAβ (sea, seg, sei, sem, sen, seo, splA, slpB, slpE, lukD, lukE), saPI3 (seb, sek, seq), and plB485 (sej, sed) [110]. The major toxin Panton-Valentine leukocidin (lukF-PV, lukS-PV) which heightens the virulence of MRSA [111] is sparsely present in modern Swiss MRSA but completely absent from our dataset prior to 2009 (Additional File 1, Figure I). This lies in contrast with the high rates of leukocidin reported in the US [112] and its presence in Switzerland dated to at least 1994 [113]. A limitation of this virulome analysis is its reliance on WGS and gene presence/absence instead of diagnostic tests. Still, this approach has shown high concordance with phenotype-based methods [114].

SCCmec I is strongly represented within ST247 and ST250 among historic isolates [27]. SCCmec type IV (2B) is common among historic and modern ST22, ST5, and ST8 isolates (Additional File 1, Figure K). Three of the five major ST8-related epidemic MRSA clones are represented in this dataset [27]. ST247-MRSA-I, ST250-MRSA-I (accounting for 80% (442/554) of the historic genomes), and ST8-MRSA-IV. Another strongly represented epidemic MRSA, although it rose to prominence later than ST250-MRSA-I and ST247-MRSA-I, is ST239-MRSA-III [27], which has the same cassette as the SLV ST7844-MRSA-III. In modern Swiss MRSA, SCCmec type IV in ST5, ST8, and ST22 isolates is the dominant cassette type and a marker for CA-MRSA [115]. These lineages have often been reported as dominant in many countries spanning the globe [116,117,118,119,120,121]. The wide presence of type IV in different MRSA lineages and other staphylococcal species might suggest an improved horizontal genetic transfer rate and/or lower fitness cost [40, 122, 123].

SNP phylogeny and Bayesian analysis of cgMLST clusters show multiple introductions of international lineages into German-speaking Switzerland, as these historic international strains are closely clustered to modern Swiss isolates by Bayesian analysis. Examples of this are historic ST239 MRSA from Singapore with Swiss ST368 and ST241, Historic Swiss, British, and Belgian ST250s with modern Swiss ST247 and Historic British ST22 together with modern Swiss ST22. Other clustered isolates suggest the historic presence of internationally introduced lineages which subsequently disappeared at both the Swiss and International level. This is the case for Swiss, German, and British ST239 closely related novel Swiss ST7844. Overall, Bayesian clustering provides evidence as to which of the old international MRSA lineages appeared briefly in Switzerland before being displaced, which contributed to the MRSA diversity we see in modern German-speaking Switzerland.

Conclusions

Since the 1960s, MRSA has been a challenging bacterial pathogen faced by clinicians worldwide. This study sheds light on the spread and relationships of major early MRSA clones. Our genome collection includes 451 MRSA samples from CCoS isolated between 1965 and 1987s in the greater Zurich area with a convenience sampling strategy, alongside 103 historic MRSA genomes from public repositories and 1207 modern MRSA isolated in Swiss-German hospitals. Despite being a sample set which is unrepresentative of the heterogenous MRSA epidemiology of Switzerland, our data reveal an historic epidemiological landscape in the investigated regions which is similar to that in the rest of Europe at the time. MRSA lineages which played an important role across European and Swiss-German hospitals from the 1960s to the 1990s, such as ST247-MRSA-I, ST250-MRSA-I, and the earliest ST239-MRSA-III are represented in the CCoS. Today, these clones appear to have been displaced in Switzerland, with international lineages from the last quarter of the twentieth century, including ST5-MRSA-IV, ST8-MRSA-IV, and ST22-MRSA-IV now being the major players in Swiss-German hospitals. Interestingly, we see little overlap between the different Swiss language regions, with major endemic clones such as ST228 MRSA in Geneva [91] being barely present in the German-speaking region. An analysis of the AMR and virulence profiles showed how different STs are associated with different AMR and virulence encoding genes. Discrepancies between phenotypes and genotype-based predictions were investigated but only a few could be resolved, and 9% (41/451) of the isolates lack genotype-phenotype concordance for one or more antibiotics. There are important limitations of our study. Firstly is the lack of sequenced Swiss isolates between 1988 and 2008, which prevents us from fully understanding the epidemiological changes which happened over the turn of the century. Secondly is the lack of phenotypic resistance data for the current Swiss MRSA, which would make our lineage-resistance association more robust and precise. The lack of both historical and contemporary isolates from other Swiss regions, notably French- and Italian-speaking Switzerland, which reported the highest MRSA prevalence means that the study is not necessarily representative on a national scale and cannot be integrated with older studies due to alternative typing methods used. Lastly, the lack of relevant clinical and epidemiological attributes limits the contextual understanding of their collection and interpretation in terms of molecular epidemiological investigation. Despite this, these data from CCoS hold high scientific interest, as the collection contains some of the first MRSA ever isolated and whose whole genome sequences we now present. We have thus significantly increased the public available genomes from the early period. The volume of isolates and phenotypic characterization make them an important addition to the pool of MRSA genomic data of samples isolated in the second half of the twentieth century.