Evolution and Global Transmission of a Multidrug-Resistant, Community-Associated Methicillin-Resistant Staphylococcus aureus Lineage from the Indian Subcontinent

The Bengal Bay clone (ST772) is a community-associated and multidrug-resistant Staphylococcus aureus lineage first isolated from Bangladesh and India in 2004. In this study, we showed that the Bengal Bay clone emerged from a virulent progenitor circulating on the Indian subcontinent. Its subsequent global transmission was associated with travel or family contact in the region. ST772 progressively acquired specific resistance elements at limited cost to its fitness and continues to be exported globally, resulting in small-scale community and health care outbreaks. The Bengal Bay clone therefore combines the virulence potential and epidemiology of community-associated clones with the multidrug resistance of health care-associated S. aureus lineages. This study demonstrates the importance of whole-genome sequencing for the surveillance of highly antibiotic-resistant pathogens, which may emerge in the community setting of regions with poor antibiotic stewardship and rapidly spread into hospitals and communities across the world.

M ethicillin-resistant Staphylococcus aureus (MRSA) is a major human pathogen with a propensity to develop antibiotic resistance, complicating treatment and allowing persistence in environments where there is antibiotic selection pressure. While multidrug resistance has traditionally been the domain of health care-associated strains, the emergence of strains in the community setting that are also resistant to multiple antibiotics poses a significant challenge to infection control and public health (1). Given the heavy burden and costs associated with MRSA infections (2,3), there is an urgent need to elucidate the patterns and drivers behind the emergence of drug-resistant community-associated MRSA lineages.
Over the last few years, several population genomic studies have started to unravel the evolutionary history of community-associated S. aureus lineages emerging in specific regions of the world. The prototype of these clones is the diverse USA300 lineage (sequence type 8 [ST8]), forming distinct genetic lineages in North America and South America (4,5), including distinct clades in Europe and Africa (6). The East Asia clone (ST59) has diverged into two distinct lineages, with evidence of establishment in Taiwan and North America (7), while the ST80 lineage originated in North Africa but went through a notable population expansion to become the dominant communityassociated lineage in North Africa, the Middle East, and Europe (8). On the Australian continent, ST93 emerged in indigenous communities of Western Australia and the Northern Territory and spread to the eastern seaboard and sporadically overseas, while also forming a clade associated with Pacific Islander and Maori populations in New Zealand (9).
This diversity of regional evolutionary histories is reflected in the various factors that have been suggested to contribute to the emergence and establishment of community-associated clones. For instance, acquisition of Panton-Valentine leucocidin (PVL), a mutation in the capsule gene cap5D, and acquisition of the SCCmec-IV and ACME/COMER elements played a defining role in the regional evolution of USA300 (6). While acquisition of a typical SCCmec-IV element was also associated with the emergence of the eastern seaboard clade of ST93, host population factors such as household poverty, crowding, and a high burden of skin infections have been postulated to be associated with both the emergence in indigenous populations of Australia and the establishment of a Pacific Islander-associated clade in New Zealand (9). In contrast, the population expansion of the ST80 lineage was linked to the acquisition of a SCCmec-IV element and fusidic acid resistance, and enrichment of resistance determinants in the Taiwan clade of ST59 suggests a strong contribution for its emergence and persistence in East Asia. While there is some evidence from these studies that resistance acquisition can be a driving force behind the regional emergence of community-associated MRSA clones, there is a lack of data on strains from critical regions that are considered hot spots for the emergence of multidrug-resistant pathogens, such as the Indian subcontinent.
In 2004, a novel S. aureus clone designated ST772 was isolated from hospitals in Bangladesh (10) and from a community setting study in India (11), where it continued to be reported in community-and health care-associated environments (see data posted at https://doi.org/10.6084/m9.figshare.8061887.v3) (12). Similar to other S. aureus strains, ST772 primarily causes skin and soft tissue infections, but more severe manifestations, such as bacteremia and necrotizing pneumonia, have been observed. Its potential for infiltration into nosocomial environments (13)(14)(15)(16) and resistance to multiple classes of commonly used antibiotics (including aminoglycosides, ␤-lactams, fluoroquinolones, macrolides, and trimethoprim) (16)(17)(18) have made ST772 an alarming public health concern on the Indian subcontinent and elsewhere. Over the last decade, the clone has been isolated from community and hospital environments in Asia, Australasia, Africa, the Middle East, and Europe (see data posted at https://doi.org/10 .6084/m9.figshare.8061887.v3). As a consequence of its discovery, distribution, and epidemiology, the lineage has been informally dubbed the Bengal Bay clone (19). Despite clinical and epidemiological hints for a recent and widespread dissemination of ST772, a unified perspective on the global evolutionary history and emergence of the clone is lacking.
In this study, we analyzed whole-genome sequences from a globally representative collection of 340 ST772 strains to elucidate the key events associated with the emergence and global spread of a multidrug-resistant community-associated MRSA clone. Our analysis suggests that the clone originated on the Indian subcontinent in the 1960s and rapidly expanded through the 1990s and early 2000s. We found that international travel and family connections to the region (India, Bangladesh, Nepal, and Pakistan) were closely linked with the global spread of the lineage. Genome integration of a multidrug resistance plasmid appeared to be a driver in the emergence of a dominant clade (ST772-A) in the early 1990s.

RESULTS
We generated whole-genome sequence data for 354 S. aureus ST772 isolates collected across Australasia, South Asia, Hong Kong, the Middle East, and Europe between 2004 and 2013 (see data posted at https://doi.org/10.6084/m9.figshare .8061887.v3). Fourteen isolates were excluded after initial quality control due to contamination, and the remainder mapped with 165ϫ average coverage against the PacBio reference genome DAR4145 (18) from Mumbai. From a core genome of 2,545,215 bp, there were 7,063 single nucleotide polymorphisms (SNPs). Phylogenetic analysis using core genome SNPs revealed little geographic structure within the lineage (Fig. 1a). Eleven ST772 methicillin-susceptible S. aureus (MSSA) and MRSA strains were basal to a single globally distributed clade (ST772-A; n ϭ 329) that harbored an integrated resistance plasmid (IRP) described for the reference genome DAR4145 (18) ( Fig. 1a and b). Population network analysis distinguished three distinct subgroups within ST772-A ( Fig. 1a and c): an early-branching subgroup harboring multiple subtypes of the staphylococcal chromosome cassette (SCCmec-V) (A1; n ϭ 81), a dominant subgroup (A2; n ϭ 153), and an emerging subgroup (A3; n ϭ 56); the last two exclusively harbor a short variant of SCCmec-V.
Emergence and global spread from the Indian subcontinent. Epidemiological and genomic characteristics of ST772 were consistent with an evolutionary origin from the Indian subcontinent. Sixty percent of isolates in this study were collected from patients with family or travel background in Bangladesh, India, Nepal, or Pakistan; the other isolates had unknown origins (19%) or were from other countries (21%) (Fig. 2a). We found more isolates from India and Bangladesh among the basal strains than clade ST772-A (Fisher's exact test, 5/11 versus 47/291, P ϭ 0.026), suggesting an association of the emergence of the clone with the Indian subcontinent. In particular, three isolates from India and Bangladesh were basal in the (outgroup-rooted) maximum likelihood (ML) phylogeny ( Fig. 1b and Fig. S1A), including two MSSA samples from the original isolations in 2004 (RG28 and NKD122). Isolates recovered from South Asia were genetically more diverse than isolates from Australasia and Europe, supporting an origin from the Indian subcontinent ( Fig. 2b and Fig. S1B).
Consistent with a methicillin-susceptible progenitor, a higher proportion of MSSA isolates was found among the basal isolates (Fisher's exact test, 4/11 versus 31/291, P ϭ 0.028), and MSSA isolates demonstrated a lower patristic distance to the root of the maximum likelihood phylogeny than MRSA (Fig. S1C). Although it appears that MSSA is proportionately more common in South Asia, it is also possible that the observed distribution may be related to nonstructured sampling (Fig. S1C).
Phylogenetic dating suggests an initial divergence of the ancestral ST772 population in 1962 (age of root node, 1961.94; confidence interval [CI], 1942[CI], .75 to 1977, with a core genome substitution rate of 1.18 ϫ 10 Ϫ6 substitutions/site/year after removing recombination ( Fig. 3 and Fig. S1D). This was followed by the emergence of the dominant clade ST772-A and its population subgroups in the early 1990s (ST772-A divergence, 1990.02; 95% CI, 197795% CI, .5 to 1995. The geographic pattern of dissemination is heterogeneous (Fig. 1a). There was no evidence of widespread endemic dissemination of the clone following intercontinental transmission, although localized health care-associated outbreak clusters occurred in neonatal intensive care units (NICUs) in Ireland (NICU-1 and NICU-2 [ Fig. 1a and Fig. S2]) (20) and have been reported from other countries in Europe (16) and South Asia (13)(14)(15). While some localized spread in the community was observed among our isolates, patients in local transmission clusters often had traveled to or had family in South Asia (19/27 clusters [ Fig. S2]).
Antibiotic resistance acquisition is associated with emergence and dissemination. We examined the distribution of virulence factors, antibiotic resistance determinants, and mutations in coding regions to identify the genomic drivers in the emergence and dissemination of ST772. Nearly all isolates (336/340) carried the Panton-Valentine leucocidin (PVL) genes lukS/F, most isolates (326/340) carried the associated enterotoxin A gene (sea), and all isolates carried scn. This indicates a nearly universal carriage, across all clades, of both a truncated hlb-converting prophage (the typically associated staphylokinase gene sak was present in only one isolate) and the PVL/sea prophage -IND772 (21). Among other virulence factors, the enterotoxin genes sec and sel, the gamma-hemolysin locus, egc cluster enterotoxins, and the enterotoxin homologue ORF CM14 were ubiquitous in ST772. We detected no statistically significant difference between core virulence factors present in the basal group and ST772-A (Fig. S3).
MRSA isolates predominantly harbored one of two subtypes of SCCmec-V: a short variant (5C2) or a composite cassette (5C2&5), which contains a type 5 ccr complex including ccrC1 (allele 8) between the mec gene complex and orfX (22) (Fig. S4). Integration of the Tn4001 transposon carrying aminoglycoside resistance gene aadA-aphD occurred across isolates with different SCCmec types (260/267) but not in MSSA (0/35). Of the 11 isolates in the basal group, 4 were MSSA (three lacked SCCmec and one had a remnant SCCmec-IV element), and the 7 MRSA isolates all carried the larger composite cassette SCCmec-V (5C2&5), with 3 of these 7 strains having a variant of SCCmec-V (5C2&5) (Fig. 1a).
The diversity of SCCmec types decreased as ST772-A diverged into subgroups ( Fig. 1a and c). ST772-A1 included MSSA (n ϭ 30) as well as SCCmec-V (5C2 [n ϭ 22] and 5C2&5 [n ϭ 18]) strains. Four isolates harbored a putative composite SCC element that included SCCmec-V (5C2), as well as pls and the kdp operon previously known from SCCmec-II. One isolate from the United Kingdom harbored a composite SCCmec-V (5C2&5) with a cadmium and zinc resistance locus (czrC), known from the European livestock-associated CC398-MRSA (23). Another six isolates yielded irregular and/or composite SCC elements. In contrast, the dominant subgroups ST772-A2 and -A3 exclusively carried the short SCCmec-V (5C2) element. Eleven of these isolates (including all isolates in NICU-2) lacked ccrC and two isolates carried additional recombinase genes (ccrA/B2 and ccrA2). Considering the distribution of SCCmec across the lineage, SCCmec integration appears to have occurred on multiple occasions in the basal strains and ST772-A1, with subsequent modifications of SCCmec-V (5C2 and 5C2&2) (Fig. 1a). In contrast, integration appears to have occurred only on a single occasion at the divergence of ST772-A2 and -A3, followed by the occasional complete or partial loss of the element in these subgroups.
Predicted resistance to erythromycin was uniquely found in ST772-A and not in any of the basal strains (Fisher's exact test, 289/291 versus 0/11, P Ͻ 0.001 [ Fig. 5a]), characterized by the acquisition of an integrated multidrug resistance plasmid (IRP) (Fig. 5b), including the macrolide resistance locus msrA/mphC, as well as determinants

Australasia
Europe South Asia showed that the element is rare and predominantly plasmid associated across ST8 genomes (6/274), with one chromosomal integration in the ST772 reference genome and the SCCmec integration in the ST80 reference genome. Three basal strains were not multidrug resistant and included two isolates from the original collections in India (RG28) and Bangladesh (NKD122) ( Fig. 1a and 4a). These two strains lacked the trimethoprim determinant dfrG and the fluoroquinolone mutations in grlA or gyrA, including only a penicillin resistance determinant blaZ on a Tn554-like transposon. However, seven of the strains more closely related to ST772-A did harbor mobile elements and mutations conferring trimethoprim (dfrG) and quinolone (grlA and blue, n ϭ 244; South Asia: purple, n ϭ 52). Error bars indicate 95% confidence intervals using nonparametric bootstrapping. Isolates from the Arabian Peninsula (n ϭ 2) and Hong Kong (n ϭ 6) were excluded from the diversity analysis due to the small number of samples from these regions. Evolution and Global Transmission of ST772-MRSA ® gyrA mutations) resistance. Interestingly, we observed a shift from the quinolone resistance grlA S80F mutation in basal strains and ST772-A1 to the grlA S80Y mutation in ST772-A2 and -A3 (Fig. 4a). Thus, the phylogenetic distribution of the key resistance elements suggests acquisition of the IRP by a PVL-positive MSSA strain in the early 1990s (ST772-A1 divergence, Canonical mutations and phenotypic comparison of basal strains and ST772-A. We found three other mutations of interest that were present exclusively in ST772-A strains. The first mutation caused a nonsynonymous change in fbpA (L55P), encoding a fibrinogen-binding protein that mediates surface adhesion in S. aureus (26). The second comprised a nonsynonymous change (L67V) in the plc gene, encoding a phospholipase associated with survival in human blood cells and abscess environments in USA300 (27). The third encoded a nonsynonymous mutation (S273G) in tet (38), an efflux pump that promotes resistance to tetracyclines as well as survival in abscess environments and skin colonization (28). The functional implication of genes harboring these canonical mutations might suggest a modification of the clone's ability to colonize and cause skin and soft tissue infections.
In light of these canonical SNPs, we selected 5 basal strains and 10 strains from ST772-A to perform a preliminary screen for potential phenotypic differences that may contribute to the success of ST772-A. We assessed in vitro growth, biofilm formation, cellular toxicity, and lipase activity (Fig. S6). We found no statistically significant differences between the basal strains and ST772-A in these phenotypic assays, apart from significantly lower lipase activity among ST772-A strains (Welch's two-sided t test, t ϭ 3.4441, degrees of freedom ϭ 6.0004, and P ϭ 0.0137), which may be related to the canonical nonsynonymous mutation in plc. However, it is increased rather than decreased lipase activity that has been associated with viability of S. aureus USA300 in human blood and neutrophils (27). We found no difference in the median growth rate of ST772-A compared to the basal strains (Mann-Whitney, P ϭ 0.8537), although there were two ST772-A strains that grew more slowly, suggesting the possibility of some strain-to-strain variability.

DISCUSSION
In this study, we used whole-genome sequencing in combination with epidemiological and phenotypic data to investigate the drivers behind the emergence and spread of a multidrug-resistant community-associated MRSA lineage from the Indian subcontinent. Our data suggest that the Bengal Bay clone has acquired the multidrug resistance phenotype of traditional health care-associated MRSA but retains the epidemiological characteristics of community-associated MRSA.
Emergence of a basal population of ST772 appears to have occurred on the Indian subcontinent in the early 1960s and included strains from the original isolations of ST772 in Bangladesh and India in 2004. While genomic surveillance studies from India are rare, two recent studies have detected ST772-MSSA and -MRSA in Nepal (29) and ST772-MRSA in Pakistan (30), but it is unclear whether the lineage was endemic in these countries prior to its emergence in India. Deeper genomic surveillance of ST772-MSSA and -MRSA in the region will be necessary to understand the local epidemiology and evolutionary history of the clone on the Indian subcontinent.
Establishment and expansion of a single dominant clade (ST772-A) occurred in the early 1990s and was associated with the acquisition of an integrated multidrug resistance mobile genome element (MGE). The element is similar to a previously described extrachromosomal plasmid of USA300 (24) and a partially integrated element in the SCCmec of an ST80 reference genome (25). While the element was found only once in the ST80 lineage (8) and occurs predominantly on plasmids in closed ST8 (USA300) genomes, its distribution and contribution to the emergence of resistance in the ST8 lineage have so far not been addressed (4,6). In contrast, the ubiquitous occurrence and retention of the element in ST772-A suggest that it was instrumental in the emergence of the dominant clade of the Bengal Bay clone.
Furthermore, we observed a replacement of the long composite SCCmec-V (5C2&5) element with the shorter SCCmec-V (5C2) and fixation of the quinolone resistance mutation from grlA S80F to the grlA S80Y as ST772-A2 and -A3 became the dominant population subgroups in the 1990s. In light of earlier studies demonstrating a fitness advantage in having a smaller SCCmec element (31)(32)(33), the fixation of the shorter SCCmec-V (5C2) may be a contributing factor to the success of ST772.
We observed a lack of significant differences in growth between basal strains and ST772-A. This may suggest that acquisition of drug resistance was not accompanied by a major fitness cost to ST772-A and raises the possibility that members of this clade will survive in environments where antibiotics are heavily used, such as hospitals or in communities with poor antibiotic stewardship, but may also be at little disadvantage in environments where there is less antibiotic use. Notably, competitive-fitness experiments involving the community-associated MRSA lineages USA300 and ST80 revealed that the biological cost of resistance to methicillin, fusidic acid, and fluoroquinolones is reversed in the presence of trace amount of antibiotics (34). In regions such as the Indian subcontinent, where community use of antibiotics is not well regulated, it is plausible that lineages such as ST772-A could thrive in such environments. While our results are consistent with those of Gustave et al. (34), it should be noted that our phenotypic assays assessed only a small number of isolates and the growth assays were conducted in vitro and under nutrient-rich conditions that are unlikely to capture anything but very high fitness costs. Future studies may build on this preliminary work and consider a more in-depth analysis of a larger number of isolates.
Given the available epidemiological data, phylogeographic heterogeneity, and the clone's limited success to establish itself in regions outside its range of endemicity in South Asia (Fig. 1a), there appears to be ongoing exportation of ST772 from the Indian subcontinent, associated with travel and family background in the region. This is supported by reports of MRSA importation in travelers, including direct observations of ST772 importation by returnees from India (35). Our data suggest nonendemic spread within households and the community, including short-term outbreaks at two NICUs in Ireland. This pattern of limited endemic transmission is consistent with reports of small transmission clusters in hospitals and households during a comprehensive surveillance study of ST772 in Norway (16).
There have been small numbers of ST772 isolates identified from other countries and regions of the world not included in this study, including Austria, Finland and Slovenia, Taiwan, Japan, China, and Nigeria (36)(37)(38)(39)(40)(41)(42)(43). These studies have generally been broader descriptions of circulating MRSA genotypes, of which only a small minority have been found to be ST772. These reports would therefore support similar patterns of either sporadic importation or short-lived chains of local transmission of ST772 as observed for isolates in our study.
Overall, the pattern of spread mirrors those of other community-associated MRSA lineages, such as USA300 (5, 44), ST80-MRSA (8), and ST59 (7), where clones emerge within a particular geographic region, are exported elsewhere, but rarely become established and endemic outside their place of origin. In contrast, health careassociated MRSA clones such as CC22-MRSA-IV (EMRSA-15) (45) and ST239-MRSA-III (46,47) demonstrate much stronger patterns of phylogeographic structure, consistent with importation into a country followed by local dissemination through the health care system. While there are indications for resistance acquisition driving regional community-asscoiated lineages, such as the Taiwan clade of ST59 (7), we found strong indications in our study that the acquisition of mutations and mobile elements associated with multidrug resistance was the dominant driver behind the emergence of the Bengal Bay clone on the Indian subcontinent and its subsequent intercontinental transmission. Moreover, we observed an unusual, near complete lack of phylogeographic structure in the population, compared to other previously investigated community-associated clones, providing evidence for ongoing circulation and exportation from the Indian subcontinent followed by limited endemic transmission.
Our data trace the evolution of the Bengal Bay clone on the Indian subcontinent, where it emerged in the 1960s and diverged into a single dominant clade in the 1990s. Its rapid emergence may have been driven by the dissemination of mobile genetic elements, particularly those that confer drug resistance, such as the acquisition of a multidrug resistance integrated plasmid and variants of SCCmec. Patient epidemiology and phylogenetic heterogeneity suggest a pattern of ongoing exportation from the Indian subcontinent and limited endemic transmission after importation. The Bengal Bay clone therefore appears to combine the epidemiological characteristics of community-associated MRSA lineages with an unusually resistant genotype traditionally seen in health care-associated MRSA.
Considering the widespread use of antibiotics and associated poor antibiotic regulation, limited public health infrastructure, and high population density in parts of South Asia, the emergence and global dissemination of multidrug-resistant bacterial clones (both Gram positive and Gram negative) are alarming and perhaps not surprising. Global initiatives and funding to monitor the occurrence of emerging clones and resistance mechanisms, and support for initiatives in antimicrobial stewardship at community, health care, and agricultural levels, are urgently needed.

MATERIALS AND METHODS
Isolates. Isolates were obtained from Australia (21 isolates The collection was supplemented with six previously published genome sequences from India (21,22,48). Notable samples include the initial isolates from Bangladesh and India (10, 11), two hospital-associated (NICU) clusters from Ireland (20), and isolates from a single health care worker at a veterinary clinic sampled over two consecutive weeks (VET) (49). Geographic regions were designated Australasia (Australia and New Zealand), East Asia (Hong Kong), South Asia (India, Bangladesh), Arabian Peninsula (Saudi Arabia and United Arab Emirates), and Europe (Denmark, England, Germany, Ireland, Italy, the Netherlands, Norway, and Scotland).
Clinical data and epidemiology. Anonymized patient data were obtained for the date of collection, clinical symptoms, geographic location, epidemiological connections based on family or travel history, and acquisition in nosocomial or community environments, where available. Clinical symptoms were summarized as skin and soft tissue infections (abscesses, boils, ulcers, exudates, pus, and ear and eye infections), urogenital (vaginal swabs and urine), bloodstream (bacteremia), or respiratory infections (pneumonia and lung abscesses), and colonization (swabs from ear, nose, throat, perineum, or the environment) (Fig. S5). Literature and sample maps (see data posted at https://doi.org/10.6084/m9 .figshare.8061887.v3) were constructed with geonet, a wrapper for geographic projections with Leaflet in R (https://github.com/esteinig/geonet).
Where available, acquisition in community or health care environments was recorded in accordance with guidelines from the CDC. Community-associated MRSA is therein classified as an infection in a person who has none of the following established risk factors for MRSA infection: isolation of MRSA more than 48 h after hospital admission; history of hospitalization, surgery, dialysis, or residence in a longterm-care facility within 1 year of the MRSA culture date; the presence of an indwelling catheter or a percutaneous device at the time of culture; or previous isolation of MRSA (50,51) (Fig. S5).
A valid epidemiological link to South Asia was declared if either travel or family background could be reliably traced to Bangladesh, India, Nepal, or Pakistan. If both categories (travel and family) were unknown or where data were available for only one category and did not show a link to the region, we conservatively declared the link as unknown or absent, respectively. The longitudinal collection (n ϭ 39) from a staff member at a veterinary hospital in England was treated as a single patient sample.
MLST and SCC typing. In silico multilocus sequence typing (MLST) was conducted using mlst (https://github.com/tseemann/mlst) on the assembled genomes with the S. aureus database from PubMLST (https://pubmlst.org/saureus/). Three single-locus variants (SLVs) of ST772 were detected and retained for the analysis, describing ST1573, ST3362, and ST3857. Sequences of experimentally verified sets of probes for SCC-related and other S. aureus-specific markers (62,63) were BLAST searched against SPAdes assemblies (in silico microarray typing), allowing prediction of the presence or absence of these markers and detailed typing of SCC elements. We assigned MRSA to four isolates that failed precise SCC Variant calling. Samples passing quality control (n ϭ 340) were aligned to the PacBio reference genome DAR4145 from Mumbai, and variants were called with the pipeline Snippy (available at https://github.com/tseemann/snippy), which wraps BWA MEM, SAMtools, SnpEff (65), and Freebayes (66). Core SNPs were defined as being present in all samples (ignoring insertions and deletions, n ϭ 7,063) and were extracted with snippy-core at default settings. We assigned canonical SNPs for ST772-A as those present exclusively in all isolates of ST772-A but not in the basal strains. Annotations of variants were based on the reference genome DAR4145.
Phylogenetics and recombination. A maximum likelihood (ML) tree under the general time reversible (GTR) model of nucleotide substitution with among-site rate heterogeneity across 4 categories (GTR ϩ ⌫), ascertainment bias correction (Lewis), and 100 bootstrap (BS) replicates was generated based on 7,063 variant sites (core genome SNPs) in RaxML-NG 0.5.0 (available at https://github.com/amkozlov/ raxml-ng), which implements the core functionality of RAxML (67). The tree with the highest likelihood out of 10 replicates was midpoint rooted and visualized with interactive Tree of Life (ITOL) (Fig. 1a and  2a and Fig. S2 and S7B) (68). In all phylogenies (Fig. 1a, Fig. 2a, and Fig. 3 and Fig. S2, S7, and S8), samples from the veterinary staff member were collapsed for clarity.
A confirmation alignment (n ϭ 351) was computed as described above for resolving the pattern of divergence in the basal strains of ST772. The alignment included the CC1 strain MW2 as an outgroup, as well as another known SLV of CC1, sequence type 573 (n ϭ 10). The resulting subset of core SNPs (n ϭ 25,701) was used to construct an ML phylogeny with RaxML-NG (GTR ϩ ⌫) and 100 bootstrap replicates (Fig. S1A). We also confirmed the general topology of our main phylogeny as described above using the core genome alignment of 2,545,215 nucleotides generated by Snippy, masking sites if they contained missing (Ϫ) or uncertain (N) characters across ST772. Phylogenies were visualized using ITOL, ape (69), phytools (70), ggtree (71), or plotTree (https://github.com/katholt/plotTree). Patristic distances to the root of the phylogeny (Fig. S1C) were computed in the adephylo (72) function distRoot.
Dating analysis. We removed sites with missing data (Ϫ) from the core alignment generated by Snippy and then ran Gubbins (73) to detect homologous recombination events, using a maximum of five iterations and the GTR ϩ ⌫ model in RaxML. Following the removal of SNPs related to recombination there were 6,907 SNPs in the alignment. We used Least Squares Dating (LSD) v0.3 (74) to obtain a time-scaled phylogenetic tree. This method fits a strict molecular clock to the data using a least-squares approach. Importantly, LSD does not explicitly model rate variation among lineages, and it does not directly account for phylogenetic uncertainty. However, its accuracy is similar to that obtained using more sophisticated Bayesian approaches (75), with the advantage of being computationally less demanding.
LSD typically requires a phylogenetic tree with branch lengths in substitutions per site and calibrating information for internal nodes or for the tips of the tree. We used the phylogenetic tree inferred using maximum likelihood in PhyML (76) (after removing recombination with Gubbins, as described above) using the GTR ϩ ⌫ substitution model with 4 categories for the ⌫ distribution. We used a combination of nearest-neighbor interchange and subtree-prune-regraft to search tree space. Because PhyML uses a stochastic algorithm, we repeated the analyses 10 times and selected that with the highest phylogenetic likelihood. To calibrate the molecular clock in LSD, we used the collection dates of the samples (i.e., heterochronous data). The position of the root can be specified a priori, using an outgroup or by optimizing over all branches. We chose the latter approach. To obtain uncertainty around node ages and evolutionary rates, we used the parametric bootstrap approach with 100 replicates implemented in LSD.
An important aspect of analyzing heterochronous data is that the reliability of estimates of evolutionary rates and timescales is contingent on whether the data have temporal structure. In particular, a sufficient amount of genetic change should have accumulated over the sampling time. We investigated the temporal structure of the data by conducting a regression of the root-to-tip distances of the maximum likelihood tree as a function of sampling time (77) and performing a date randomization test (78). Under the regression method, the slope of the line is a crude estimate of the evolutionary rate, and the extent to which the points deviate from the regression line determines the degree of clocklike behavior, typically measured using R (79). The date randomization test consists of randomizing the sampling times of the sequences and reestimating the rate each time. The randomizations correspond to the distribution of rate estimates under no temporal structure. Consequently, the data have strong temporal structure if the rate estimate using the correct sampling times is not within the range of those obtained from the randomizations (80). We conducted 100 randomizations, which suggested strong temporal structure for our data (Fig. S1D). We also verified that the data did not display phylogenetictemporal clustering, a pattern which sometimes misleads the date randomization test (81).
Nucleotide diversity. Pairwise nucleotide diversity and SNP distance distributions for each region with Ͼ10 samples (Australasia, Europe, South Asia) were calculated as outlined by Stucki et al. (82). Pairwise SNP distances were computed using the SNP alignment from Snippy (n ϭ 7,063) and the dist.dna function from ape, with raw counts and deletion of missing sites in a pairwise fashion. An estimate of average pairwise nucleotide diversity per site () within each geographic region was calculated from the SNP alignments using raw counts divided by the alignment length. Confidence intervals for each region were estimated using 1,000 bootstrap replicates across nucleotide sites in the original alignment via the sample function (with replacement) and the 2.5% to 97.5% quantile range (Fig. 2b).
Population structure. We used the network analysis and visualization tool NetView (83, 84) (available at http://github.com/esteinig/netview) to delineate population subgroups in ST772. Pairwise Hamming distances were computed from the core SNP alignment derived from Snippy. The distance matrix was used to construct mutual k-nearest-neighbor networks from k ϭ 1 to k ϭ 100. We ran three commonly used community detection algorithms as implemented in igraph to limit the parameter choice to an appropriate range for detecting fine-scale population structure: fast-greedy modularity optimization (85), Infomap (86), and Walktrap (87). We thereby accounted for differences in the mode of operation and resolution of algorithms. Plotting the number of detected communities against k, we were able to select a parameter value at which the results from the community detection were approximately congruent (Fig. S7A).
Since we were interested in the large-scale population structure of ST772, we selected a k value of 40 and used the low-resolution fast-greedy modularity optimization to delineate final population subgroups. Community assignments were mapped back to the ML phylogeny of ST772 (Fig. 1a). All subgroups agreed with the phylogenetic tree structure and were supported by Ն99% bootstrap values (Fig. S7B). One exception was isolate HW_M2760, located within ST772-A2 by phylogenetic analysis but assigned to ST772-A3 by network analysis (Fig. S7A and B). This appeared to be an artifact of the algorithm, as its location and connectivity in the network representation matched its phylogenetic position within ST772-A2. The network and communities were visualized using the Fruchtermann-Reingold algorithm (Fig. 1c), excluding samples from the veterinary staff member in Fig. 1c (Fig. S7A).
Local transmission clusters. We obtained approximate transmission clusters by employing a network approach supplemented with the ML topology and patient data, including date of collection, location of collection, and patient family links and travel or family links to South Asia. We used pairwise SNP distances to define a threshold of 4 SNPs, corresponding to the maximum possible SNP distance obtained within 1 year under a core genome substitution rate of 1.61 ϫ 10 Ϫ6 nucleotide substitutions/ site/year. We then constructed the adjacency matrix for a graph, in which isolates were connected by an undirected edge, if they had a distance of less than or equal to 4 SNPs. All other isolates were removed from the graph, and we mapped the resulting connected components to the ML phylogeny, showing that in each case the clusters were also reconstructed in the phylogeny, where isolates diverged from a recent common ancestor (gray highlights in Fig. 2a and Fig. S2). We then traced the identity of the connected components in the patient metadata and added this information to each cluster. NICU clusters were reconstructed under these conditions ( Fig. 2a and Fig. S2).
Antimicrobial resistance, virulence factors, and pan-genome analysis. Mykrobe Predictor was employed for antibiotic susceptibility prediction and detection of associated resistance determinants and mutations. Mykrobe Predictor has a demonstrated sensitivity and specificity Ͼ99% for calling phenotypic resistance and is comparable to gold-standard phenotyping in S. aureus (64). Predicted phenotypes were therefore taken as a strong indication for actual resistance phenotypes in ST772. Genotype predictions also reflect multidrug resistance profiles (aminoglycosides, ␤-lactams, fluoroquinolones, macrolideslincosamides-streptogramin B [MLS], and trimethoprim) reported for this clone in the literature (16)(17)(18)(19)(20)88). As most resistance-associated MGEs in the complete reference genome DAR4145 are mosaic-like and flanked by repetitive elements (18), we used specific diagnostic genes present as complete single copies in the reference annotation of DAR4145 (18) to define the presence of the IRP (msrA) and Tn4001 (aacA-aphD). Mykrobe Predictor simultaneously called the grlA mutations S80F and S80Y for quinolone resistance phenotypes. However, in all cases one of the variants was covered at an extremely low median k-mer depth (Ͻ20), and we consequently assigned the variant with a higher median k-mer depth at grlA (see data posted at https://doi.org/10.6084/m9.figshare.8061887.v3).
ARIBA (89) with default settings and the core Virulence Factor database were used to detect the complement of virulence factors in ST772. We corroborated and extended our results with detailed in silico microarray typing, including the presence of the egc gene cluster or S. aureus-specific virulence factors such as the enterotoxin homologue ORF CM14. Differences in detection of relevant virulence factors between the in silico typing and ARIBA included, among others, lukS/F-PVL (337 versus 336), sea carried on the -IND772 prophage (336 versus 326), sec (333 versus 328), and sak (1 versus 2). Since in silico microarray typing was based on assembled genomes and may therefore be prone to assembly errors, we used results from the read-based typing with ARIBA to assess statistical significance of virulence factors present in basal strains and ST772-A (Fig. S3).
We searched for the three resistance regions which aligned to the 11819-07 and the 11809-03 plasmid (DAR4145 reference genome; R1, 1,456,024 to 1,459,959 bp; R2, 1,443,096 to 1,448,589 bp; and R3, 1,449,679 to 1,453,291 bp) in all completed S. aureus genomes (including plasmids) in RefSeq (NCBI) and the NCTC3000 project (http://www.sanger.ac.uk/resources/downloads/bacteria/nctc/) using nctctools (https://github.com/esteinig/nctc-tools) and nucleotide BLAST with a minimum of 90% coverage and identity (n ϭ 273). Since the IRP is mosaic-like and composed of several mobile regions, we retained query results only if all three of the regions were detected. We then traced the integration sites in the accessions, determining whether integrations occurred in the chromosome or plasmids. Multilocus sequence types were assigned using mlst (https://github.com/tseemann/mlst).
Growth curves. S. aureus strains were grown overnight in 5 ml of tryptic soy broth (TSB; Fluka) with shaking (180 rpm) at 37°C. Overnight cultures were diluted 1:1,000 in fresh TSB, and 200 l was added to a 96-well plate (Costar) in triplicate. Growth was measured 37°C with shaking (300 rpm) using a FLUOROstar fluorimeter (BMG Labtech) at an absorbance wavelength of 600 nm. Growth curves represent the means of triplicate results.
Cell culture conditions. The monocyte/macrophage THP-1 cell line was maintained in suspension in 30 ml of RPMI 1640 medium supplemented with 10% heat-inactivated fetal bovine serum (FBS), 1 M L-glutamine, 200 U/ml of penicillin, and 0.1 mg/ml of streptomycin at 37°C in a humidified incubator with 5% CO 2 . Cells were harvested by centrifugation at 700 ϫ g for 10 min at room temperature and resuspended to a final density of 1 ϫ 10 6 to 1.2 ϫ 10 6 cells/ml in tissue-grade phosphate-buffered saline (PBS), typically yielding Ͼ95% viable cells as determined by easyCyte flow cytometry (Millipore).
Human erythrocytes were harvested from 10 ml of human blood following treatment in sodium heparin tubes (BD). Whole blood was centrifuged at 500 ϫ g for 10 min at 4°C. Supernatant (plasma) was aspirated and cells were washed twice in 0.9% NaCl and centrifuged at 700 ϫ g for 10 min. The cell pellet was gently resuspended in 0.9% NaCl and diluted to 1% (vol/vol).
Cytotoxicity assay. To monitor S. aureus toxicity, S. aureus strains were grown overnight in TSB, diluted 1:1,000 in 5 ml of fresh TSB, and grown for 18 h at 37°C with shaking (180 rpm). Bacterial supernatants were prepared by centrifugation of 1 ml of bacterial culture at 20,000 ϫ g for 10 min. For assessing toxicity to THP-1 cells, 20 l of cells was incubated with 20 l of bacterial supernatant and incubated for 12 min at 37°C. Both neat and 30% diluted supernatant (in TSB) were used, as certain S. aureus strains were considerably more toxic than others. Cell death was quantified using easyCyte flow cytometry using the Guava viability stain according to the manufacturer's instructions. Experiments were done in triplicate. For assessing hemolysis, 150 l of 1% (vol/vol) erythrocytes were incubated with 50 l of either neat or 30% supernatant in a 96-well plate for 30 min at 37°C. Plates were centrifuged for 5 min at 300 ϫ g, 75 l of supernatant was transferred to a new plate, and absorbance was measured at 404 nm using a FLUOROstar fluorimeter (BMG Labtech). Normalized fluorescence was achieved using the equation (A t -A 0 )/(A m /A 0 ), where A t is the hemolysis absorbance value of a strain, A 0 is the minimum absorbance value (negative control of 0.9% NaCl), and A m is the maximum absorbance value (positive control of 1% Triton X-100).
Lipase assay. Bacterial supernatants used in the above-described cytotoxicity assays were also used to assess lipase activity, using the protocol published by Cadieux et al. (92), with modifications. Briefly, 8 mM para-nitrophenyl butyrate (pNPB), the short-chain substrate, or para-nitrophenyl palmitate (pNPP), the long-chain substrate (Sigma), was mixed with a buffer (50 mM Tris-HCl [pH 8.0], 1 mg/ml of gum arabic, and 0.005% Triton X-100) in a 1:9 ratio to create assay mixtures. A standard curve using these assay mixtures and para-nitrophenyl (pNP) (Sigma) was created, and 200 l of each dilution was pipetted into 1 well of a 96-well plate (Costar). A total of 180 l of each assay mixture was pipetted into the remaining wells of a 96-well plate, and 20 l of the harvested bacterial supernatant was mixed into the wells. The plate was placed in a FLUOstar Omega microplate reader (BMG Labtech) at 37°C, and a reading at 410 nm was taken every 5 min for 1 h. The absorbance readings were converted to micromolar concentration of pNP released per minute using the standard curve.
Biofilm formation. Semiquantitative measurements of biofilm formation on 96-well, round-bottom polystyrene plates (Costar) were made based on the classical, crystal violet method of Ziebuhr et al. (93). Eighteen-hour bacterial cultures grown in TSB were diluted 1:40 into 100 l of TSB containing 0.5% glucose. Perimeter wells of the 96-well plate were filled with sterile H 2 O, and plates were placed in a separate plastic container inside a 37°C incubator and grown for 24 h under static conditions. Following 24 h of growth, plates were washed five times in PBS, dried, and stained with 150 l of 1% crystal violet for 30 min at room temperature. Following five washes of PBS, wells were resuspended in 200 l of 7% acetic acid, and optical density at 595 nm was measured using a FLUOROstar fluorimeter (BMG Labtech). To control for day-to-day variability, a control strain (E-MRSA15) was included on each plate in triplicate, and absorbance values were normalized against this. Experiments were done using six technical repeats from 2 different experiments.
Statistical analysis. All statistical analyses were carried out in R or python and considered significant at a P value of Ͻ0.05, except for comparisons of proportions across the multiple virulence and resistance elements, which we considered be statistically significant at a P value of Ͻ0.01. Veterinary samples (n ϭ 39) were restricted to one isolate (one patient, Staff_E1A) for statistical comparison of region of isolation, proportion of resistance, virulence, and MSSA between basal strains and ST772-A (n ϭ 302; see above [ Fig. 5a and Fig. S3]). Differences in pairwise SNP distance and nucleotide diversity between all regions were assessed using nonparametric Kruskal-Wallis test and post hoc Dunn's test for multiple comparisons with Bonferroni correction, as distributions were assumed to be not normal (n ϭ 340 [ Fig. 2b and Fig. S1B]). Phenotypic differences were assessed for normality with Shapiro-Wilk tests. We consequently used either Welch's two-sided t test or the nonparametric two-sided Wilcoxon rank sum test (Fig. S6).
Data availability. Supplemental data tables and interactive files include a summary of publication and isolates used in this study, raw data from various analyses, and measurement of phenotype experiments, as well as interactive maps and gene comparisons. Supplemental data tables and files, with the exception of supplemental figures, are hosted on Figshare (https://doi.org/10.6084/m9.figshare .8061887.v3). Core analyses, including parameter settings, cluster resource configurations, and versioned software distributions, are reproducible through the bengal-bay-0.1 workflow, which can be found along with other scripts and data files at our GitHub repository (https://github.com/esteinig/ST772). The workflow is implemented in Snakemake (94) and runs in virtual environments that include software distributed in the Bioconda (95) channel. Analyses were conducted on the Cheetah cluster at the Menzies School of Health Research, Darwin, Australia.