Methicillin-resistant Staphylococcus aureus in Saudi Arabia: genomic evidence of recent clonal expansion and plasmid-driven resistance dissemination

Abstract

Objectives: Staphylococcus aureus is a leading cause of hospital-acquired infections worldwide. Over recent decades, methicillin-resistant Staphylococcus aureus (MRSA), which is resistant to multiple antimicrobials, has emerged as a significant pathogenic strain in both hospital and community settings. The rapid emergence and dissemination of MRSA clones are driven by a dynamic and evolving population, spreading swiftly across regions on epidemiological time scales. Despite the vast geographical expanse and diverse demographics of the Kingdom of Saudi Arabia and the broader West Asia region, the population diversity of MRSA in hospitals in these areas remains underexplored.

Methods: We conducted a large-scale genomic analysis of a systematic Staphylococcus aureus collection obtained from 34 hospitals across all provinces of KSA, from diverse body sites between 2022 and 2024. The dataset comprised 581 MRSA and 31 methicillin-susceptible Staphylococcus aureus (MSSA) isolates, all subjected to whole-genome sequencing. A combination of phylogenetic and population genomics approaches was utilized to analyze the genomic data. Hybrid sequencing approach was employed to retrieve the complete plasmid content.

Results: The population displayed remarkable diversity, comprising 48 distinct sequence types (STs), with the majority harboring community-associated SCCmec loci (types IVa, V/VII, and VI). Virulence factors associated with community-acquired MRSA (CA-MRSA), including Panton-Valentine Leukocidin (PVL) genes, were identified in 12 distinct STs. Dominant clones, including ST8-t008 (USA300), ST88-t690, ST672-t3841, ST6-t304, and ST5-t311, were associated with infections at various body sites and were widely disseminated across the country. Linezolid and vancomycin resistance were mediated by cfr-carrying plasmids and mutations in the vraR gene (involved in cell-wall stress response) and the murF gene (involved in peptidoglycan biosynthesis) in five isolates, respectively. Phylodynamic analysis revealed rapid expansion of the dominant clones, with their emergence estimated to have occurred 10-20 years ago. Plasmidome analysis uncovered a diverse repertoire of blaZ-containing plasmids and the sharing of erm(C)-encoding plasmids among major clades. The acquisition of plasmids coincided with clonal expansion.

Conclusions: Our results highlight the recent concurrent expansion and geographical dissemination of CA-MRSA clones across hospitals. These findings also underscore the interplay between clonal spread and horizontal gene transfer in shaping the resistance landscape of MRSA.

Keywords: Kingdom of Saudi Arabia; MRSA; MSSA; ST; bloodstream infection; precision epidemiology.

Figure 1

The distribution of major ST clones (STs with more than five representative isolates) across sources and provinces. (A) Each pie chart represents the frequencies of major STs in the provinces of the KSA, with the chart size corresponding to the collection size from each province. (B) The relative (top) and absolute (bottom) frequencies of STs across body sites of isolation for major STs. Asterisks (*) and (**) indicate significance levels of < 0.05 and < 0.01, respectively, based on a one-sided proportion test comparing the frequencies of blood (red) and wound (yellow) isolates. (C) Population diversity of the collection across provinces. Each boxplot represents pairwise SNP distances between genomes within a province. The ** symbol denotes a significance level of < 0.01 based on a one-sided Wilcoxon rank test, comparing the SNP mean distances of each province to those of the other provinces.

Figure 2

The frequency of resistance determinants (antimicrobial resistance genes/mutations) counts across the major STs. Each boxplot in (A) denotes the count of resistance determinants identified by AMRFinder. The ** and *** signs in (A) corresponds to the significance levels of < 0.01 and < 0.001, respectively, from the one-sided proportion test, indicating whether the mean of the frequency of resistance determinants count were higher in the ST clone compared to the rest of the collection. The green circles on top genes in (B) show the determinants significantly linked with the resistance phenotype from the GWAS analysis. The gray squares in (B) shows genes/mutations which had a significantly higher frequency in each ST compared to the rest of the collection (p-value < 0.01 from one-sided proportion test).

Figure 3

The relative frequencies of the spa types SCCmec types across the major STs. Spa types with a relative frequency of < 0.1 for each clone were aggregated. Bars corresponding to the sub-clones that have at least one isolate carrying the PVL and tsst-1 genes are marked.

Figure 4

Phylodynamic analysis of clones with the most prevalent spa type. (A) The estimated clock rate and age values of the clones from Bayesian analysis. The error bars are the credible intervals [95% Highest Posterior Density (HPD)]. (B) The skyline growth for the clones in (A), based on the Bayesian tree for the clones. The shaded area corresponds to 95% confidence interval.

Figure 5

The phylogenetic trees for the clones with key resistance determinants, virulence factor genes, SCCmec type, antimicrobial resistance phenotype and the presence of key plasmids for the ST8-t008 clones in (A) and other expanding clones (B). Only resistance determinants, virulence gene factors and antimicrobial resistance phenotypes that were variably present in each clone according to Figure 2, Supplementary Figures S1, S3 are shown. The plasmid band shows the coverage percentage after mapping the short reads against the plasmid with the resistance genes retrieved from long-read sequencing data. Color codes in (B) are the same as in (A). The genomic map of the plasmids is provided in Supplementary Figure S2.