Genomic rearrangements uncovered by genome-wide co-evolution analysis of a major nosocomial pathogen, Enterococcus faecium.

Enterococcus faecium is a gut commensal of the gastro-digestive tract, but also known as nosocomial pathogen among hospitalized patients. Population genetics based on whole-genome sequencing has revealed that E. faecium strains from hospitalized patients form a distinct clade, designated clade A1, and that plasmids are major contributors to the emergence of nosocomial E. faecium. Here we further explored the adaptive evolution of E. faecium using a genome-wide co-evolution study (GWES) to identify co-evolving single-nucleotide polymorphisms (SNPs). We identified three genomic regions harbouring large numbers of SNPs in tight linkage that are not proximal to each other based on the completely assembled chromosome of the clade A1 reference hospital isolate AUS0004. Close examination of these regions revealed that they are located at the borders of four different types of large-scale genomic rearrangements, insertion sites of two different genomic islands and an IS30-like transposon. In non-clade A1 isolates, these regions are adjacent to each other and they lack the insertions of the genomic islands and IS30-like transposon. Additionally, among the clade A1 isolates there is one group of pet isolates lacking the genomic rearrangement and insertion of the genomic islands, suggesting a distinct evolutionary trajectory. In silico analysis of the biological functions of the genes encoded in three regions revealed a common link to a stress response. This suggests that these rearrangements may reflect adaptation to the stringent conditions in the hospital environment, such as antibiotics and detergents, to which bacteria are exposed. In conclusion, to our knowledge, this is the first study using GWES to identify genomic rearrangements, suggesting that there is considerable untapped potential to unravel hidden evolutionary signals from population genomic data.

Data Summary

Raw core-genome alignment (1.1 MB, Harvest suite v1.1.2), including the 1644 clade A isolates and the complete AUS0004 (accession number CP003351) as a reference, is available under the following GitLab repository https://gitlab.com/sirarredondo/efm_gwes.

has emerged as an important nosocomial pathogen around the world. Population genetics revealed that clinical strains form a distinct clade, designated clade A1, and that plasmids are major contributors to the emergence of nosocomial . Here, the adaptive evolution of was further explored using an unsupervised machine learning method (SuperDCA) to identify genome-wide co-evolving single-nucleotide polymorphisms (SNPs). We identified three genomic regions harbouring large numbers of SNPs in tight linkage that are separated by a large chromosomal distance in a clinical clade A1 reference isolate, but appeared adjacent to each other in non-clade A1 isolates. We identified four different types of large-scale genomic rearrangements and in all cases we found insertion of two different genomic islands and an insertion element at the border. In contrast, no genomic rearrangement and insertions were identified among a group of clade A1 pet isolates, suggesting a distinct evolutionary trajectory. Based on the in silico predicted biological functions, we found a common link to a stress response for the genes encoded in three regions. This suggests that these rearrangements may reflect adaptation to the stringent conditions in the hospital environment, such as antibiotics and detergents, to which bacteria are exposed.

Introduction

are commensals of the gastrointestinal tract but are now recognized as a major causative agent of healthcare-associated infections [1]. The transition from commensal to nosocomial pathogen coincided with increased resistance to antibiotics [2]. Since the early 1980s, first gained high-level resistance to ampicillin, followed by resistance to aminoglycosides, fluoroquinolones and glycopeptides, particularly vancomycin [3, 4]. Previous whole-genome sequencing (WGS) studies identified that the population can be divided into two lineages, including a hospital-associated clade (clade A) and a community-related clade (clade B) [5, 6]. Later, clade A was subdivided into clade A1, representing the majority of hospital-associated isolates, and clade A2 for animal-related isolates [7], although two studies that included a larger collection of animal-related isolates suggested that these isolates clustered in polyphyletic groups and not in a distinct clade A [8, 9]. In a recent study also including 55 isolates from Latin America, phylogenomic analysis suggested that the animal isolates represent multiple lineages that diverged prior to the emergence of the clinical subclades in clade A [10].

Recently, we determined the plasmid content (plasmidome) of 1644 isolates from different sources, countries and years using short- and long-read whole-genome sequencing technologies in combination with machine-learning classifiers [9, 11]. This analysis revealed that the hospital-associated isolates carried a larger number of plasmid sequences compared to isolates from other sources and different configurations of plasmidome populations in the hospital environment. In addition, the source specificity was determined for the chromosomal, plasmid and whole-genome components. It was concluded that plasmid sequences make the greatest contribution to source specificity.

In this paper, WGS of the same 1644 isolates was used with the aim of identifying signals of selection acting to shape co-evolution of single-nucleotide polymorphisms (SNPs). Co-evolved SNPs can facilitate adaptation to different environments when sequentially selected mutations in adaptive elements decrease the actual fitness costs of individual mutations (antagonistic epistasis) and thus have a beneficial effect on fitness [12]. For this purpose, we used the Direct Coupling Analysis tool (SuperDCA) as previously described for , where previously undetected epistatic interactions related to e.g. survival of the pneumococcus at lower temperatures were discovered [13]. For this study, SuperDCA was applied on the core-genome alignment of to identify likely candidates of sequentially selected or coupled mutations not in strong linkage disequilibrium due to chromosomal proximity [13].

Results

Co-evolution of core genomic variation

SuperDCA was applied to a core-genome alignment generated using the Harvest suite, including 1644 clade A isolates [Table S1 (available in the online version of this article), accession number PRJEB28495] [9] and the complete genome of strain AUS0004 (accession number CP003351) as reference [14]. All locus tags mentioned below refer to AUS0004. SuperDCA identified 262 877 significant couplings between SNPs with a distance of >1 kbp (Table S2). We generated a top-10 list of SNP positions with the highest number of links using non-overlapping bins of 100 bp (Fig. S1a, b). To obtain an insight into the distribution on the AUS0004 core genome, all links for each of these 10 SNP positions were plotted (Fig. S1c). In this paper, we will mainly focus on the SNP positions from bins 3 and 10 (Fig. S1b) representing genes EFAU004_02176 and EFAU004_02173, because we were able to identify a possible/plausible biological explanation. For both bins, we observed a similar pattern of coupled SNPs with the genes EFAU004_00665 to EFAU004_00675 located at a large chromosomal distance of around 500 kbp in the AUS0004 core genome (Fig. S1c). Detailed examination of the total list of coupled SNPs revealed similar patterns of coupled SNPs for a larger region encompassing genes EFAU004_02173 to EFAU004_02178, referred to as region-1 (Fig. 1a). There were in total 2323 coupled SNP positions, i.e. between 181 different SNP positions from region-1 and 131 different SNP positions in genes EFAU004_00665 to EFAU004_00675, referred to as region-2 (Figs 1a and S1c, Table S3). In addition, all coupled SNPs for region-1 and region-2 were retrieved from the total list. These 17 236 linked SNPs were plotted on the core genome to obtain an insight into the distribution of the links (Fig. 1a). This revealed a cluster of a total of 142 linked SNPs, represented by 32 different SNP positions in region-2 and 25 SNP positions in genes EFAU004_02092 to EFAU004_02100, representing region-3, located close to region-1 (Fig. 1a, Tables S3 and S4).

Fig. 1.

(a) All 17 236 epistatic links of the genes contained in region-1, region-2 and region-3. (b) Genomic organization and annotation of genes from regions 1–3 in AUS0004.

In order to elucidate a possible biological explanation for the identified links between region-1, region-2 and region-3, we first determined putative biological functions for the proteins in these regions using a homology search for similarity with other proteins and investigated domains with known function.

Putative biological functions for proteins in region-1

Region-1 contains five protein-encoding genes and a transposase with no links (Fig. 1b). For two proteins a function prediction is challenging, i.e. EFAU004_02173 is annotated as a membrane protein, but there are only domains of unknown function (DUF), while EFAU004_02175 is annotated as a permease with unknown specificity. In contrast, EFAU004_02176 contained conserved domains belonging to the thioredoxin-like family (Trx-like) and EFAU004_02177 contained conserved domains belonging to the Rex (Rex-like) family of transcriptional regulators (Fig. 1b). Trx-like and Rex-like are likely involved in redox homeostasis of the bacterial cell, which was found to be critical for DNA synthesis and defence against oxidative stress [15]. The exact function of the ABC transporter (EFAU004_02178) is difficult to predict despite the presence of several domains. The protein contains two so-called AAA domains and a leucine-zipper (bZIP) domain. Proteins with AAA domains are members of a conserved family of ATP-hydrolyzing proteins with all kind of activities in many cellular pathways, including replication, DNA and protein transport, transcriptional regulation, ribosome biogenesis, membrane fusion and protein disaggregation or degradation [16]. bZIP domains are known to be involved in transcriptional regulation, e.g. in stress conditions such as heat in [17, 18], or abiotic stress in eukaryotes, e.g. plants [19]. Furthermore, homology search against the GenBank database to non- revealed co-localization of a similar ABC transporter and rexA gene in other species such as (strain 2842STDY5753625, accession number FEGY01000003), (strain DK-PW-092, accession number LBKE01000082) and (strain CFSAN060067, accession number AABAZK010000087), suggesting that this ABC transporter might also be involved in some kind of stress response.

Putative biological functions for proteins in region-2

The genes from region-2 can be split into two distinct gene clusters based on their putative biological functions (Fig. 1b). The first cluster encompasses genes EFAU004_00665 to EFAU004_00670, which, based on their annotation, are putatively involved in respiration and coenzyme A biosynthesis (Fig. 1b). Although SNPs linked with region-1 were observed with all these genes, we will focus on two specific genes, i.e. EFAU004_00665 and EFAU004_00670, as they contained the majority of linked SNPs (Tables S3 and S4). Very similar conserved domains for NAD(P)H : quinone oxidoreductase were identified in both proteins. EFAU004_00670 contains transmembrane helices and is therefore likely membrane bound, while EFAU004_00665 lacks these transmembrane helices and is therefore likely to be soluble. The conserved domains in EFAU004_00670 had a high similarity to YhdH of (overall amino acid similarity of 48 %). The conserved domains in EFAU004_00665 were very similar to QorEc of (overall amino acid similarity of 32 %) and QorTt Thermus thermophiles (overall amino acid similarity of 32 %). In , expression of a gene cluster containing two qor-like genes was induced under oxidative stress conditions [20]. These two Qor proteins (SA1988 and SA1989) are predicted to be soluble and SA1988 contained similar conserved domains with an overall amino acid similarity of 26 % with EFAU004_00665. This might suggest that this protein could also be involved in stress response.

The second gene cluster contains genes EFAU004_00671 to EFAU004_00675 and is likely organized as an operon (Fig. 1b). EFAU004_00671 contains a helix–turn–helix motif and is therefore predicted to be a transcriptional regulator. Its location upstream of EFAU004_00672 to EFAU004_00675 suggests that it regulates the expression of these genes. EFAU004_00672 to EFAU004_00675 display similarity with a polyamine transport system, as described for S. pneumoniae, including an ATP-binding protein, PotA, two permeases, PotB and PotC, and a substrate-binding protein, PotD [21]. Polyamines are polycationic molecules and are required for optimal growth in both eukaryotic and prokaryotic cells and are implicated in pathogenicity of . In polyamines are pivotal in survival strategies in the host when bacteria are confronted with stress conditions such as temperature shock, oxidative stress, or choline limitation [21].

Putative biological functions for proteins in region-3

Genes EFAU004_02092 to EFAU004_02100 are annotated to encode a membrane bound V-type ATPase (Fig. 1b). An identical V-type ATPase has been studied in detail in (100 % amino acid identity) [22]. The V-type ATPase belongs to the family of proton pumps and is involved in the translocation of Na+ or H+ over the cell membrane by using the energy of ATP. In the V-type ATPase appeared to be highly expressed under stress conditions such as high pH and plays an important role in sodium homeostasis under these conditions.

Genomic rearrangements in completed genomes

For 48 of 1644 isolates, a fully assembled circular chromosome was available (accession number PRJEB28495) as described previously [11]. This allowed us to assess the genomic organization in these isolates. Using the non-clade A1 isolate E0139 as the reference, 4 different types of large chromosomal rearrangements were observed among 38 isolates. Type-1 (n=1) involves an inversion of 0.73 Mbp; type-2 (n=34) an inversion of 1.2 Mbp; type-3 (n=2) two inversions of 0.38 Mbp; and type-4 (n=1) an inversion of 0.12 Mbp (Figs 2a–d). Except for the type-1 inversion, all these genomic rearrangements were exclusively observed in clade A1 isolates (Fig. 3). In clade A1 isolates we detected two previously described genomic islands, inserted adjacent to the genomic rearrangements that were described as being enriched among clade A1 hospital isolates (Fig. 2b–d) [23, 24]. One genomic island, putatively encoding a carbohydrate transport and metabolism pathway [23], was always found to be located downstream of the ABC transporter (EFAU004_02178) of region-1 (Fig. 4a). The other island encoding a phosphotransferase system (PTS) [24] was always found to be located downstream potD (EFAU004_00675) of the polyamine transport system of region-2 (Fig. 4a). In all cases, the genomic rearrangement was flanked by an IS30-like transposon downstream of a methionine synthase (metE) and a pyridine nucleotide-disulfide oxidoreductase (pyr) (Fig. 4a). In the 10 isolates for which a fully assembled circular chromosome was available and that lacked the genomic rearrangement, insertions of both genomic islands and IS30-like transposon, region-1 and region-2 are located adjacent to each other as schematically indicated for the non-clade A1 isolate E0139 (Fig. 4b).

Fig. 2.

Genome comparisons using strain E0139 (non-clade A1) as a reference. (a) E0139 compared to type 1 strain E1334; (b) E0139 compared to type 2 strain E8202; (c) E0139 compared to type 3 strain E6055; (d) E0139 compared to type 4 strain E7098. ‘n’ indicates the number of isolates with a similar genomic organization and for which a complete chromosome was available. Same-strand DNA similarity is shaded red, while reverse similarity is shaded blue. Blue bar, position of the pot operon; green bar, position of the ABC transporter; red arrow, insertion site for the phosphotransferase system (PTS)-encoding genomic island [24]; brown arrow, insertion site for the carbohydrate transport system-encoding genomic island [23].

Fig. 3.

Core-genome tree based on 1644 strains with the 3 metadata panels: (i) distribution of genomic rearrangements types among 38 complete chromosomes, (ii) indication of strains lacking a genomic rearrangement among 10 complete chromosomes and draft genomes and (iii) indication of strains with insertion of a carbohydrate transport system encoding genomic island, orange [23], or a phosphotransferase system (PTS) encoding genomic island, purple [24]. Black cross, indication of clade A1 dog isolates that lack the insertions/inversion or non-clade A1 dog isolates that do contain insertions/inversion.

Fig. 4.

(a) Genomic organization of the most predominant type 2 genomic rearrangement in strain E8202 with indication of insertion and inversion sites. (b) Genomic organization for reference strain E0139.

We determined the presence of the two genomic islands and their insertion sites in the draft genomes of the remaining 1596 isolates and confirmed that the genomic islands were inserted at the same position in 93 % of the clade A1 isolates (Fig. 3). In addition, the insertions were also identified in 57 (14 %) non-clade A1 isolates, which mainly represented dog isolates (n=43) (Fig. 3). In contrast, a branch in clade A1 consisting of 68 genetically closely related isolates, including 55 pet, 9 isolates from hospitalized patients and 4 from non-hospitalized persons, lacked both genomic island insertions (Fig. 3). Based on these draft genomes, we cannot determine whether or which genomic rearrangement could have occurred, though we always observed a contig break downstream of metE and pyr, strongly suggesting the presence of an IS element on that position. In isolates that lacked insertions of the two genomic islands the genomes were organized in a manner comparable to the non-clade A1 isolate E0139, i.e. without inversion (Figs 3 and 4b).

Discussion

In this study, we identified large genomic rearrangements, enriched in hospital-associated clade A1 strains, which were uncovered by the appearance of co-evolved SNPs in three chromosomal regions; region-1, region-2 and region-3. These chromosomal regions were located on the borders of the genomic rearrangement and therefore adjacent to each other in the majority of non-clade A1 isolates, but at a large distance in clade A1 isolates. In all cases, the genomic rearrangement is flanked by an IS30-like transposon downstream from metE and pyr, which suggests that the rearrangement is the result of recombination between two IS30-like elements. In addition, upstream of metE and pyr, we identified the insertion of two, clade A1 enriched, genomic islands, encoding a PTS and carbohydrate transport system, respectively [23, 24]. For both genomic islands, it has been suggested that they may provide a fitness advantage, which was confirmed for PTS in an in vivo mouse colonization model [24]. Recently, Yan et al. reconstructed the evolution of the marine photosynthetic micro-organism based on genome rearrangements and identified rearrangement hotspots that were all in the vicinity of genomic islands [25]. In fact, the authors suggest that the genomic islands serve as hotspots that induce genome rearrangement. In addition, they observed that different clades shared a conserved backbone, but also contained clade-specific regions, which were associated with ecological adaptations. Further, for , a Gram-negative bacterium that can be found in different environments, it was shown that horizontal gene transfer played a key role in adaptation process, as many of the niche-specific functions were found to be encoded on clearly defined genomic islands and were linked to genomic rearrangements [26]. The observation that among the clade A1 isolates only the pet isolates do not contain the two genomic islands and genomic rearrangement may suggest that these elements are not advantageous in the niche represented by the gut of pet animals.

The predominant chromosomal rearrangement, type 2, identified in 71 % of the fully assembled genomes from the current study, is also present in other publicly available completed genomes, such as DO (NC_017960.1) [27], Aus0085 (NC_021994.1) [28] or E39 (NZ_CP011281.1). For the isolates in our study for which we only have a draft genome it is difficult to determine whether genomic rearrangements have occurred. However, in isolates in which metE or pyr were not adjacent to each other, as in the non-clade A1 configuration, we always identified an IS30-like element downstream from the metE or pyr gene at the border of a contig, suggesting a potential genomic rearrangement. A chromosomal rearrangement in was also observed in the first completely sequenced isolate AUS0004, where it was hypothesized that it occurred between two phage elements, resulting in a replichore imbalance [14]. However, we did not observe a replichore imbalance in any of the 38 completely assembled genomes with large chromosomal rearrangements (data not shown), which suggests that AUS0004 is an exception in that regard.

Genomic rearrangements or inversions are not rare among prokaryotic genomes, as was investigated by Repar et al. [29]. The authors developed an alignment-based approach to systematically assess inversion symmetry between pairs of prokaryotic genomes, including the genera and , which belong to the phylum Firmicutes like E. faecium. Bacillus and are dominated by X-shaped inversions that are symmetric around the ori–ter axis, similar to the type 1–3 inversions from this study.

The gene clusters located in the chromosomal regions-1, -2 and -3 that harboured co-evolved SNPs clearly encode common functions. The thioredoxin pathway, Rex-like (region-1) and the NAD(P)H-dependent oxidoreductase/epimerase proteins (region-2) are all dependent on the oxidation/reduction of NAD+/NADH and involved in redox homeostasis, while the polyamine transport system (region-2) and V-type ATPase (region-3) are both involved in ATP hydrolysis/synthesis. The fact that in non-clade A1 isolates chromosomal regions-1,-2 and -3 are located adjacent to each other explains the high number of statistically significantly linked SNPs between the genes in these two regions, as they likely have co-evolved in fairly tight linkage disequilibrium (LD). Conversely, this finding illustrates that the co-evolutionary analysis of SNPs can uncover inversions in the chromosome. Enterococci are rich in mobile genetic elements harbouring IS elements and transposons, many of which can integrate into the chromosome [3, 30] and may drive the observed rearrangements. These rearrangements may play a role in adaptation to new environments, more particularly adaptation to stress conditions, like the hospital environment, where bacteria are exposed to different stressors, including high concentrations of antibiotics, detergents and antiseptics. Furthermore, it is possible that the rearrangements are associated with the rapid adaptation of the organism to survive in the perturbed microbiota of hospitalized patients and that survival in the overall context of the intestinal microbiome may be an important driver of evolution. The lack of comparable selection pressure is a plausible explanation for the observed absence of genomic arrangements among the pet isolates from clade A1.

To our knowledge, identification of genomic rearrangements by direct coupling analysis of SNPs has not previously been described in the literature. The expected rapid increase in the availability of long-read sequences combined with large population-based collections of short-read sequenced genomes that are amenable to a statistical LD analysis such as performed here will thus facilitate the development of new methods to identify candidates for inversions associated with selective advantage.

Methods

Genomic DNA sequencing and assembly

A detailed description of Illumina and ONT sequencing has been provided previously and includes a full description of ONT selection of isolates (n=62) [9, 11] and consecutive hybrid assembly using Unicycler [31]. The R package ggtree (version 1.14.6) [32] was used to plot distinct metadata panels together with a previously described core genome tree of 1644 strains [9].

Co-evolution analysis

A core-genome alignment (1.1 MB) was generated using the Harvest suite v1.1.2, including the 1644 clade A isolates and the complete AUS0004 (accession number CP003351 [14]) as reference. SuperDCA was performed with the default settings as described by Puranen et al. [13] using the alignment of 1644 isolates, which was unfiltered for recombination events. Circos was used for the visualization of coupled SNPs on the genome [33].

Prediction of putative biological functions of proteins from regions 1–3

To predict putative biological functions of proteins, the protein sequences were compared with the non-redundant public database using the National Center for Biotechnology Information (NCBI) blastp server (www.ncbi.nlm.nih.gov/BLAST/) [34]. In addition, a graphical summary was used to determine whether putative conserved domains were identified and, if present, this graphical summary was selected to investigate the list of conserved domains in more detail in order to make a prediction regarding the putative biological function of the protein.

genome organization

To unravel rearrangements in the chromosome, we considered a non-clade A1 isolate from a non-hospitalized person (E0139) with a complete chromosome sequence as reference. We only considered isolates with a complete circular chromosome (n=48). Pairwise chromosomal comparisons were computed using blastn (version 2.7.1+) and alignments were visualized using the Artemis Comparison tool (version 17.0.1). The insertions of two genomic islands encoding a PTS [24] and carbohydrate transport system [23] and their putative co-localization with the pot operon and ABC transporter, respectively, were determined by blast. Additionally, for each chromosomal rearrangement, we generated dotplots using Gepard (version 1.40) [35] with the non-clade A1 isolate E0139 as a reference. Clone manager 9 was used to visualize the genomic organization of regions 1–3.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.