ABSTRACT
Whole genome sequencing is increasingly applied to healthcare-associated vancomycin resistant Enterococcus faecium (VREfm) outbreak investigations using single bacterial colonies from individual patients. However, within-patient diversity could complicate transmission resolution. We sequenced 14 colonies from VREfm positive rectal swabs over a 1-month period on a haematology unit. In total, 224 isolates from 11 patients were sequenced. Carriage of 2-3 sequence types was detected in 27% of patients. Presence of antimicrobial resistance genes and plasmids was variable within isolates from the same patient and sequence type. We identified two dominant sequence types (ST80 and ST1424), with two putative transmission clusters of two patients within ST80, and a single cluster of six patients within ST1424. We found transmission resolution was impaired using fewer than 14 colonies, and therefore there is a need to capture and consider VREfm within-patient diversity to ensure accurate resolution of transmission networks.
INTRODUCTION
Enterococcus faecium is a leading nosocomial pathogen causing opportunistic infections mostly in immunocompromised hosts. Antimicrobial resistance is a key concern, particularly against front-line anti Gram-positive agents amoxicillin and vancomycin (1). Vancomycin resistant E. faecium (VREfm) infections cause significant issues due to increased length of stay, cost an estimated USD200 (GBP ∼167) per case per day, and confer mortality of 23-47% (2–7). In 2020, vancomycin resistance of 45.6% was reported among all E. faecium bloodstream isolates in Scotland, among the highest rates in Europe (8).
VREfm does not cause symptomatic gastrointestinal disease but asymptomatic intestinal carriage is common. In healthcare settings, this can increase the risk of shedding into the environment and transfer to other patients or staff, hence challenging efforts to limit the incidence of nosocomial infections in high-risk settings (9). Whole genome sequencing is increasingly applied to investigate transmission networks and identify control measures (10,11). Most whole genome sequencing (WGS) based analyses of bacterial outbreaks, however, rely on analysing single colony picks from clinical samples making the assumption that this represents the entire infecting or colonizing population within individual patients (12). However, it is increasingly recognised that within-patient diversity of bacterial populations can be significant from a clinical and infection prevention control point of view, and can influence transmission network resolution (13–19). Several studies have identified that individual patients can carry multiple strains of E. faecium concurrently, but few have applied this to transmission resolution (20–24).
In this study we aimed to identify within-patient diversity of VREfm from rectal screening swabs and determine how this impacts transmission inference in a 1-month snapshot on a haematology unit. We designed a sampling strategy to reliably detect within-patient diversity, and supplemented short read and long read sequencing to generate high-quality reference genomes to identify genomic variants in the isolate collection.
RESULTS
Epidemiological Context
This study was performed over a single month in 2017 on an 18-bed unit for patients undergoing management of haematological malignancies. The unit is split into two wards (A and B). Due to the high risk of invasive infection in this patient group, VREfm rectal screening was performed on all new febrile admissions and any inpatients with febrile episodes to inform patient placement and antimicrobial administration. Given the short study timeframe there was significant overlap between patient stays with some patients moving between the two study wards or to other wards in the hospital (Figure 1). There was a policy in place for patients to be cohorted or placed in side rooms based on known colonisation or infection with high risk pathogens (including Candida auris, Clostridioides difficile, methicillin resistant Staphylococcus aureus, multiresistant Gram negative organisms, norovirus, respiratory viruses, varicella zoster virus, and VREfm). However, as not all rooms had ensuite facilities, risk of VREfm transmission remained. At the time of the study, pathogen surveillance systems in the hospital had not detected any suspected VREfm outbreak within the study population.
Patient timeline. Each row denotes the location of a patient during admission, blocks denote hospital stay, circles denote VREfm cultures, stars denote bloodstream isolates, dotted lines indicate the start and end of prospective collection of screening isolates for this study. This study was undertaken mainly within Wards A and B, although patients were moved to different wards within the hospital during their stay and were often admitted through the assessment unit.
Results of VREfm Screening
In total, 45 rectal swabs from 27 patients were screened for VREfm. Of these, 18 samples from 13 patients were VREfm positive (Table 1). Additionally, three (23.1%) colonised patients developed VREfm bacteraemia 9, 24, or 46 days after being identified as VREfm carriers. We applied our sampling strategy to 16 rectal screens generating 224 isolates from 11 patients, and five blood cultures (five isolates) from two patients yielding a total of 229 isolates. Two rectal swabs and one blood culture were not available for further study. Most patients were female, and the median age was 66 years, a range of primary diagnoses were represented in the patient cohort (Table 1). Most colonised patients had received antibiotics in the preceding six months, although only 30% had received vancomycin (Table 1).
Characteristics of patients with rectal VREfm colonisation (n = 13)
Simultaneous carriage of multiple VREfm strains
In silico MLST typing using short reads from all 229 isolates showed four sequence types (ST) from the hospital-associated clade A1 (25) were present in the sample cohort, dominated by ST80 (n=130) and ST1424 (n=97) (Table 2). ST789 and ST1659 were identified in single isolates only. Multiple STs were detected in three (27%) samples. Sample VRED06 from patient P49 carried 10 (71.4%) ST80, three (21.4%) ST1424, and one (7.1%) ST789 isolate; sample VRED07 from P14 carried 10 (71.4%) ST1424 and four (28.6%) ST80 isolates; sample VRED11 from P50 carried 13 (92.9%) ST1424 and one (7.1%) ST1659 isolate. A further rectal swab sample from P49 collected two days after VRED06 contained only ST1424, and a blood culture collected nine days later also contained ST1424. P9 had three rectal swab samples collected over an 11 day period and had positive blood cultures one month later, all samples contained ST80 only.
MLST Sequence Types Detected
Generation of high quality reference genomes
As 99% of isolates grouped into only two STs we next sought to identify relationships within these two sequence types. To generate strain-specific reference genomes for further analysis one ST80 (VRED06-10) isolate and one ST1424 (VRED06-02) isolate were selected for long read sequencing. Both isolates were from P49 and were from the same swab sample collected on study day 15. Complete genomes were assembled and polished using the Trycycler pipeline generating a single circular chromosome and 5-7 plasmids for each isolate (Supplementary Table 1). BUSCO showed 99.5% complete core genes and indel analysis showed <4% coding sequences were significantly shorter than the reference match, within the expected range of truncated genes due to evolutionary processes (26). These results show the assemblies were complete and contained very few errors.
Genomic population structure of VREfm isolates suggests recent transmission events
The chromosomes of the two strain-specific genome assemblies were used as references for short read mapping within each sequence type. Within-patient diversity was low when isolates of the same MLST ST were compared, generally differing by zero SNPs and a maximum pairwise difference of 4 SNPs (Table 2). Where within-patient isolates were from different STs, SNP differences were significantly higher as expected (490-4,400 SNPs, Table 2). Similarly, detection of insertions, deletions, and plasmids were usually shared in isolates from the same patient. However, the presence of DEL3 (12 bp non-coding deletion) and DEL4 (11 bp deletion in a solute binding protein) were variable within 24 ST80 isolates from P20 with 0-2 differentiating SNPs (Figure 2). In isolates from P9 p1_VRED06-10 and p3_VRED06-10 were variably detected despite most genomes having no differentiating SNPs (Figure 2). Where multiple samples from the same patient were collected over time we found low (0-3 SNPs) accumulation of SNPs and no pattern in the prevalence of other genomic variants, reflecting the short study time period.
Phylogeny of ST80 isolates showing structured population with three patient specific clusters and two clusters indicating putative patient-patient transmission of VREfm. All ST80 isolates mapped to VRED06-10 chromosome and phylogeny built on SNP sites (n=96) after removal of putative transposable and recombinant regions.
The ST80 genomes formed a well-structured population with five clear clusters each separated by >10 SNPs (Figure 2). Within the clusters isolates were closely related, differing by 0-2 SNPs, and were mostly derived from individual patients. However, two clusters included isolates from two different patients (patients P7 and P33, and P2 and P9), suggesting two recent transmission events. All the reference plasmids were detected in the P7 and P33 isolates and there was a difference of two insertions in three isolates. As mentioned above there was significant variation in detection of p1_VRED06-10 and p3_VRED06-10 plasmids within P9 isolates, although in P2 isolates all plasmids were detected.
Mapping of the ST1424 isolates showed a much more homogeneous population than in ST80 (Figure 3). Of the 97 ST1424 isolates, 69 had no SNPs and the remaining 28 had 1-2 SNPs differentiating them from the rest of the collection. The SNPs that were detected did not lead to any clear clustering of isolates, except for the 14 isolates from P6 which all carried a SNP in a penicillin binding protein which differentiated them from the other ST1424 isolates. Two of the P6 isolates had further independent SNPs (one each) and another one isolate had lost p1_VRED06-02. Two isolates from P49 collected on day 17 shared a SNP in a penicillin binding protein as well as a 13 bp deletion in a non-coding region. The same variants were detected in a blood culture isolate collected 13 days later. Of note, the specific SNP in these three P49 isolates was in the same penicillin binding protein gene as the P6 isolates, but a different region. No insertions were detected in the ST1424 collection, and of the six deletions found five were only in isolates from P49. Apart from the loss of p6_VRED06-02 mentioned above, the only other plasmid not detected was p1_VRED06-02 in five isolates from three patients.
Phylogeny of ST1424 isolates showing homogeneous population suggestive of recent transmission outbreak. ST1424 isolates (n=98) mapped to VRED06-02 chromosome and phylogeny built on SNP sites (n=13) after removal of putative transposable and recombinant regions.
Analysis of multiple VREfm colonies supports transmission resolution
A transmission network was then constructed taking account of the phylogenetic placement of all 14 colony picks in each sample (Figure 4). The network supports transmission of ST80 between P2 and P9, and between P7 and P33, although P20 appears to be a singleton not linked to transmission within this patient cohort. Considering patient movements, P9 was admitted to Ward B prior to the start of the study period and screened positive for VREfm on day 13 of their admission, P2 screened positive on the same day although this was their first day of admission to Ward B (Figure 1). Both P2 and P9 were admitted through Ward A prior to being placed on Ward B, although no shared rooms or bed spaces were identified (Supplementary Figure 1). The directionality or location of this putative transmission event cannot be ascertained. P7 screened positive on admission to Ward B. P33 was admitted from an external Health Board to Ward B, screened negative on day 2 of their admission, then screened positive on day 10 of admission (6 days after P7). No shared rooms or bed spaces were identified for P7 and P33, although given the timings of positive screens it is likely that ST80 VREfm was transmitted from P7 to P33 on Ward B (Figure 1, Supplementary Figure 1). P20 screened negative towards the end of a 32 day stay on Ward B, was discharged for 13 days and then readmitted to Ward B and screened positive on admission, suggesting they may have become colonised outside of the hospital.
Phyloscanner transmission network. Each patient is represented by a node coloured by detection of the two outbreak STs. Edge thickness corresponds to fraction of Phyloscanner trees with given relationship, relationship fraction is printed alongside each edge, and edge colour based on type of relationship (orange, direct transmission; blue, transmission but direction unclear).
Patients carrying ST1424 all clustered together with P34 having strong links to all other patients, including likely direct transmission to P6 (Figure 4). P34 was admitted to Ward B from an external Health Board and screened positive with ST1424 on day 6 of admission, 6 days later P49 screened positive. P49 was on Ward B two days prior to P34 but screened negative 3 days before P34 was found to be colonised, suggesting P34 was the source of ST1424 transmission to P49. Soon after, P14 was found to be colonised with ST1424 and ST80 on day 2 of admission to Ward B but transmission from P49 is unlikely as ST80 isolates in both patients were unrelated and P49 had an additional ST789 strain. It is possible that P14 and P49 had superinfection of ST1424 on top of existing VREfm colonisation. Nine days later, P24 screened positive with ST1424 on Ward B on their seventh day of admission. P24 shared time on the ward with all three previously identified carriers, based on genetic similarity of all sequenced isolates P34 is the most parsimonious source of transmission (Figure 3). P6 and P50 were both admitted to Ward A on the same day late in the study period, P50 screened positive with ST1424 and ST1659 on day 2 and P6 screened positive with ST1424 only on day 6. From the epidemiology it may be assumed that P6 was infected by P50 however, sequencing suggests the two ST1424 populations may be derived from different hosts, with P6 isolates having all acquired a single SNP and P50 isolates having multiple different SNPs and all lacking the p6_VRED06-02 plasmid (Figures 3 and 4). P6 and P50 were the only patients on Ward A to acquire ST1424, the route of transmission between the two wards is unclear. None of the patients with ST1424 shared a room or used the same bed space previously used by an identified ST1424 positive carrier during their stay (Supplementary Figure 1).
Performing the analysis again with three, five, or ten genomes per screening sample did not change the transmission inference for ST80 isolates (although higher confidence was achieved with five or more genomes), but there was significant variability in the ST1424 cluster due to low diversity in the genomes (Supplementary Figure 2). Within ST1424 it was only possible to reliably link all patients and give high confidence that P34 was the main driver in the cluster with 14 isolates.
Plasmids were mostly ST-specific
Focusing on the plasmids identified in the reference isolates, VRED06-02 (ST1424) contained seven plasmids, and VRED06-10 (ST80) contained five plasmids. Plasmids in the two isolates were not particularly similar, suggesting there was limited sharing of plasmids between the two isolates within P49 at the time of sampling (Supplementary Table 2).
We then sought to identify carriage of similar plasmids in the entire isolate collection by mapping the short reads to each assembled plasmid (Supplementary Table 3). This showed that for the most part plasmids were ST-specific as there were few examples of ST1424 isolates carrying plasmids from the ST80 reference, and vice versa. However, all ST80 isolates from P7 and P33 carried p7_VRED06-02 from ST1424, and almost all isolates appeared to carry p4_VRED06-10. p4_VRED06-02 and p4_VRED06-10 had the closest pairwise identity (Supplementary Table 2), with mash distance of 0.05 and average nucleotide identity of 89% over 4.3 kb. As p4_VRED06-10 is slightly smaller, we believe the hits against the ST1424 isolates here are due to cross-mapping of reads from the related p4_VRED06-02. P7_VRED06-02 carriage in P7 and P33 does not look like a false positive as it is unrelated to others in the collection (Supplementary Table 2), but no close links to any ST1424 positive patients were identified for P7 and P33 (Figure 4), so the source of this plasmid is unclear. We also found the single ST1659 isolate from P50 had an identical plasmid profile to 12 ST1424 isolates in the same sample, and the single ST789 isolate from P49 had the same plasmid profile as 10 ST80 isolates in the same patient sample.
AMR gene load differs between closely related isolates
We next sought to determine the variability of AMR genes within the collection (Figure 5). In total 13 AMR genes were detected with three (23.1%) present in all isolates, five (38.5%) genes (ant(9)-Ia, dfrG, erm(A), and tet(M)) only in ST1424 isolates and two (15.4%) genes (ant(6)-Ia and tet(S)) found only in ST80 isolates. Tetracycline resistance gene tet(M) was found in 62.2% of ST1424 and ST1659 genomes and the aminoglycoside resistance gene aac(6’)-aph(2’’) was found in 69.9% of all genomes (Table 3 and Figure 5).
Detection of antimicrobial resistance genes. Panels represent different patients; resistance genes are plotted on the y-axis and isolates on the x-axis. Presence of a gene is represented by a filled square and coloured based on the MLST ST of the isolate.
Presence of AMR Genes
tet(M) was identified on the chromosome of VRED06-02, as part of Tn6944 first identified in C. difficile (Supplementary Figure 3A). Excision of Tn6944 is likely responsible for variable presence of tet(M).
aac(6’)-aph(2’’) was present on p1_VRED06-02 (ST80) and p1_VRED06-10 (ST1424). In p1_VRED06-02, two copies of aac(6’)-aph(2’’) were surrounded by insertion sequences IS256, IS1216, and IS3, potentially providing multiple mechanisms of excision. aac(6’)-aph(2’’) was not detected in any ST80 isolates that were p1_VRED06-02 negative, although only 60.2% (n=65) of isolates that carried this plasmid also carried aac(6’)-aph(2’’), highlighting that although the cassette was linked to p1_VRED06-02 it may have been excised in some cases. In ST1424 p1_VRED06-10, aac(6’)-aph(2’’) was surrounded by two copies of IS256 and carried genes including erm(B) with similarity to Tn6218 first identified in C. difficile, although the transposition machinery was missing (Supplementary Figure 3B) (27). In ST1424 aac(6’)-aph(2’’) was detected in 97.8% (n=90) isolates with p1_VRED06-10 but was also found in four isolates that lacked this plasmid (Table 3 and Figure 5). Short read assemblies were too fragmented to fully resolve the location of the gene in these four cases, but three of the assemblies did show aac(6’)-aph(2’’) in association with an IS3 gene suggesting it had mobilised to another genetic element within the cell.
The tetracycline resistance gene tet(L) was identified in a single ST1659 isolate, the gene was co-located with tet(M) on a 30 kb contig that was similar to Tn6248 from E. faecium over ∼19 kb (Supplementary Figure 3C).
DISCUSSION
By taking account of within-patient diversity in VREfm carriage populations we identified unknown links between patients that could supplement efforts to control transmission within hospitals. We also show that diversity exists not just at the level of SNPs – AMR gene presence/absence, indels, and plasmid presence all vary within and between patients. The relative importance of each of these layers of diversity for understanding bacterial transmission and evolution are key for future WGS based analyses.
Identified transmission clusters
We identified an ST1424 transmission cluster involving six patients and two ST80 clusters involving two patients (Figure 4). The ST80 population appeared relatively stable otherwise, with three clusters of isolates belonging to individual patients. However, in P14 and P49 where ST1424 and ST80 isolates were both present in carriage samples we did not detect transmission of ST80 suggesting these individuals were infected by an ST1424 on top of existing VREfm carriage, or a wide transmission bottleneck and acquisition of both lineages at the same time from an undetermined source.
Power calculation determined 14 colonies to be the minimum needed to sequence per sample to ensure reliable detection of a variant making up 20% of the within-patient population. This successfully identified minor population variants in 27% of patients. We also performed transmission inference with three, five, and ten isolates per rectal sample and found that the strength and accuracy of transmission inference increased as the number of isolates per patient increased (Figure 4 and Supplementary Figure 2). Sequencing more than 14 colonies would further improve the detection of minor variants but would also increase the costs of sequencing and turnaround time. Our findings show that sequencing single colonies without considering within patient diversity in suspected VREfm outbreaks will miss some putative transmission events, but the implementation of sequencing multiple colonies into routine use will require careful thought and validation as there was low within-patient diversity in most cases in our collection. Alternatively, strain-resolved metagenomics directly on clinical samples or sequencing sweeps of selective culture growth may be more feasible in future (28–30). Further work is required to determine the optimum sampling strategy to support infection prevention and control investigations in healthcare settings.
We identified three (27%) patients that carried two or three different strains at the same time (Table 2). This in line with recent studies that identified up to half of patients carry 2-4 different E. faecium strains and within-patient diversity varies over time, with identical AMR elements passing between unrelated strains within patients (20,24,31,32). These findings highlight the challenges of tracking and controlling highly variable E. faecium populations within healthcare settings, particularly when relying on single colony picks and only short read sequencing which cannot fully resolve most plasmids.
We found very little SNP-based diversification over time in our study, with ST1424 isolates often being identical to isolates collected weeks earlier, although we did identify a defining SNP in the P6 isolates collected late in the study (Figure 3). Estimates of diversification rates in E. faecium from national single colony sampling suggest 7 mutations per year (33), although other studies of longitudinal within-patient diversification have estimated higher rates of 12.6 – 128 mutations per year (20,22,23). The low SNP diversity identified in our one-month collection of carriage isolates is in keeping with these mutation rates, although it is perhaps unexpected that in the four P9 bloodstream isolates only one SNP was detected in one isolate compared to rectal samples collected 35-37 days previously (Figure 2). The low SNP level diversity leads to difficulty in reliably inferring specific transmission events within ST1424, as multiple potential transmission events are plausible, as reflected in the transmission network (Figure 4).
Variability in AMR and Mobile Elements
We identified some variability of AMR gene carriage within patients, probably related to loss of mobile elements (Table 3 and Figure 5). However, the impact on phenotype is unclear – all isolates carried aac(6’)-Ii and aph(3’)-III which together confer high level resistance to the clinically relevant aminoglycosides amikacin and gentamicin, so the loss of aac(6’)-aph(2’’) may be more efficient for the cell without an overt phenotypic change in minimum inhibitory concentration. tet(M) was present in 62% of ST1424 genomes but no other tetracycline resistance genes were found in these isolates, so this variability could have impacts on clinical reporting depending on which colony is selected for antibiotic susceptibility testing. However, tetracyclines are not generally used for direct treatment of enterococcal human infections so the clinical impact of this variability will be limited. Similar variable presence of the vancomycin resistance element within patients has been described elsewhere and could lead to inadvertent use of vancomycin even when the patient harbours a resistant population (20,21,34,35).
We identified a linear plasmid (p2_VRED06-10) carrying vanA in our ST80 reference genome. Boumasmoud et al (36) recently described the linear plasmid pELF_USZ in VREfm from Switzerland that carried an operon conferring the ability to utilise the human gut mucin N-acetyl-galactosamine; however, p2_VRED06-02 did not carry this operon. The advent of long read sequencing has allowed the increasing identification of linear plasmids, and their impact on niche adaptation remains to be fully elucidated (36–40).
We recognise some limitations. Around 60% of E. faecium carriers in hospitals can be linked to nosocomial transmission either from other patients or reservoirs in the hospital environment (31,41–44). Our study did not include environmental samples, although patients were mostly located in individual rooms. Bathroom facilities were shared posing a significant environmental reservoir for VREfm. Also, we relied on direct plating to solid VREfm screening agar for inclusion in our study. Previous studies have shown a sensitivity of 58-96% for this approach, rising to 97-100% with a pre-enrichment step (45–47). In our study, P49 carried ST1424, ST80, and ST789 on their first positive screen having been negative for VREfm 11 days previously, suggesting either infection with a mixed population from another host not detected by our sampling strategy or an initial false negative VREfm screen. Our focus on VREfm is another potential limitation as vancomycin sensitive E. faecium may contribute to transmission clusters causing significant disease and could gain/lose the resistance mechanism during transmission as described by Raven et al (33) – although isolates in the identified clusters were collected over years rather than weeks in our study. Despite these limitations, we believe our study are applicable to current practice as surveillance systems are mainly focussed on significant antimicrobial resistance phenotypes such as vancomycin resistance. Routine, prospective, genomic surveillance over longer periods than our study may be more suited to a resistance agnostic approach to understanding nosocomial spread of pathogenic bacteria (31).
Our findings suggest that sequence-based surveillance of VREfm screening samples would be useful in identifying outbreaks proactively. A proactive approach would avoid large outbreaks of infection, reduce costs of ward closure and the clinical impact of invasive disease (48–51). In our setting an outbreak of VREfm was suspected 3 weeks after the study collection period when P9 and P49 developed bloodstream infection concurrently but this was many weeks after VREfm transmission had likely occurred (Figure 1 and Figure 4). Due to the retrospective nature of our study, we were not able to use the findings from sequencing to directly influence patient care.
CONCLUSION
We have shown that sequencing 14 VREfm carriage isolates per patient improved resolution of possible transmission pathways. This has implications for outbreak investigations which may miss possible transmission events if only single colonies are investigated. We also show the utility of long read sequencing to supplement routine short read sequencing to generate high-quality, closed reference genomes and allow identification of non-SNP variation such as indels, AMR genes, and plasmid loss/gain events which can be used for fine-level discrimination within suspected bacterial outbreaks.
MATERIAL AND METHODS
Bacterial Isolates
Screening for VREfm intestinal carriage was performed by taking rectal swabs as part of admission sampling pre-chemotherapy/transplant, and on all inpatients on the haematology unit developing febrile neutropenia (defined as neutrophils <0.9×109/l or <1.0×109/l and falling after chemotherapy, plus body temperature ≥38°C). Swabs were plated to CHROMID® VRE agar (bioMérieux, Marcy-l’Étoile, France) and incubated at 37°C in air. Plates were examined for growth at 24 and 48 hours. Growth of purple colonies were presumed VREfm and a single colony was subcultured to Columbia blood agar plates (Oxoid, Basingstoke, UK) with a 5 μg vancomycin disc. Species identification was confirmed by MALDI-TOF on a Microflex instrument (Bruker, Billerica, USA), and the first positive isolate per patient had antimicrobial susceptibility determined on the VITEK2 system AST-P607 card (bioMérieux). Once confirmed as VREfm, all purple colonies were removed from the initial screening plate within 24 hours and stored at - 80°C in a Microbank cryovial (Pro-Lab Diagnostics, Birkenhead, UK). Any VREfm isolated from clinical samples (blood, wounds, or urine samples) within 60 days of the patient being identified as a VREfm carrier were also stored. Electronic records were reviewed for demographic and inpatient stay metadata. Patient movements were visualised with HAIviz v0.3 (https://haiviz.beatsonlab.com/). This work used excess diagnostic material and was approved by the NHS Scotland BioRepository Network (Ref TR000126) and University of St Andrews Research Ethics Committee (Ref MD12651).
DNA Extraction
Two cryobank beads were used to inoculate a CHROMID® VRE agar plate and incubated for 20-24 hours at 37°C in air. Fourteen purple colonies were selected at random, each inoculated into 5 ml of brain heart infusion broth (Oxoid) and incubated at 37°C for 20-24 hours. Bacterial cells were pelleted by centrifugation (3000 g, 10 min), resuspended in 400 μl Buffer P1 (Qiagen, Hilden, Germany), and 200 μl used for DNA extraction.
Cells were lysed with 20 μl lysozyme (100 mg/ml) (Sigma-Aldrich, Dorset, England) at 37°C for 60 min, followed by 30 μl Proteinase K (Qiagen) at 56°C for 60 min, and RNA was digested with 4 μl RNase A (Qiagen) at 37°C for 45 min. DNA was extracted from lysates on the QiaSymphony instrument using the DNA Mini kit (Qiagen) and protocol Tissue_HC_200_V7_DSP. DNA concentration was determined from 2 μl extract with the Qubit BR Kit (ThermoFisher Scientific, Basingstoke, UK) on a Qubit 3.0 fluorometer, extract purity was assessed with a Nanodrop 1000 (ThermoFisher Scientific).
Isolates were processed in batches of 24, with every 8th position being 5 ml of uninoculated BHI as a negative control to identify possible carry-over contamination. Negative controls were processed as above and assessed for purity using the Qubit HS Kit (ThermoFisher Scientific) in triplicate. Negative controls were accepted if DNA concentration was <0.5 ng/μl in all replicates.
Genome Sequencing
Short read libraries were prepared using the Nextera XT kit (Illumina, San Diego, USA) and sequenced with a MiSeq instrument (Illumina) using 300 bp paired-end reads on a 600-cycle v3 reagent kit.
Extracts from two isolates representative of the main clades identified in the collection were selected for long read sequencing. Long read libraries were generated with the LSK109 Ligation Sequencing Kit (Oxford Nanopore Technologies, Oxford, UK) and sequenced for 8h using an R9.4 flow cell on a GridION sequencer (Oxford Nanopore Technologies). High accuracy real time basecalling was performed in MinKNOW v19.12.6 (Oxford Nanopore Technologies).
Sequence data from this study has been deposited in the NCBI under BioProject accession number PRJNA877253 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA877253).
Sequence Assembly and Mapping
Long reads were filtered to remove <1000 bp fragments using Nanofilt v2.7.1 (52), then trimmed with Porechop v0.3.2 (https://github.com/rrwick/Porechop) to remove adapters and split chimeric reads. Genome assemblies were generated with Trycycler v0.3.3 (53). Trycycler makes multiple long read assemblies from read subsets and produces a consensus long read assembly which is then polished with long and short reads. Four assemblies were generated using Flye v2.8.1, Redbean v2.5, Raven v1.1.10, and Miniasm v0.1.3 (54–57), giving 12 assemblies in total. The final consensus assembly was polished with Medaka v0.11.5 (https://github.com/nanoporetech/medaka) and Pilon v1.23 (58), the latter for 2-3 cycles until no further changes were made. Assembly quality was assessed with assembly-stats v1.0.1 (https://github.com/sanger-pathogens/assembly-stats), Ideel (https://github.com/phiweger/ideel) and Busco v4.1.4 (59).
Hybrid assemblies were annotated with Prokka v1.14.6 (60). Abricate v1.0.1 (https://github.com/tseemann/abricate) was used to identify matches to the ResFinder, VirulenceFinder, and PlasmidFinder databases in the assembled contigs (61–64). Putative prophages were identified with PHASTER (65), and recombination blocks were identified with Gubbins v2.4.1 (66).
Short reads were quality trimmed with Trimmomatic v0.32 (67) and settings MAXINFO:200:0.4. MLSTs were determined from the short reads using SRST2 v0.2.0 (68) and the E. faecium pubMLST database (69). For mapping to reference genomes a core alignment was generated by mapping short reads with Snippy v4.6.0 (https://github.com/tseemann/snippy) and masking all putative transposases, complete prophage regions, and recombination blocks. For mapping to reference genomes, a core alignment was generated by mapping short reads with Snippy v4.6.0 (https://github.com/tseemann/snippy) and masking all putative transposases, complete prophage regions, and recombination blocks. Non-ACGT bases were converted to N with snippy-clean and then the variant sites were retained with snp-sites v2.5.1 (70). Maximum-likelihood phylogenies were constructed with IQTree v2.0.3 from the core SNP alignment with automatic model selection and 1000 ultrafast bootstraps (71–73). Phylogenies were visualised with iTOL (74). Short read assemblies were generated with Unicycler v0.4.8 (75) and searched for antimicrobial resistance genes using Abricate and the ResFinder database.
Transmission Network Inference
Short reads were mapped to the VRED06-10 ST80 reference chromosome with Snippy, the V24 E. faecium ST80 genome (Accession CP036151) was included as an outgroup. A core alignment was generated as above then a posterior set of phylogenies were generated with MrBayes v3.2.7 (76). Two MCMC runs consisting of four coupled chains were run for 5,000,000 generations and sampled every 5000 generations. At the completion of the run the standard deviation of split frequencies was 0.013, the log likelihood was stable, and the effective sample size of the prior, likelihood, and all numerical parameters was >800. A random sample of 100 trees from the MrBayes posterior set were used as input for Phyloscanner v1.6.6 (77). Parsimony reconstruction was performed in Phyloscanner using the Sankoff algorithm and k parameter of 281494.5, equivalent to a within patient diversity threshold of 10 SNPs. The host relationship summary was then used to construct a transmission network in Cytoscape v3.9.0 (78), only edges with complex or transmission state and >0.5 probability were visualised.
Plasmid Detection
As described by Pinholt et al (37), to detect plasmids in the entire collection, short reads were mapped to plasmids from the two complete assemblies using Snippy. Plasmids were considered present if ≥85% bases were called with <20 SNPs/1000 bp.
Statistical Analysis
To determine the optimal number of colonies to analyse for within-sample diversity a power calculation was performed as described by Huebner er al (79):
Where q = 1 – concentration of organisms, ^ = exponentiation operator, n = number of colonies sequenced, and P = probability of finding one or more variants.
Population variants were considered distinct if they differed by >10 chromosomal SNPs. Moradigaravand et al (20) have shown within-patient populations of VREfm can harbour minority variants at 20-50% of the total population based on sequencing ∼5 colonies. To detect a variant at 20% of the population with 95% confidence would require 14 colonies per sample so this was selected for analysis of faecal population diversity in this study. However, more in depth sequencing of blood culture isolates showed much lower diversity – no minor population variants from ∼10 colonies (20) – therefore single colonies from blood culture samples were sequenced in this study to identify the sequence type causing invasive disease.
Presence/absence matrices of AMR genes were generated in R v4.0.5 using ggplot2 and patchwork packages (80–82).
Data Availability
Sequence data have been uploaded to the NCBI under BioProject PRJNA877253
AUTHOR CONTRIBUTIONS
M.P.M., S.H.G., K.E.T., and M.T.G.H. conceived the project. S.H.G., K.E.T., and M.T.G.H. supervised the project. M.P.M. and K.A.P. performed genome sequencing. M.P.M. performed data analysis. M.P.M. and M.T.G.H. wrote the manuscript. M.T.G.H., S.H.G., T.J.E., and A.L. secured funding for the project. All authors reviewed and approved the manuscript.
COMPETING INTERESTS
The authors declare no competing interests
ETHICS DECLARATIONS
Access to isolates and clinical data was approved by the NHS Scotland Biorepository Network (Ref TR000126) and University of St Andrews Research Ethics Committee (Ref MD12651).
ACKNOWLEDGEMENTS
This work was supported by the Chief Scientist Office (Scotland) through the Scottish Healthcare Associated Infection Prevention Institute (Reference SIRN/10). The authors acknowledge the Research/Scientific Computing teams at The James Hutton Institute and NIAB for providing computational resources and technical support for the “UK’s Crop Diversity Bioinformatics HPC” (BBSRC grant BB/S019669/1), use of which has contributed to the results reported within this paper. Bioinformatics and Computational Biology analyses were further supported by the University of St Andrews Bioinformatics Unit which is funded by a Wellcome Trust ISSF award [grant 105621/Z/14/Z].