Exploring the mobilome and resistome of Enterococcus faecium in a One Health context across two continents

Abstract

Enterococcus faecium is a ubiquitous opportunistic pathogen that is exhibiting increasing levels of antimicrobial resistance (AMR). Many of the genes that confer resistance and pathogenic functions are localized on mobile genetic elements (MGEs), which facilitate their transfer between lineages. Here, features including resistance determinants, virulence factors and MGEs were profiled in a set of 1273 E. faecium genomes from two disparate geographic locations (in the UK and Canada) from a range of agricultural, clinical and associated habitats. Neither lineages of E. faecium, type A and B, nor MGEs are constrained by geographic proximity, but our results show evidence of a strong association of many profiled genes and MGEs with habitat. Many features were associated with a group of clinical and municipal wastewater genomes that are likely forming a new human-associated ecotype within type A. The evolutionary dynamics of E. faecium make it a highly versatile emerging pathogen, and its ability to acquire, transmit and lose features presents a high risk for the emergence of new pathogenic variants and novel resistance combinations. This study provides a workflow for MGE-centric surveillance of AMR in Enterococcus that can be adapted to other pathogens.

Fig. 1.

Overview of data-processing workflow. (a) Short reads were trimmed using fastp, with quality checks by FastQC before and after trimming. The reads were then assembled using Unicycler and quality checked using quast. (b) The quality-checked assemblies were annotated with Prokka, MOB-suite, RGI and DIAMOND. DIAMOND annotation was performed using the VFDB and BacMet databases. The Prokka annotations were passed to IslandCompare to infer probable GIs. The outputs from IslandCompare were passed to DBSCAN-SWA to infer probable phages. (c) Annotated contigs from Prokka were passed to Roary for pangenome calculation and core-genome alignment. The core-genome alignment’s static sites were removed using SNP-sites and the resultant alignment was passed to IQ-Tree to calculate a maximum-likelihood phylogenetic tree. (d) All annotated target genes and MGEs were tabulated and passed to BayesTraits for co-evolutionary analysis, performing hypothesis testing for correlated evolution between pairs of features. Acronyms used in the study are summarized in Table S8.

Fig. 2.

Size and count distributions of genomic features. First row: genome size (a) and Pan-genome distribution predicted by Roary (b). (c–h) Histograms showing the abundance of predicted genes and MGEs across the set of all genomes. Second row: frequency distribution of AMR genes (c), HMR genes (d) and VFs (e). Third row: frequency distribution of plasmids (f), GIs (g) and phages (h). Multiple occurrences of a feature in a given genome were counted only once. For clarity, only features detected in at least five genomes were plotted. Plot annotations indicate the number of features plotted and the number of total features detected.

Fig. 3.

Abundance of features by habitat type and geographic location. ‘AB’ indicates genomes sampled from Alberta, Canada. ‘UK’ indicates genomes sampled from the UK. Counts indicate the number of unique features of a given category found per genome. Bars indicate quartiles. Points/diamonds are considered to be outliers if they fall outside 1.5×the interquartile range. Grey bars indicate mean values.

Fig. 4.

Maximum-likelihood core-genome phylogenetic tree of 1273 E. faecium genomes with E. hirae ATCC9790 as the outgroup and E. faecium DO ASM17439v2 as the reference genome. The tree was constructed with 1854991 nucleotide sites, 79440 of which were parsimony informative, using the general time reversible substitution model with invariant sites and four Gamma rate categories. Branch lengths are log-transformed and scaled down to 13 % length for improved readability. Nodes are coloured by sampling location, with hue indicating habitat and saturation indicating geography.

Fig. 5.

Statistical associations and physical localization of vanA (a–d) and vanB (e–h) genes. (a, e) Phylogenetic distribution of van genes and other features with an associated likelihood ratio ≥100. (b, f) Statistical association network of vanA/vanB genes with other features. Gene and MGE colours are consistent with those in Fig. 2. (c, g) Example of gene order on an annotated plasmid (c) and GI (g). Green genes correspond to ‘Perfect’ matches with reference genes in the CARD database, yellow genes are ‘Strict’ hits. (d–h) Distribution of genes by habitat. Bar colours correspond to their habitats as per the legend in (a, e).

Fig. 6.

Heatmap showing the presence of AMR determinant genes detected in 1273 E. faecium genomes analysed in this study. The y-axis indicates genomes (colour coded by habitat, geography and type) sorted by topology of the core-genome maximum-likelihood tree. AMR determinants (x-axis) are sorted by drug class. * denotes variant versions of intrinsic genes conferring AMR.