Unsuspected clonal spread of Methicillin-resistant Staphylococcus aureus causing bloodstream infections in hospitalized adults detected using whole genome sequencing

Patient cohort selection

We identified all patients diagnosed with a MRSA BSI between July 2018 and June 2019, admitted to either of two hospitals of the University of Pennsylvania hospital system. The Hospital of the University of Pennsylvania (HUP) is a 625-bed academic tertiary and quaternary care medical center in West Philadelphia with approximately 32,000 patient admissions, 633,000 outpatient visits, and 40,000 Emergency Department visits annually. The Penn Presbyterian Medical Center (PMC) is a 324-bed urban community hospital in West Philadelphia with 12,000 admissions, 130,000 outpatient visits, and 26,000 Emergency Department visits annually. A single case of MRSA BSI was defined as any instance of a MRSA isolate collected from blood of any patient at HUP or PMC during the study period. Each subject was only included once. The study was approved by the University of Pennsylvania IRB and given a waiver of consent, as the study was retrospective and no data or samples were collected specifically for research purposes.

Isolate selection and DNA sequencing

Isolates were obtained from a biobank that collected all clinical MRSA isolates cultured from routine diagnostic studies in the HUP Clinical Microbiology Laboratory during the study period. They included all sequential MRSA BSI isolates from the study population. Isolates were screened for phenotypic resistance to twenty antibiotics using the Vitek 2 automated system, and susceptibility and resistance were assigned in accordance with CLSI protocols [17]. A 1μL loopful of each frozen isolate was streaked onto blood agar and incubated overnight at 37°C. A single isolated colony was selected from each plate and grown again on blood agar overnight at 37°C. A 10μL loopful of each isolate was then frozen in a bead beating tube and underwent whole genome sequencing (WGS) at the Penn/Children’s Hospital of Philadelphia Microbiome Center. Sequencing libraries were prepared using the Illumina Nextera library preparation kit, and sequenced on the Illumina MiSeq. Sequences were made publicly available through the NCBI Sequence Read Archive (Bioproject PRJNA751847).

Bioinformatic pipelines

Paired-end 150 bp FASTQ files were passed through the “Bactopia” workflow to assess data quality, assemble contigs, call MLST, SCCmec type, antibiotic resistance and virulence genes [18]. To compare SNP-based core-genome multiple-sequence alignments, the total number of assembled contigs or subsets grouped by clonal complex (CC) were passed through two reference-based core SNP alignment tools, Snippy (v.4.6.0) [19] and Parsnp (v.1.5.6) [20], or concatenated and aligned using tools within the PIRATE workflow [21]. ClonalFrameML [22] and maskrc-svg (v0.5)(https://github.com/kwongj/maskrc-svg) were used to identify and mask possible recombinant regions with in the core-genome alignment. Pairwise SNP distances for each alignment were calculated using snp-dists [23]. Maximum likelihood trees were created using IQ-Tree (v2.1.2) [24] and visualized using the ggTree package [25]. Bootstrap values of the most recent common ancestral node were calculated for 1000 repetitions. We used N315 strain genome (GCF_000009645.1) for total isolate alignments and CC5-specific alignments; for the CC8 specific alignment we used NCTC 8325 as reference (GCF_000013425.1). Phylogenetic similarity across alignment strategies was measured by calculating cophenetic correlation [26] between SNP distance matrices and estimated phylogeny tip distance and Robinson-Foulds distances [27] between trees generated with the different alignments as well as randomly generated trees using ape version 5.5 [28]. The pangenome-mapping of the 104 isolates was used because the tree had the strongest correlation between the estimated phylogeny, high tree structure similarity to other alignment tools, and conservation of SNP distances between alignment strategies (Table 1; Fig 1A-B).

• Download figure
• Open in new tab

Figure 1. Frequencies and distribution of single nucleotide polymorphism distances between isolates varies by alignment tool.

Whole-genome sequence alignments were created from MRSA alignments using either all 104 isolate sequences (blue) or isolates categorized only as clonal complex (CC) 8 and 5 (red), and using three alignment strategies: Pangenome-pipeline (A,B), Assembly-pipeline (C,D), and Pseudoread-pipeline (E,F). Pairwise single nucleotide polymorphism (SNP) distances were plotted for each alignment strategy.

Epidemiological Investigation of Clustered Isolates

We analyzed core genomes with pairwise SNP differences less than or equal to two thresholds to detect patient-patient transmission: 35 SNPs and 15 SNPs. The threshold of 35 SNPs was based on the approximate cutoff found in within-patient versus between-patient BSI MRSA lineages in a hospital setting [6,29]. The more stringent threshold of 15 SNPs for core genome comparisons was proposed as a conservative cutoff for recent inter-patient MRSA transmission [13]. Demographic data, comorbidities, Pitt Bacteremia Score, source of BSI, and in-patient mortality were collected for each study subject by abstraction from the electronic medical record (EMR). A BSI was considered to be healthcare-associated (HA) if the index blood culture was drawn >48 hours after hospital admission; healthcare-associated, community-onset (HACO) if the index culture was obtained <48 hours after admission or in the community setting, and if the subject had one or more previous healthcare risk factors (hospitalization, surgery, hemodialysis, or stay in a nursing home or another residential medical facility in the previous year; or presence of an indwelling intravascular catheter at time of culture); and community-associated (CA) if the index culture was obtained <48 hours after admission or in the community setting and the subject lacked these healthcare risk factors. Demographic and clinical characteristics were tabulated for all subjects. The EMR was examined for evidence of overlap or sequential hospital/unit stays among cluster subjects, defined as cases in which subjects’ MRSA isolates from >1 subject clustered at the 35 SNP threshold, and visualized using the vistime package (https://github.com/shosaco/vistime). Dates of admission and discharge for all hospital stays at any of the four networked hospitals in the one year before the first collected BSI isolate in a cluster and one year after the last collected BSI isolate in the cluster were recorded for each cluster subject. In addition to HUP and PMC, any relevant hospitalizations were recorded at Pennsylvania Hospital (PH) and at a single, University of Pennsylvania long-term acute care hospital in Philadelphia during the year prior to their index BSI culture. PH is a 481-bed urban community hospital located in the Society Hill district of Philadelphia with over 27,000 hospital admissions, over 24,000 Emergency Department visits, and 201,000 outpatient visits annually.

Logistic regression modelling was used to assess the predictive power of pairwise SNP distances and likelihood of patient temporal hospitalization overlaps. Goodness of fit was assessed using a receiver operating characteristic (ROC) curve and measuring the area under the curve. All analyses were conducted in R studio (v1.4.1106)[30] run with R version 4.0.4, and final figures labelled in InkScape (v0.92.5) [31]. Analysis code is available at https://github.com/Read-Lab-Confederation/MRSA_bloodstream_clusters.

Distribution of patient demographics and isolate characteristics

We screened all patients diagnosed with a MRSA BSI at three academic hospitals between July 2018 and June 2019, identifying 106 qualifying subjects (see Materials and Methods). Of the BSI source sites that could be identified from EMR, skin site infections made up 19% and central venous catheter infections made up 14% (Table 2). Among included subjects, 16% died while in the hospital.

From each individual, single MRSA samples were isolated and sequenced, of which 105 had sufficient coverage for further analysis and 104 isolates were determined to be Staphylococcus aureus. One isolate was identified by WGS as Staphylococcus argenteus and was excluded from further analysis. Among the 104 genomes, 55 were assigned to CC8, 49 of which were USA300 strains; 40 were assigned to CC5; and the remaining nine isolates were in CC30, CC72, and CC78.

Assessment of sequence alignment strategies

We generated multiple alignments of the MRSA isolates using three approaches to determine their effect on pairwise SNP distances between isolates (Fig. 1). The three approaches were (1) randomly fragmenting assembled genomes to create “pseudoreads” and mapping these to a reference genome using Snippy (“pseudoread-pipeline”); (2) mapping assembled genomes to the reference using Parsnp (“assembly-pipeline”); and (3) aligning concatenated core genes that were identified using the PIRATE tool (“pangenome-pipeline”). Alignments generated with all 104 isolates included had lower pairwise distances compared to alignments generated only with isolates of the same CC. This was because the alignments with larger numbers of more diverse strains always covered a smaller amount of the chromosome (Table 2). Pairwise SNP distances produced by the pangenome-mapping method were consistent between the CC groups and whole species alignments (Fig 1A-B), whereas the pairwise SNP distances produced by pseudoread-mapping and assembly-mapping methods were greater when smaller groups of more closely-related strains were the input (Fig. 1 C-F).

Suspected transmission clusters

Using alignments containing all 104 isolates for each pipeline, 29 isolates grouped into nine clusters (C1-C9) where isolates had less than a 35-SNP threshold difference from at least one other subject isolate (Table 3). The psuedoread-pipeline clustered 29 isolates, the assembly-pipeline clustered 21, and the pangenome-pipeline clustered 19. Five clusters contained isolates from CC5, three clusters were CC8, and one cluster was CC30. Three clusters also contained sub-cluster groups of isolates which fell within the SNP threshold for only some of the tools. One cluster (C7) was only detected by assembly-mapping. The median cluster size was three isolates (range 2 - 6). Blood culture collection date differences between clustering isolates was wide, with the longest period being 265 days (C1), and the shortest being 12 days (C6). Median SNP differences were diverse across clusters, and smaller differences did not correlate with shorter collection date differences.

Phylogenetic analysis of patient isolates

BSI isolates, and those isolates falling into possible transmission clusters, spanned across significantly divergent clades of CCs (SH-aRLT and uBS values >70) (Fig 2A). Candidate transmission clusters formed distinct sub lineages within the clades. The largest cluster, C5, diverged significantly from other CC8 isolates and were identified as part of the CC8c lineages, a strain of USA500 as defined by the phylogenetic subtyping scheme by Bowers et al. [32]. Cluster and non-cluster isolates had a similar distribution for infection setting, with most BSI isolates categorized as HACO (68%). At a 15-SNP threshold, only isolates in clusters C1, C2, C5(a,c), and C8 remained clustered. All isolates were susceptible to vancomycin and daptomycin, but isolates in both the CC5 and CC8 clades contained multi-antibiotic resistant phenotypes to beta-lactams (as expected for MRSA) and quinolones. Together, these data suggested that multiple lineages of MRSA associated with BSI had the potential to transmit between patients, and those clusters had a higher likelihood of harboring multi-class resistance phenotypes.

• Download figure
• Open in new tab

Figure 2. Suspected transmission clusters fall into distinct clonal groups.

Maximum likelihood trees were generated from the PIRATE alignment of 104 isolates and visualized using ggtree. (A) Tree indicating clades containing individual clonal complexes (CCs). (B) Subtrees from the complete maximum likelihood trees for the two most abundant CCs. Nodes with bootstrap values >= 70 are marked in red. Heat maps show strain type, SCCmec element type, and resistance phenotype for indicated antibiotics per sequence, infection setting (Healthcare-associated (HA), Community associated (CA) and Healthcare-associated community-onset (HACO)), admission hospital, and transmission cluster at a threshold of 35 SNPs or 15 SNPs.

Genomic similarity predicts overlapping hospital stay in transmission clusters

To compare common healthcare exposures among clustered isolates, we examined the hospitalization history for every subject within a cluster at four networked hospitals in the University of Pennsylvania system one year before the first index BSI isolate and one year after the last patient index isolate. Six had subjects with overlapping hospital stays (Fig. 3), of which three had median SNP distances between 1-16 core SNPs with corresponding hospital unit overlaps (Table 3; Fig. 3). Cluster C5c had a median SNP difference of seven (range 6-10 SNPs across tools) with no common hospital overlap. In comparison, cluster C4 also had no subjects with overlapping hospital admissions prior to their index BSI cultures, but the median SNP distance range of 20 - 26 SNPs across tools.

• Download figure
• Open in new tab

Figure 3. Hospitalization history among patients in genomic BSI clusters.

Hospitalization history at 4 study hospitals (A, B, C, and D) up to 365 days before the date of the earliest MRSA bloodstream isolate culture in each cluster (relative Day 0) and up to 365 days after the latest MRSA bloodstream isolate in the cluster. Note that bloodstream infections were only included at hospitals A and B. Rows represent the hospitalization history of each patient associated with a sequenced cluster isolate. Colored rectangles and circular marks represent individual hospitalization durations (rectangles) or one-day admissions (circles); the color indicates Hospital, A, B, C, or D. Black outlined boxes represent areas where two or more patients overlapped in the same hospital at the same time. Red stars indicate the date of collection of the sequenced BSI isolate for each patient. Yellow triangles indicate a hospitalization where two or more patients overlapped in the same hospital unit.

The SNP threshold is a critical component for identifying clusters for further investigation. We therefore performed a logistic regression to determine if there was an association between likely hospital exposure and SNP threshold (Fig. 4). The log odds of clustered patient pairs overlapping in the same hospital decreases by 0.065 with every increase of one SNP (p=0.05). Our logistic model showed that with 13 SNPs the likelihood that any two patients in a cluster overlapped in a hospital was 50 percent (Fig. 3A). Pairwise SNP distances that were greater than 13 SNPs trended toward patients with no overlapping hospital stays in the year prior to their MRSA bacteremia infection. The area under the curve of the receiver operating characteristic analysis indicates that we could classify true positive exposures (i.e., known prior overlapping hospitalizations) 66% of the time from the SNP threshold (Fig. 3B).

• Download figure
• Open in new tab

Figure 4. Higher SNP distances trend toward ruling out hospital overlaps between clustering patients.

(A) Logistic regression model indicating the relationship between patient pairs overlapping in the same hospital at the same time (prior to the diagnosis of an index MRSA bloodstream infection) and the pairwise SNP distance. Points indicate the true result for each pair as overlapping (1.0) or not overlapping (0). The color of the points indicates whether hospital overlap patient pairs also overlapped (black) or did not overlap (gray) in the same hospital unit. Gray ribbon indicates the 95% confidence interval. (B) Receiving operating characteristic (ROC) curve of the logistic model in A. Area under the curve (AUC) = 0.662.