Methicillin-resistant Staphylococcus aureus in China: a multicentre longitudinal study and whole-genome sequencing

ABSTRACT The aim of this study was to investigate the genomic epidemiology of MRSA in China to identify predominant lineages and their associated genomic and phenotypic characteristics. In this study, we conducted whole-genome sequencing on 565 MRSA isolates from 7 provinces and municipalities of China between 2014 and 2020. MRSA isolates were subjected to MLST, spa typing, SCCmec typing, analysis of virulence determinants and antimicrobial susceptibility testing. Among 565 MRSA isolates tested, clonal complex (CC) 59 (31.2%), CC5 (23.4%) and CC8 (13.63%) were the major lineages, and the clonal structure was dominated by ST59-t437-IV (14.9%), ST239-t030-III (6.4%) and ST5-t2460-II (6.0%), respectively. Of note, CC8, the predominant lineage in 2014–2015, was replaced by CC59 after 2016. Interestingly, the extension and unstable structure of the CC5 population was observed, with ST5-t311-II, ST764-t1084-II, ST5-t2460-II and ST764-t002-II existing complex competition. Further analysis revealed that virulence determinant profiles and antibiograms were closely associated with the clonal lineage. The CC59 MRSA was less resistant to most tested antimicrobials and carried fewer resistance determinants. But rifampicin resistance and mupirocin resistance were closely linked with CC8 and CC5, respectively. MRSA isolates conservatively carried multiple virulence genes involved in various functions. PVL encoding genes were more common in ST338, CC30, CC398, ST8 and CC22, while tsst-1 was associated with ST5. In conclusion, the community-associated CC59-ST59-t437-IV lineage was predominant in China, with diverse clonal isolates alternately circulating in various geographical locations. Our study highlights the need for MRSA surveillance in China to monitor changes in MRSA epidemiology.


Introduction
Methicillin-resistant Staphylococcus aureus (MRSA) is a notorious multidrug-resistant bacterium, which can cause a series of infectious diseases, such as septic shock, endocarditis, and severe pneumonia [1]. MRSA infections are complicated to manage since these pathogenic bacteria are resistant to multiple antibiotics. During the past few decades, the spread of MRSA, including several clonal lineages of hospital-associated MRSA (HA-MRSA) and communityassociated MRSA (CA-MRSA) and more recently livestock-associated MRSA(LA-MRSA) had been a major health problem worldwide [2]. ST239 and ST5 are the classic HA-MRSA clones prevalent in China and other Asia countries, while heterogenetic CA-MRSA clones such as ST59, ST338, ST30, ST72 and ST8 are disseminated in different regions [2,3].
Although there have been several epidemiological studies of MRSA in mainland China, they mainly focus on a subset of MRSA, such as specific infection type source or community origin [4,5]. However, MRSA infections due to lineages that were initially exclusively associated with CA-MRSA have been observed in hospital settings with increasing frequency, which may cause the overlap of the definition [2,6]. Knowledge of the entire MRSA population in hospital settings in China was still limited.
Analysing the genotypic characteristics of MRSA clones is valuable for understanding MRSA evolution and dissemination [1,7]. Whole-genome sequencing allows for a comprehensive analysis of the epidemiology and genomic repertoire of S. aureus and provides in-depth insights into the evolution of particular MRSA clones [6,8]. While few studies have characterized extensive collections of MRSA isolates from China using whole-genome sequencing [9,10].
Untargeted profiling of the entire MRSA population in a specific area is crucial in monitoring its emergence and spread and informing prevention strategies [11]. To produce a more comprehensive national description of the molecular epidemiology and phenotypic properties of MRSA in China, we characterized 565 MRSA isolates collected from several geographically dispersed Chinese hospitals over the past seven years.

MRSA isolates
For this study, non-duplicated MRSA clinical isolates were collected between 2014 and 2020 from 7 provinces and municipalities in China, including Shanghai, Zhejiang, Guangdong, Sichuan, Hubei, Jiangxi, and Inner Mongolia. These regions represented varying prevalence levels of MRSA in China [12] ( Figure  1). One representative tertiary care hospital was selected in each region for further investigation. In all the stored MRSA from 2004 to 2020, chosen randomly MRSA using random selection module in Microsoft excel. Information about these isolates was extracted from the laboratory information system. All isolates were reidentified species by the MALDI-TOF MS (Bruker Daltonics GmbH, Bremen, Germany), and cefoxitin resistance was determined by the disc diffusion method [13].

Genome sequencing and bioinformatic analyses
Total DNA from the MRSA isolates was extracted using QIAamp DNA Mini Kit (Qiagen) and sequenced on the Illumina NovaSeq platform using the 2 × 150-base pair paired-end mode. The raw reads were trimmed, and de novo assembled into contigs using CLC Genomics Workbench software (version 12.0; CLCbio). The average sequencing depth was 224.3 ± 57.92. An average of 4.34 ± 1.18M reads count was obtained from each MRSA isolate, and the GC rate was 32.72 ± 0.05. The contig number was 49.73 ± 59.57, with 217.9 ± 159.8Kbp of N50.
The assembled contigs were used for molecular typing, including MLST, spa-typing and SCC mec typing. STs were inferred using the S. aureus MLST database (https://pubmlst.org/organisms/staphylococcusaureus). Moreover, the spa and SCCmec-types were predicted using the Centre for Genomic Epidemiology website (https://cge.cbs.dtu.dk/services). In addition, the phylogenetic relationship of all MRSA genomes was analysed using recombination-filtered core SNPs by ParSNP [15].

Comparison of antimicrobial resistance pheno-and genotypes of MRSA major clonal complexes
The dominant CCs MRSA isolates, including 176 CC59, 132 CC5 and 77 CC8 isolates, displayed various antimicrobial resistance profiles, as shown in Figure 4 and Table S2,. Overall, CC59 displayed low resistance for most tested antimicrobial agents, except for erythromycin and clindamycin, as for the less content of numerous resistance determinants. Specifically, the resistance rate of the CC59 isolates to tetracycline was significantly lower than the CC5 and CC8 isolates (31.8% vs 78.0% and 58.4%, p = 0.000). Resistance to gentamycin, fusidic acid and mupirocin were rarely identified in the CC59 isolates (0.0%∼0.6%), whose genomes tended to lack aac (6')-Ie/aph (2'')-Ia, fusA_L461 K, fusB, mupA, ileS_V631F and ileS_V588F (Figure 4(A,C)). Moreover, unlike CC5 or CC8 with high cefoxitin resistance level (the most frequently observed MICs were ≥256 mg/L), the distribution of cefoxitin MICs of CC59 was mainly concentrated in 16 mg/L or 32 mg/L (Figure 4(B)). Notably, CC8 isolates tend to be resistant to rifampicin, of which resistance mostly conferring by rpoB_L466S and rpoB_H481N; CC5 isolates more likely to be resistant to mupirocin, and higher carriage rates of mupA and ileS mutations (Figure 4(A,C)).
Additionally, clonal replacement of predominant MRSA strains was observed in this study ( Figure 5 (A)). The CC8 lineage (mainly ST239-t030-III) was overwhelmingly predominant in 2014-2015 with up to 66.7% frequency, while subsequently losing its superiority and dramatically declining to 8.9% in 2020. Meanwhile, the frequency of CC59 (mainly ST59-t437-IV) and CC5 (mainly ST5-II and ST764-II) exhibited an escalating trend with from 27.8% in 2014-2015 to 35.6% in 2020 and 5.6% in 2014-2015 to 27.4% in 2020, respectively ( Figure 5(A,B)). Also it is worthy to notice that CC398, a notorious livestock-associated lineage in European countries, Australia, North and South America, increased from 0% in 2014-2015 to 8.9% in 2020 [19]. MIC of methicillin was defined as methicillin resistance. (C) The resistome analysis. Only resistance determinants found in more than 10 isolates and carried by CC59, CC5, or CC8 are specifically displayed. For more details, see Table S2.

Dynamic population structure of CC5
There was an association between CC5 and the high MRSA morbidity (Figure 1). In the current study, CC5, primarily composed of New York/Japan clone ST5 and its variant ST764, displayed a jagged increase trend over the past seven years ( Figure 5(A,B)). Interestingly, we also noticed that the multiple lineages in this CC population involved replacement events along with this trend, including ST5-t311-II, ST764-t1084-II, ST5-t2460-II and ST764-t002-II ( Figure 5  (B)). This temporal dynamic of heterogeneous clones indicated the potential for more complex competitive interactions in CC5. Moreover, although the geographical distribution of CC5 was scattered, ST5-II was more frequent in Hubei (45.2%, 28/62), and ST764-II exclusively prevalent in coastal cities, including Shanghai (47.8%, 22/46) and Guangdong (41.3%, 19/46) (Table S1).

Phylogenetic analysis based on core genome SNPs of STs
The phylogenetic trees of specific STs, including ST59, ST5, ST239, ST764 and ST398, were generated to investigate their heterogeneity further ( Figure 6). Interestingly, all these STs were clustered into two or more clades, indicating their high diversity. ST59 frequently transmit between sampled provinces, while ST239, ST5 and ST764 clones often spread locally. This was supported by the fact that ST59 clones from various provinces were interspersed in phylogenetic branches of ST59, while ST239, ST5 and ST764 within a single hospital setting or province tend to cluster into tightly linked clades.
We expanded these datasets to include genomes from publicly available worldwide collections ( Figure  S2). The global phylogenetic analysis of strains showed that ST59, ST5, and ST239 in China and European countries or North America tend to fall into different evolutionary clades. But, this pattern was not observed for the ST398 clone. China (n = 90), France (n = 20) and USA (n = 6) ST398 MRSA and MSSA strains clustered in the same clade, which also includes ST398 strains from Latin America and other European countries. Among them, the global pandemic ST5 clone represented the most heterogeneous cluster with at least fifteen clades. Notably, as a single-locus variant of ST5, ST764 has almost no publicly available sequenced genome, but it has established a local endemic in China.

Distribution of virulence determinants
In this study, 81 virulence factors were analysed for all these 565 isolates. Half of the analysed genes (40/81) associated with capsule, exoenzymes, secretion, iron uptake and metabolism, biofilm formation and haemolysis were highly conserved and prevalent in almost all strains (≥99%) ( Table S3).
The MRSA strains showed a diverse repertoire of virulence determinants across different lineages ( Figure 7). Typical human IEC (immune evasion cluster) was usually carried by ΦSa3 integrated into the hlb gene, and often a loss in host shift from human to animals [20,21]. IEC associated scn, sak, chp were detected in most isolates of STs in our cohort but absent in CC9 and ST630, indicating their potential livestock association. It's worth noting that the hlb and IEC (scn) were alternately identified in most isolates due to the insertion of ΦSa3 in hlb, except CC59 harboured both of them, which might indicate the uniqueness of IEC in CC59.
The incidence rates of toxin encoding genes were only 0.2-37.2%, while 446 out of 565 isolates possess at least one toxin gene (Figure 7 and Table S3). Regarding Panton-Valentine leucocidin (PVL), only 14.3% of isolates harboured pvl (lukF-PV and lukS-PV). Most CA-MRSA ST59 clones lacked pvl, while a high frequency of this gene was observed in ST338, ST8, ST1232, CC22, and CC30 (>66.7%). Toxic shock syndrome toxin encoding gene tsst-1 was identified in most ST5 isolates (96.8%) and some strains of ST508, ST1, and CC22 with a low frequency of 40.0-53.8%. Additionally, the adhesion associated with fnbA and sdr was more frequently carried by CC5 than CC59 and CC8. Exfoliative toxin encoding genes eta/etb were extremely rare and only detected in CC121. Additionally, patterns of enterotoxin genes are highly variable among MRSA isolated, while they are often associated with particular clonal lineage. Most STs of CC59 harboured seb-selk-selq, the majority of ST239 carried sea-selk-selq, while ST5 often had sec-sell. None of the enterotoxin encoding gene was found in ST950, CC2, CC9, CC398, CC30, CC121, ST630, ST188 and CC7.

Discussion
Our study aimed to gain epidemiology data on the genetic and phenotypic traits of MRSA circulating in China over a period of seven years. This study presents an extensive collection of 565 MRSA isolates collected by seven provinces and municipalities covering a substantial geographical zone of the country.
In this study, CC59-ST59, CC5-ST5/ST764 and CC8-ST239 were the key lineages responsible for mediating the development of the methicillin resistance in S. aureus in clinical environment in China [9]. The HA-MRSA CC8-ST239 has been a pandemic linage circulating in many countries worldwide since the 1970s [22]. Previous studies and our cohort identified ST59 as the major MRSA clones across China for decades [4,[23][24][25][26][27]. However, we confirmed the clonal shift of MRSA has occurred in Chinese hospital settings. After 2015, CC8-ST239 has lost its predominance and was replaced by CC59-ST59 and CC5. Not coincidentally, this similar clonal shift was also documented in America, Malaysia, Singaporean, Hungary, Portuguese, Czech Republic and some Latin American Countries in recent decades [2,[28][29][30], indicating it appears to be a common phenomenon.
The finding of this work confirmed that ST59-IV became the most dominant lineage in Chinese hospitals in recent years. ST59 was one of the most successful and persistent CA-MRSA clones in Asia, initially evolved independently of hospital strains and associated with fatal infections of outside hospitals [31,32]. Our data supported previous analyses that CA-MRSA clones have switched their epidemiological niches and are the most dominant lineage in nosocomial settings [33][34][35].
Most CC59 in the present study lacked pvl and exhibited relatively low resistance to common antimicrobials [36,37]. Researchers have proved that ST59 with fewer resistance fitness costs displays more competitive, which might facilitate CC8-ST239 replaced by CC59-ST59 [38]. Moreover, several CA-MRSA clones have smaller cassette chromosome mec elements such as SCCmec IV and SCCmec V, which provide only low-level methicillin resistance, supporting the idea that the extended evolution of CA-MRSA lineages is generally characterized by a very low methicillin resistance level, instead of harvesting more virulence determinants [39].
Another important finding in our study is the complex competition in the CC5 population structure. We found that multiple CC5 clones, including ST5-t311-II, ST5-t2460-II, ST764-t1084-II and ST764-t002-II, were alternatively predominant many times over the past seven years. Importantly, we observed that multi-drug resistant CC5 seems to be associated with high MRSA morbidity, and most of CC5 MRSA isolates were recovered from sputum, but most epidemiological studies were solely focused on the bacteraemia associated MRSA [23,40,41]. Therefore, we speculated that the importance of CC5 in clinical settings might be underestimated. Especially, the Figure 7. Heatmap illustrating virulence genotype of different STs (CCs). The evolutionary tree based on MLST was generated in MEGA X using the Maximum Likelihood method and Tamura-Nei model with 1000 bootstrap replicates. These virulence factor genes with more than 99% carriage were not displayed in this figure. For more details, see Table S3. refined temporal analysis of the different clones in our study showed that, with CC239-ST239 clones losing their dominance, the most obvious change is the marked increase in the incidence of CC5. We also found that CC5 harboured more adhesion-associated genes (fnbA and sdr) and superantigenic toxin gene tsst-1. However, it is still unclear whether the CC5 with an unstable clonal structure plays some potential role in the clonal replacement of CC8-ST239 in China.
Our analysis also highlights the growing public health threat of LA-MRSA clones. CC398 is associated with various animals, particularly pigs across European countries and North America. Occasionally, it can cause severe infections in humans. In our study, the incidence of CC398 persistently increased over time and displayed a high heterogeneous in the global geographical phylogenetic evaluation, suggesting that the LA-MRSA ST398 clone might be introduced in China and gain further adapt to the clinical environment.
CC9 was the livestock-associated MRSA (LA-MRSA) that prevailed prosperously in livestock farms and animals in China and other Asian nations [42]. And ST630 was also reported associated with bovine mastitis in China [43]. In our study, four ST9 and nineteen ST630 clinical isolates all had an intact hlb gene and no IEC genes, supported their animal origin [44,45].
Consistent with previous reports, most MRSA isolates in our study displayed multiple-drug resistance [46]. Fortunately, we confirmed that several antibiotics, such as vancomycin, teicoplanin, dalbavancin, quinupristin-dalfopristin and ceftaroline, still possessed superior antimicrobial activity and may be available in selections for anti-MRSA infections. However, intermediated resistance to vancomycin has sporadically been reported in China [38,47], alerting potential challenges in antibiotic treatment. Furthermore, our data indicated that the CC59 tends to be susceptible to multiple antibiotics in vitro, while the CC8 exclusively enriches the resistance of rifampicin. According to the surveillance data from the CHINET, resistance to quinolones, aminoglycosides and rifamycin in MRSA isolates was dramatically decreased in China recent years [48]. Our findings suggested that this shift of resistance phenotype is largely a result of clonal displacement (from CC8-ST239 to CC59-ST59).
The producer of many toxins collectively contributes to S. aureus virulence potential. Most of the isolates in our study carried at least one toxin gene. pvl was only carried by a small proportion of isolates, and there was no significant association with CA-MRSA lineage ST59 in our data. The presence of PVL has previously been strongly associated with CA-MRSA in many studies [49,50], while the findings of our study do not support this notion. Most of ST338, CC30, CC398, ST8 and CC22 carried pvl in the present study, demonstrating the pathogenic potential of these MRSA lineages. Although unusual in other STs, 96.8% of ST5 isolates carried tsst-1 gene in this study, which is consistent with the report from Suzhou, with a similarly high positive percentage (96.3%) [51]. Additionally, they emphasized that tsst-1 positive CC5 isolates were associated with higher mortality rates.
The limitation of this study is the absence of adequate clinical information, makes us fail to analyse the types of infection, the mortality, the baseline of patients and association with specific lineages. This is the first longitudinal large-scale surveillance of MRSA in China, covering 7 provinces and municipalities, providing important insight into the development of MRSA. Our major findings were as follows: Although MRSA exhibiting a range of strain types were detected in China, CC59, CC5 and CC8 lineages are responsible for the spread of MRSA in China. The epidemiology of MRSA varied temporally and geographically, while CC59-ST59-t427-IV lineage continually dominance in China hospitals in recent years. Different lineages are markedly associated with specific antibiotic resistance profiles and virulence patterns.