Genomic diversity of prevalent Staphylococcus epidermidis multidrug-resistant strains isolated from a Children’s Hospital in México City in an eight-years survey

Ethical considerations

This study was carried out following the recommendations of the ethics review committee of the Facultad de Medicina-UNAM. The SE strains used in this work were obtained by a donation from the microbiology collection of the Instituto Nacional de Perinatología “Isidro Espinosa de los Reyes” (INPer). A consent form was not required. Original identification keys and clinical data concerning the isolates are maintained under the control of INPer. In the present work, new strain identifiers were assigned to the INPer strains. Authors do not have access in any form to the specific clinical information of strains and patients.

Clinical isolation of staphylococci and characterization

Staphylococci strains were recovered from primus isolates taken from patients at INPer. They were conserved frozen (−80  °C in an ILSHIN ultra-freezer). Characterization of staphylococci was carried out by streaking the cultures on tryptic soy agar plates (BD Diagnostic Systems, Germany). Plates were incubated for 18–24 h at 37 ° C to obtain isolated colonies. Individual colonies were used to inoculate tryptic soy broth and grown overnight after which genomic DNA was extracted (see below). Identification of Staphylococcus species was carried out using the VITEK^®2 equipment (bioMérieux SA, 376 Chemin de l‘ Ome. France). Briefly, staphylococci colonies were resuspended in 0.45% saline solution and adjusted to a turbidity of 0.5 in a McFarland nephelometer. The bacterial suspension was transferred into the VITEK^®2 GP ID card testing 64 biochemical properties for gram-positive bacteria, and eight specific tests for SE species (phosphatase and urease production, growth in 6% NaCl, resistance to novobiocin and polymyxin B, arginine hydrolysis, catalase, D-mannitol fermentation). Additional growth characteristics, biochemical (D-mannitol fermentation, coagulase, and catalase), and molecular tests (PCR detection of coa and mecA; see below) were performed at the arrival of the staphylococci collection to the Pathogenicity Laboratory at Faculty of Medicine, UNAM.

Molecular tests

To confirm the presence/absence of the coagulase gene (coa), and the methicillin-coding gene (mecA) in SE and S. aureus strains, we did a duplex PCR according to Hookey (Hookey et al., 1999). In brief, a PCR amplification of a 177 bp fragment of the coa gene, with primers coa-forward 5′-AACTTGAAAT-AAAACCACAAGG-3′, and coa-reverse 5′TACCTGTACCAGCATCTCTA-3′. In the same reaction, a 458 bp fragment of the mecA gene was amplified using primers mecA-forward 5′-ATGGCAAAGATATTCAACTAAC-3′and mecA-reverse 5′-GAGTGCTACTCTAGCAAAGA-3′. PCR reactions were performed on a BIO-RAD C1000 Thermal Cycler (BIO-RAD, USA). The amplification conditions were as follows: 94 °C for 10 min and 30 cycles at 94 °C for 30 s, 56.2 °C for 30 s and 72 °C for 30 s and a final amplification at 72 °C for 10 min; PCR products were separated by electrophoresis in a 1% agarose gel stained with ethidium bromide and visualized with UV light.

Antibiotic resistance

The antibiotypes for 17 antibiotics were determined using automatic VITEK^®2 equipment with panel card AST-GP577 (VITEK^®2 Biomerieux, Francia). Resistance and susceptibility patterns to 17 antimicrobials were determined according to the M100S Performance Standards for Antimicrobial Susceptibility Testing (Clinical and Laboratory Standards Institute C, 2019). The following antibiotics were tested using the respective reference MIC: Amoxicillin/Clavulanic Acid: ≤4/2 µg/mL, Ampicillin: >8 µg/mL, Cefazolin: ≤2 µg/mL, Ciprofloxacin: ≥4 µg/mL, Clindamycin: >4 µg/mL, Erythromycin: >8 µg/mL, Gentamicin: >16 µg/mL, Imipenem: 8 µg/mL, Levofloxacin: ≥4 µg/mL, Linezolid: ≥8 µg/mL, Oxacillin: ≥0.5 µg/mL, Penicillin: ≥0.25 µg/mL, Rifampicin: ≥4 µg/mL, Synercid: > µg/mL, Tetracycline: ≥16 µg/mL, Trimethoprim/Sulfamethoxazole: >4/76 µg/mL, and Vancomycin: ≥32 µg/mL.

Methicillin resistance was evaluated by the disc diffusion method for cefoxitin on Mueller-Hinton agar plates (OXOID Ltd, Basingstoke, Hampshire, England), and with the plate method in Mueller-Hinton agar plates supplemented with 6 µg/ml of oxacillin, according to the standards from the Center for Disease Control and Prevention, USA (Centers for Disease Control and Prevention C, 2019)

Genome sequencing

Genomic DNA was isolated by using the Lysostaphin-Lysozyme method of Hookey (Hookey et al., 1999), only modifying the lysis step. Briefly, bacterial cells were resuspended in 250 µl of lysis solution (lysozyme 50 mg/ml; lysostaphyn 1 mg/ml in 25 nM Tris-HCl- 10 mM EDTA; all reagents from SIGMA-ALDRICH, St. Louis MO USA), and incubated at 37  °C by 90 min. The sample was then warmed up at 95  °C for 10 min and centrifuged at 8,000× g, after which the sample was processed as indicated in the method.

DNA integrity and purity were assessed by agarose gel electrophoresis and spectrophotometry (NanoDrop 2000; ThermoFisher, Waltham, MA, USA); final concentrations were assessed using the Qubit high sensitivity dsDNA assay (ThermoFisher). Libraries were generated with Nextera XT DNA library preparation kit (Illumina, San Diego, CA, USA), according to the manufacturer’s protocol, using 1 ng of DNA. Libraries quality were evaluated on an Agilent 4200 TapeStation system (Santa Clara, CA, USA) and sequenced using a NextSeq 500/550 mid output kit (2 × 150 bp) (Illumina), at an estimated depth of coverage of 100x.

Genome assembly and annotation

Assemblies were made with Spades 3.6.0 in genomes sequenced in high coverage (50x to 100x). Genome annotations were obtained from the PATRIC server (Wattam et al., 2018) (https://www.patricbrc.org/), through the automated bioinformatic method RAST (Rapid Annotation Subsystem Technologies) (Antonopoulos et al., 2019). Gene annotation of antibiotic resistance and virulence-related genes were obtained from the section of Special Genes of PATRIC. Annotations for virulome were supported with other sources of information such as PATRIC VF, VFDB, and Victors. For the prediction of antibiotic resistance genes, the PATRIC server was used to search specialty genes using the databases CARD (Comprehensive Antibiotic Resistance Database) and NDARO (National Database of Antibiotic Resistant Organism) (Davis et al., 2016).

Genome comparisons and pangenome modeling

Average Nucleotide Identity by MUMmer (ANIm) and Genomic Coverage (G _cov) were calculated with the JSpecies program, with MUMmer used as a pairwise comparison tool for pairs of SE genomes (Richter et al., 2016). The pangenome was modeled using the GET_HOMOLOGUES package and its default values (Contreras-Moreira & Vinuesa, 2013). Briefly, groups of orthologues proteins were computed using the Ortho-MCL integrated into GET_HOMOLOGUES. Accessory genes (unique genes and genes present in at least in two genomes but not in all genomes), were obtained from the cloud, shell, and soft-core pangenome components according to GET_HOMOLOGUES (Contreras-Moreira & Vinuesa, 2013). The distribution of accessory genes in SE genomes was illustrated with the heatmap.2 function of the R’s ggplot2 package.

Phylogeny

A consensus core genome calculated from the pangenome model of 29 SE genomes was obtained with GET_HOMOLOGUES (Contreras-Moreira & Vinuesa, 2013). The resulting 1575 core protein clusters were subject to multiple alignments with MUSCLE (Edgar, 2004), gaps removed with TrimAl v2.1, and were concatenated using homemade Perl scripts (Capella-Gutiérrez & Gavaldón, 2013). Phylogenetic trees were constructed by the maximum likelihood (ML) method based on the substitution matrix of WAG (Whelan and Goldman model), with 1000 bootstrap replicates, using RaxML program (Stamatakis, 2015). To draw and edit the phylogenetic trees, we used the iTool program (Letunic & Bork, 2019).

Recombination

The inference of recombination was performed with ClonalFrameML (Didelot & Wilson, 2015). First, we ran GET-HOMOLOGUES to obtain the common protein clusters encoded by the genomes of each set of SE strains (Contreras-Moreira & Vinuesa, 2013). Second, they were converted to nucleotide sequences and concatenated using homemade Perl scripts. Multiple alignments were made as described above for phylogeny construction. Third, a RAxML (Stamatakis, 2015) tree was done to obtain the Newick format and the transition/transversion parameter κ for running ClonalFrameML under default parameters (Didelot & Wilson, 2015). Z-score statistics were obtained for all the sets of SE genomes and p-values using the web application Z Score Calculator (https://www.socscistatistics.com/tests/ztest/zscorecalculator.aspx). BoxPlots were performed with the R ggplot2 system.

Mobile elements identification

Prophages were identified with PHAST (Arndt et al., 2019). Only predictions ranked as “intact prophages” were considered for analysis. IS, CRISPR-Cas elements and spacer sequences were obtained from PATRIC annotations (Antonopoulos et al., 2019). Then, IS were classified into families by BLASTx comparison with the ISfinder (Siguier et al., 2006). The phage identity of the spacers within the CRISPR-Cas elements was determined by BLASTn comparisons with the NCBI virus database. Only identical matches with a phage sequence in the database were recorded (Páez-Espino et al., 2016).

Identification and classification of SCC mec and ACME proteins

All the predicted proteins that compose the 12 types of SCCmec and the six ACME types were downloaded from GenBank (Table S1) and compiled in a local protein database (SE-17db). Then, all encoded proteins of the 17 SE genomes were used as a query in BlastP similarity searches with the SE-17db. Best-blast hits with a minimum similarity of 80% and 60% of coverture over the length of the smaller protein were taken as evidence of homology with the respective SCCmec or ACME type. To classify the SCCmec types, we used the phylogenies of the recombinases ccr A, B, and C, performed according to the method described above. ACME types were classified by the presence/absence of the genes of the operons arg, opp, ars, and kpd according to a reference table (Table S2).

GenBank accession numbers

S. epidermidis of the INPer collection used in this work were uploaded in GenBank with the following Biosample identifiers: SAMN11086744, SAMN11086745, SAMN11086746, SAMN11086747, SAMN11086748, SAMN11086749, SAMN11086750, SAMN11086751, SAMN11086752, SAMN11086753, SAMN11086754, SAMN11086755, SAMN11086756, SAMN11086757, SAMN11086758, SAMN11086759, SAMN11086760. The accession numbers for the genomes of reference S. epidermidis strains are listed in Table S3.

Survey of antibiotic multidrug-resistant intrahospital staphylococci in eight years period

We studied the incidence of staphylococci species and strains multi-resistant to antibiotics in the INPer children hospital from 2006 to 2013. A total of 822 staphylococci strains were recovered from distinct infection sites of newborns and adults. Staphylococcus species were identified using the VITEK^®2 system and standard clinical methods, including tests for biofilm formation, coagulase reaction, presence of coa and mecA gene, and resistance to 17 antibiotics (see methods). The analysis showed that SE strains were the most abundant (573 strains), followed by S. aureus (146 strains), and other coagulase-negative Staphylococcus (CoNS) species in low proportion (81 strains) (Fig. S1). The SE group was composed mainly of strains recovered from blood samples and catheters (Fig. S1).

The number of antibiotic multi-resistant SE strains per year showed similar profiles through the studied period (Figs. 1A–1H). Methicillin resistance was found in about 80% of the SE isolates. SE strains multi-resistant to nine antibiotics were the most frequently found throughout the eight years. In contrast, S. aureus strains showed resistance to about four antibiotics as the most common profile (Fig. S2). Similarly, the resistance to each antibiotic remained without change, during the eight-year study (Figs. 1I–1P). These results suggest the persistence of a multidrug-resistant clone or few clones of SE at the hospital during the eight years, or frequent horizontal exchange of genes encoding antibiotic resistance between SE clones. To study these alternatives, we choose one or more SE strains per year to totalize 17 clinical SE coagulase-negative (CoNS) strains (Fig. 1). These particular strains came from nosocomial infections of 14 newborns and three adults, isolated from three different infection sites: blood (8 strains), catheters (7 strains), cerebrospinal fluid (1 strain), and soft tissue (1 strain) (Table 1).

Broad gene catalog from draft genomes

We obtained the whole genome sequence of 17 nosocomial SE strains (see ‘Material & Methods’). After testing different parameters with the assembler programs Velvet and Spades (Bankevich et al., 2012; Zerbino, 2010), we got draft genomes assemblies consisting of about 92 up to 432 contigs with a 60–70×average sequence coverture per genome (Table S4). To assert that the 17 assemblies represent a substantial part of the SE genomes, we compared the total genome length, and the number and length of the predicted ORFs, with 11 complete genomes of SE downloaded from GenBank (Table S5). There were no differences between the genome length of the SE INPer genomes and the complete genomes from the GenBank (unpaired t-test = 2.33; p-value = 0.022), in ORFs number (unpaired t-test = −1.68; p-value = 0.117), or ORFs length (unpaired t-test = 2.60; p-value = 0.014) (Table S5). Therefore, we can conclude that the SE INPer genome sequences obtained here provide a broad catalog of genes per genome useful for comparative genomics.

Genomic, pangenomic, and phylogenetic relationships among SE isolates

To define the genomic similarity between the 17 clinical SE strains, we performed whole-genome pairwise nucleotide identity estimates (Average Nucleotide Identity by Mummer, ANIm) (Richter et al., 2016). The 17 SE showed high genomic ANIm values of about 99%, covering more than 90% of the genome length (Fig. S3). One exception is strain S10, which showed ANIm values of about 97% concerning the rest of the strains.

Besides the genomic identity between SE isolates, we wanted to investigate the extent of their genetic variability by performing a pangenome model. To this end, the SE collection was complemented by the inclusion of 12 complete genomes of SE strains downloaded from GenBank (Table S3, NCBI reference complete genomes). The model obtained using the GET_HOMOLOGUES software package (Contreras-Moreira & Vinuesa, 2013), indicated an open pangenome (Fig. S4). The core genome component for the 29 SE strains was predicted to consist of about 1,575 gene clusters, whereas the sum of genes unevenly distributed in the 29 SE genomes (accessory component) contains 4,360 gene clusters.

To know the phylogenetic relationship of the SE from INPer in the context of reference SE strains, we did an un-rooted ML phylogenetic tree using the predicted 1,575 concatenated core proteins (Fig. 2). There were three clades separated by the largest branches in the tree that comprise most of the SE local strains, and one or more SE strains isolated worldwide (Fig. 2, clades B, C, D). The clade marked as D in the tree consisted of two groups, one of which includes only reference SE strains, whereas the other had most of the SE strains of the collections studied. Strain S10, the most different strain by ANIm in the SE collection, was grouped in the clade A with the commensal strains SE ATCC12228 (2) and SE 14.1.R1 isolated from USA and Denmark.

The distribution of accessory genes in individual SE strains coincided with the phylogenetic clades. The SE strains grouped in three phylogenetic clades are closely related by their similitude in the gene presence/absence profile (Fig. 3). Despite the high identity of the SE strains, there is still considerable individual variation that may account for adaptations to local milieu.

Clonal structure

To investigate to which clonal ST complex the SE strains belong, we looked for the seven proteins of the S. epidermidis MLST scheme, and compared them with their respective alleles in the Staphylococcus epidermidis MLST database (https://pubmlst.org/sepidermidis/; Table 1; see methods) (Feil et al., 2004; Thomas et al., 2007). The analysis showed a total of 8 different STs; seven of them already recorded in the database. The S10 strain had an unassigned ST in the database and only differed by a single amino acid substitution in the YqiL protein (L370C). The ST2, ST5, and ST23 are worldwide distributed and are the most represented in our sample (Lee et al., 2018; Miragaia et al., 2007). In agreement with global data, ST35, ST59, ST81, ST89 are less frequently represented. The clonal relationships among STs determined by eBURST, indicate that founder clones are ST2 and ST5, whereas the other four STs (ST 59, 81, 89 and 23), are peripheral clones mostly related to each than to the primary founder clones (Fig. S5).

Virulence genes

Several known virulence genes of SE are shared by commensal and pathogenic strains (Otto, 2009). Among them, we found the cluster icaADBCR (biofilm formation) in nine out of the 17 SE genomes analyzed, except for icaC, absent in strain S05 (Table 1). The phenol-soluble modulins (PSMs), involved in inflammatory response and lysis of leukocytes were present in the genomes of all SE strains (Peschel & Otto, 2013). The novel ESAT 6 (esaAB, essABC, and esxAB genes) secretion system implicated in immune system evasion and neutrophil elimination was only present in S10 strain (Burts et al., 2005; Dai et al., 2017; Wang et al., 2016).

Frequently, SE isolates contain the genomic island ACME, composed by arginine deaminase catabolic (arg) and an oligopeptide permease operons (opp) (Miragaia et al., 2009). ACME has been associated with successful adaptation to the skin and mucosal surfaces outcompeting other related bacteria (Planet et al., 2013). The diversity of ACME has been organized in six distinct types according to their gene composition (Table S2 and references therein). We looked for the presence of ACME genes in the 17 SE genomes of INPer collection. By means of BlastP comparisons we identified genes belonging to the ACME Type I (arg⁺ opp⁺) in genomes of SE strains S03, S17, and S21, and the ACME Type V that contain the arg⁺ opp^+, and ars⁺ (arsenate resistance operon), and kpd+, a potassium transporting operon in S24 (Centers for Disease Control and Prevention C, 2019; O’Connor, McManus & Coleman, 2018) (Table S2). Even though some other SE strains present genes of the arg and ars operon, they lacked the genes of the opp operon; therefore, they were not assigned to ACME known types.

Antibiotic resistance genotype and phenotype

The SE strains were tested for their susceptibility to methicillin and other β-lactams antibiotics as described in methods. All the 17 SE strains have the β-lactamase gene (blaZ) and their regulators (blaR and I), which are probably responsible for the broad resistance spectrum to penicillin, carbapenems, and cephalosporins determined by the VITEK^®2 system (Table 2). Methicillin resistance was evaluated by the disc diffusion method for cefoxitin and corroborated by the presence of mecA in the genomes (Table 1 and Table 2). 12 out of 17 SE strains showed agreement between the cefoxitin resistance phenotype and the mecA genotype. Although strains S03, S05, S14, S15, and S21, have the mecA gene in the genome, they were scored as susceptible for cefoxitin (disc diffusion ≥ 25 millimeters). To support that the 15 mecA positive strains are Methicillin-Resistance SE (MRSE) strains, we evaluated their resistance to oxacillin. Commonly, the Methicillin-Resistance S. aureus (MRSA) strains are evaluated by the resistance to oxacillin (Centers for Disease Control and Prevention C, 2019). In the SE collection, only the S07 strain was susceptible to oxacillin (6 µg/ml) in plate assays, despite the detection of mecA in the genome. Altogether, these results indicate that most of the SE strains (15) studied here are MRSE, identified by the presence of the mecA gene; any phenotypic inconsistencies could be due to lack of expression of the mecA, heteroresistance, or the limitation of the current antibiotic susceptibility testing methods (Band et al., 2019; Harrison et al., 2019; Nicoloff et al., 2019).

The gene mecA encodes for a penicillin-binding protein (PBP) carried in a mobile element known as Staphylococcal Chromosomal Cassette or SCCmec (International Working Group on the Classification of Staphylococcal Cassette Chromosome E, 2009). By localizing the recombinases ccrA, B, and C, as well the mecA genes in the contigs of the respective genomes, and then constructing phylogenies including known SCCmec recombinase genes, we classified the SCCmec types of the SE INPer strains (Table 1; Fig. S6). The analysis demonstrated the presence of the community-acquired SCCmec type IV cassette in 13 out of 14 methicillin-resistant strains. The SE INPer strains S10 and S16 lack the SCCmec cassette, and no mecA gene was detected. Although the S21 strain has a mecA gene, we were unable to find other gene elements to show the presence of a mec cassette. Moreover, strains S07, S09, S13 strains contain an additional SCCmec type VIII cassette in tandem with the SCCmec IV cassette. The S07 strain carries a contig of 36 535 bps of the SCCmec type IV and VIII, suggesting the probable structure of the recombined cassette (Fig. S7).

The genomic analysis revealed that some SE strains studied here included genes for resistance to fluoroquinolones, macrolides, sulfonamides, aminoglycosides, tetracycline, and other antibiotics not used as the first choice in clinical therapy (Table 2). The results given in Table 2 corroborate that in most of the cases, the probable gene responsible for the resistance is present in the genomes. Besides, non-synonymous mutations in the antibiotic target proteins GyrA and RpoB were identified in some SE strains resistant to quinolones and rifampicin. Despite other SE strains lack these mutations, they still were resistant to these antibiotics. Then, other mutations in the antibiotic target protein or other genetic mechanisms not yet known would be responsible for these resistances.

Mobile genetic elements

The genomic variability observed in the SE strains suggests active processes of recombination and gene exchange. To study this concern, we first look for prophages and CRISPR-Cas related systems in the genomes. Prophage footprints were found in 10 out of 17 genomes of the SE strains. The most significant prophage hits detected by the PHAST program, were for genomic regions spanning about 28 to 95 kb that include an attachment site, a signature of lysogenic phages (Arndt et al., 2019) (Table S6). In this analysis, prophages were found integrated into the genomes of some SE strains, such as CNPH82 found in the SE strains S14, S17, and S18 (Daniel, Bonnen & Fischetti, 2007). Some other prophages such as StB20 and SpBeta were present in S03 and S16 strains respectively and, the prophages IPLA5 in strain S07 and IPLA7 in S12 and S15 strains (Gutiérrez et al., 2012). In the remaining SE strains, prophages sequences were not detected.

SE strains have also acquired defense mechanisms against phage infection. The search for CRISPR-Cas immune systems identified nine out of the 17 SE INPer strains carrying a CRISPR-Cas Type III system. The system III is composed by the Cas1 and Cas2, responsible for spacer processing and insertion; the ribonuclease Cas6, and the cascade proteins Csm1 to Csm6, involved in the processing of the target transcript (Table S7). The CRISPR-Cas Type III system has been already reported to confer immunity to phages as well as to conjugative plasmids in SE (Marraffini, 2015; Marraffini & Sontheimer, 2008). Although the enzymatic organization of the CRISPR-Cas systems is remarkably conserved in SE, there are variations in the array of repeats and spacers in the CRISPR loci. Three distinct types of identical repeated units of 30 or 36 nucleotides, associated with specific sequence spacers have been described (Marraffini, 2015). Three spacers that correspond to the CRISPR loci found in the strains S02, S05, and S24, match precisely with a sequence in the Staphylococcus phage PH15 genome for the first two strains and the Staphylococcus phage 6ec genome for the last (Aswani et al., 2014; Daniel, Bonnen & Fischetti, 2007).

SE strains harbor many ISs, belonging to different families (Table S8). The presence of IS256 has been found within pathogenic SE, associated with biofilm formation and virulence (Kozitskaya et al., 2004; Murugesan et al., 2018). Among the SE strains of our collection, the isolates having the IS256 also contain the ica operon, confirming previous observations. Exceptionally, only the strain S21 has the ica genes but lacks IS256.

Recombination

The above results suggest that high frequency of Horizontal Gene Transfer (HGT) and recombination might promote diversification of local SE populations. To evaluate this, we measured the ratio of recombination to mutation (r/m), using the ClonalFrameML program (Didelot & Wilson, 2015) in the 17 INPer SE genomes and 16 combined sets of SE genomes downloaded from GenBank (Fig. 2; Table S3). A median average r/m rate of about 6.9 was calculated when the 17 SE were tested, suggesting that nucleotide substitutions by recombination are more frequent than random point mutations (Vos & Didelot, 2009). Every COG class shows r/m values equal or higher than the estimate for the complete set of 17 SE genomes. Indeed, the r/m values on virulence or antibiotic resistance gene class result similar to the other COGs involved in housekeeping functions.

To determine whether or not the recombination estimates were affected by the sample composition of SE strains, we design several control tests, with distinct groups of genomes. First, we discarded the most divergent S10 strain of the SE collection and ran the ClonalFrameML test only with the 16 SE most related SE strains of the collection. As shown in Fig. 4, the r/m rate for this set reduced to a median of four, and this value is not significantly different respect to the r/m of SE of the complete genome of SE strains obtained from the GenBank (z-score = 1.4, P-value 0.135) (Fig. 4, boxplot 25). Second, we computed the r/m rate in eight different sets randomly selected among 260 complete and draft genomes of SE strains from the GenBank (Fig. 4, 26–33; Table S3). Some sets (Ctr2 and 8) display the lowest r/m values, whereas the rest control sets have r/m upper than two up to four. These results indicate that the strains sample composition influence the recombination estimates.

The RaxML nucleotide phylogenetic tree used as a reference to estimate recombination looks similar to the core protein phylogeny presented in Fig. 2; SE strains within clades maintain cohesive relationships (Fig. 5). However, multiple recombination events were detected in the ancestral nodes leading to the SE strains. The most prominent branch (red dot line in Fig. 5) divided the SE strains into two large clades, one including the clade D and the other constituted by clade B and C. Within different divergent lineages recombination introduces much more nucleotide variants than mutation as presented before (r/m). These results indicate that local hospital settings SE strains may contain enough genomic diversity despite their close relationship with the main clonal ST complexes of worldwide distribution.