Population genomics and antimicrobial resistance dynamics of Escherichia coli in wastewater and river environments

Escherichia coli is amongst the dominant pan-aminoglycoside-resistant bacteria from wastewater and river water

A total of 168 bacteria highly resistant to aminoglycosides were obtained by growing collected samples in MacConkey media supplemented with high aminoglycoside concentrations, from both WWTPs (66 isolates) and rivers (102 isolates) sampled at the same period of time and located in the Barcelona region. Most of the isolates were identified as Escherichia coli, being the predominant Enterobacteriaceae species in WWTPs (50 isolates, 75.76%) and the second most prevalent in rivers (20 isolates, 19.61%). Klebsiella pneumoniae was the predominant bacterium identified in river environments (50 isolates, 49.02%), although its presence in WWTP samples was scarce (7 isolates, 10.61%). The remaining aminoglycoside-resistant species, members of genus Enterobacter, Citrobacter and Aeromonas among others, represented a low fraction in WWTPs (9 isolates, 13.64%) comparing to river environments (32 isolates, 31.37%), showing a higher bacterial diversity of natural over WWTP niches (Fig. 1, Supplementary Data 1).

Antimicrobial resistance in E. coli varies between wastewater and river water environments

Resistance to diverse antibiotic compounds was related to the sampling origin and the species of the isolates. All isolates were highly resistant to 4,6-DOS aminoglycosides, since they were selected according to their resistance to this family of antibiotics. The 16S-RMTase gene identified in most wastewater isolates was rmtB (55 isolates, 83.33%), whereas armA was the responsible in the remainder (11 isolates, 16.67%). On the contrary, 100% of isolates from river environments harbored the armA gene (Fig. 1, Supplementary Data 1). No other 16S-RMTase gene was detected neither in wastewater nor river isolates. Bacteria from wastewater samples showed higher MIC values to all tested antibiotics (e.g., 77.3% isolates resistant to cefotaxime) than bacteria from freshwater environments (e.g., 0% isolates resistant to cefotaxime), even when comparing isolates belonging to the same species. This was also the case for E. coli, where isolates from WWTPs were among the most resistant (e.g., 90% isolates resistant to cefotaxime) while those from rivers were among the most susceptible ones (e.g., 0% isolates resistant to cefotaxime), highlighting key differences between bacteria from these different ecological niches (Supplementary Data 1).

Distinct E. coli populations exist in wastewater and river water environments

Pulsed-field gel electrophoresis (PFGE) of E. coli isolates revealed the existence of two highly prevalent pulsotypes present in the two WWTPs (31 isolates, 62%), and a total of 5 characterized pulsotypes (≥90% similarity in the PFGE pattern) from wastewater samples. Interestingly, 15 isolates (30%) from both WWTPs could not be typed by PFGE. In contrast, river environments showed higher clonal diversity, with up to 10 different pulsotypes characterized, despite the lower number of E. coli isolates recovered from these samples (Supplementary Fig. 1). Of the 70 E. coli isolates, 43 were selected for further analyses and sequencing by Illumina technology (25 from WWTPs and 18 from river environments), ensuring wide representation in relation to the origins, clonal relatedness, and antibiogram resistance profiles.

Multi-locus sequence typing (MLST) results were highly correlated with PFGE profiles. Almost all E. coli isolates from wastewater samples (23 isolates, 92%) belonged to two predominant sequence types (STs) (Fig. 2). The most prevalent one, ST1196/ST632 (Warwick/Pasteur MLST scheme), encompassed the two main pulsotypes identified by PFGE. ST1196 is an increasingly prevalent ST related to OXA-48-carbapenemase production and associated with the mcr-1 colistin-resistance gene expansion in clinical settings⁹, hospital wastewater¹⁰ and companion animals¹¹. Minimal inhibitory concentration (MIC) values of the isolates belonging to this ST showed resistance to numerous antibiotic classes, including clinically critical antibiotics such as 3^rd generation cephalosporins and colistin⁶. The other predominant ST in WWTPs, ST224/ST479, comprised all the isolates which were non-typeable by PFGE. ST224 is a pandemic multi-drug resistant ST previously associated with NDM-, CTX-M- and KPC-carbapenemase production, found in both clinical human¹² and animal¹³ samples, but also in natural environments¹⁴. ST224 E. coli isolates identified exhibited a common resistance profile, including resistance to 3^rd generation cephalosporins. However, two isolates from river samples belonging to this same ST showed a different resistance pattern, revealing the influence of the niche in the bacterial resistance phenotype, even within the same ST. One of the E. coli isolates from El Prat WWTP was also identified as ST131, a clonal group present in multiple environments with a plethora of resistance mechanisms and virulence factors, representing a major public health concern¹⁵. River-related E. coli showed a higher number of E. coli STs compared to wastewater isolates, comprising up to 6 different E. coli STs (Fig. 2). All river E. coli STs presented similar MIC values for most of the tested antibiotics, which were lower than those exhibited by wastewater E. coli. ST607 (Warwick MLST scheme) was the most prevalent ST found in rivers. This ST has been scarcely reported, although it has been already detected in river sediments¹⁶. MICs of plazomicin for all isolates, both from wastewater and river environments, were ≥512 mg/L, demonstrating the high-level resistance to aminoglycosides conferred by 16S-RMTases, even to an aminoglycoside that has not yet been approved for clinical use in the EU (Fig. 2).

Resistome analysis results were highly correlated with resistance phenotypic profiles. Thus, isolates from WWTPs harbored an heterogenous antibiotic resistance gene content depending on the E. coli ST, such as the strong association between the mcr-1 gene and isolates belonging to ST1196¹⁷. This heterogenicity was even present between isolates belonging to ST1196, which exhibited different resistance levels to chloramphenicol and tetracycline depending on the presence of cmlA1/floR and tet(A), respectively. However, all wastewater isolates showed a common high-level resistance to β-lactam compounds, including third generation cephalosporins, which could be attributed to the presence of bla_CTX-M-55 (ST1196) and bla_CMY-2 (ST224) (Fig. 2). On the contrary, all isolates from river environments possessed a uniform antibiotic resistance gene content, despite comprising different STs (Fig. 2). Almost all of them (17 isolates, 94.4%) were susceptible to β-lactams, including the isolates belonging to ST224, contrasting with wastewater isolates, which have been under anthropogenic pressure, such as clinical treatments with a combination of aminoglycosides and β-lactams¹⁸, that led to this resistance associations. Likewise, plasmidome analysis revealed that the total plasmid content, based on the plasmid incompatibility groups, was heterogeneously distributed among isolates from WWTPs according to E. coli STs and closely correlated with the resistance gene content (an average of 6.36 different plasmid replicons per isolate). ST1196 isolates carried a higher plasmid content, exhibiting different plasmid profiles among them. The only plasmid incompatibility group carried by all wastewater isolates was IncFII type, specifically a pC15-1a-like plasmid (Fig. 2). Considering river E. coli STs, the total plasmid content was generally lower (an average of 3.72 different plasmid replicons per isolate) and more uniform compared to wastewater STs, similar to the pattern of the resistome (these differences are addressed in the E. coli pan-genome structure and genome, plasmid and antibiotic resistance gene diversities section). The most prevalent plasmid incompatibility group found among them was IncHI2A, present in all river isolates except ST224, which harbored a completely different plasmid content (Fig. 2).

E. coli STs have a similar genomic complexity in wastewater and river water, but the diversity of plasmids and resistance genes is higher in wastewater STs

The pan-genome of all 43 E. coli isolates, with the independence of the origin, was constituted by a total of 13,819 different genes. The pan-genome was distributed in a common core-genome of 3109 genes (22.5%) and a variable accessory-genome of 10,710 genes (77.5%), covering a considerable E. coli diversity, considering that the estimation of the global E. coli core-genome comprises around 1500 genes¹⁹ (Supplementary Fig. 2). Attending to the pan-genome configuration depending on the water source, 3340 out of a total of 9099 genes (36.71%) were included in the core-genome of wastewater isolates, and the core-genome of river isolates was formed by 3410 out of 9927 genes (34.35%), showing that river isolates presented a larger total gene pool and a smaller relative core-genome comparing to WWTP isolates. The genes conforming the pan-genomes of wastewater and river water E. coli were statistically different, considering the genes that were present and absent in each niche from the total pan-genome (Jaccard, P-value = 0.001) (Fig. 3a), indicating that the different environments, and/or upstream environments from which they have seeded, led to the selection of specific genomic populations, even between members of E. coli ST224 originating from the two different environments (Jaccard, P-value = 0.049). The genomic diversity of the whole E. coli population from river water was significantly higher than the one found in wastewater (Jaccard, P-value = 3.599 × 10⁻⁸) (Supplementary Fig. 3a). However, considering the number of different STs constituting the genomic pool of each environment, the genetic diversity of an E. coli ST from wastewater was statistically similar to the diversity of an E. coli ST from freshwater (Jaccard, P-value = 0.3123) (Fig. 3b). This model, which took into account the number of STs from each source in the diversity analysis, was previously checked by random sampling analysis, obtaining similar results. Essentially, the diversity of the total E. coli population was dependent on the number of different STs that defined the population, but the diversity of a specific ST was independent of this factor and more suitable to estimate the variability of specific bacterial clones according to the origin. Likewise, the diversity of plasmid content, taking into account the different plasmid incompatibility groups identified, was distinctive for each water type (Jaccard, P-value = 0.001) (Fig. 3c). Thus, the ecological niche also influenced the presence of specific plasmid types. The level of complexity of the E. coli plasmid pool circulating in each aquatic environment was statistically similar (Jaccard, P-value = 0.2447) (Supplementary Figure 3b). However, the number of different plasmid types carried by each E. coli ST was significantly higher in WWTPs than in rivers (Jaccard, P-value = 1.037 × 10⁻¹⁴) (Fig. 3d). Both wastewater and river E. coli populations also showed a great divergence in the antibiotic resistance gene content (Jaccard, P-value = 0.001) (Fig. 3e), showing different resistance mechanisms to particular antimicrobial classes. Furthermore, this diversity was significantly higher in the total E. coli population from wastewater comparing to the E. coli population from river water (Jaccard, P-value = 5.093 × 10⁻¹⁴) (Supplementary Fig. 3c), even when the latter possessed a higher total genomic diversity. Thus, the sum of genetic resistance determinants carried by each E. coli ST from wastewater environments was much higher than for each E. coli ST from natural effluents (Jaccard, P-value < 2.2 × 10⁻¹⁶) (Fig. 3f).

Regarding the predicted functions associated with the gene clusters from the E. coli pan-genome, 2377 different functions were identified in both wastewater and river isolates, of which 1973 (83%) were included in the functional core-genome and 404 (17%) were assigned to the functional accessory-genome (Supplementary Fig. 4). Although the accessory-genome was considerably larger than the core-genome at the gene level, most of the accessory-genome genes were involved in core-genome functions. The total genetic functions encoded in wastewater isolates was 2291, of which 2021 were shared by all of them (88.21%), similar to the 2305 genetic functions identified in river isolates, including the 2023 (87.77%) belonging to the river functional core-genome. Like the pan-genome results at genetic level, the functional diversities found in both water environments were statistically different (Jaccard, P-value = 0.001), showing that the functions are closely linked and specific for each ecological niche (Fig. 3g). Due to the greater complexity of the total genomic content, the number of different functions carried out by the entire river E. coli population was significantly higher comparing to the wastewater E. coli population (Jaccard, P-value = 2.137 × 10⁻¹²) (Supplementary Fig. 3e). Nonetheless, the degree of functional diversity, considering the different STs presented in each environment, was not statistically different (Jaccard, P-value = 0.2637) between human-related and natural water E. coli (Fig. 3h). Among functions conducted by wastewater isolates, 51 of them were specific from this niche compared to river isolates. These functions were specially related to antibiotic resistance and plasmid processes, which according to the antibiotic resistance and plasmid diversities showed higher levels in WWTPs. On the other hand, 45 functions were well-represented in river isolates but absent in wastewater environments, predominantly including functions linked to regulation processes and metal metabolism (Fig. 4). Interestingly, specific phage functions were representative for each of the ecological niches, suggesting that particular phages were associated with certain environments, as described recently²⁰.

Long- and short-read hybrid assembly allowed the complete resolution of the genomic structure of 15 representative E. coli, attending to the origin, ST, antibiotic resistance gene content, and plasmid content. The length of the chromosomes of E. coli from wastewater samples ranged between 4,674,785 bp and 4,979,626 bp, while isolates from rivers possessed chromosomes between 4,695,951 bp and 5,213,442 bp, with most of them being larger than 5 Mb. Besides, the chromosome of isolates belonging to the same E. coli ST from rivers also showed a greater genetic structural variation comparing with the chromosome of wastewater E. coli STs, which was highly conserved among their members. Most of the genes involved in this chromosome diversity of river isolates were associated with mobile genetic elements other than plasmids, specially insertion sequences and transposons from multiple families, which, in turn, could act as source of chromosome variation by favoring the integration of other genetic elements. Phage-related genes were also significant elements for chromosome variation of river E. coli STs, pointing out the abundance of this kind of elements in natural aquatic environments. Considering antibiotic resistance genes integrated in the chromosome, all isolates from both origins harbored the gene mdfA flanked by the same chromosome-related genes and no mobile genetic elements, and no other chromosome resistance mechanism was identified in almost all of the cases. MdfA is a multi-drug transporter commonly found in E. coli, which enables resistance to a broad spectrum of chemically different compounds, including antibiotics such as certain aminoglycosides and fluoroquinolones²¹. Among the 15 representative E. coli isolates, a total of 70 plasmid structures were assigned to at least one plasmid incompatibility group, of which only 8 (11.43%) showed more than one replicon type in the same structure. The majority of these hybrid plasmids were comprised of replicons belonging to the IncF type family, which are well-known for their capability to coexist in these chimera structures²². None of these hybrid plasmids harbored any 16S-RMTase gene, the genetic resistance targets of the present study. Focusing on the implications of these hybrid plasmids in the plasmidome diversity of the whole E. coli populations, the diversity analyses were based on the different replicons identified in the isolates, not the independent plasmid structures, since hybrid plasmids are the result of plasmid fusions, reflecting a higher plasmid diversity. Furthermore, the presence of more than one replicon confers a great adaptative advantage to the plasmid structure to be transferred to other bacteria and it can potentially split into independent structures²³, which supported the use of the replicons as the key marker to study plasmidome diversity.

Aminoglycoside resistance is disseminated by distinctive plasmid mechanisms in wastewater and river water

16S-RMTase genes were located in plasmids in all wastewater and river isolates. rmtB was integrated in an IncFII plasmid type, specifically a pHN7A8-like plasmid, in all isolates sequenced by long-read WGS and belonging to ST1196, which came from wastewater samples. pHN7A8 is a 76,878 bp multi-drug resistance plasmid firstly characterized in China a decade ago and closely involved in rmtB dissemination²⁴. The pHN7A8-like plasmids identified in ST1196 E. coli were 69,320 bp in length and shared 96.0% minimum pairwise nucleotide identity among them and 90.3% with the original plasmid (Fig. 5a). Focusing on the variable region of the sequenced plasmids (the region of the plasmid where genetic rearrangements are focalized because of the presence of mobile genetic elements), the gene content was identical between all pHN7A8-like plasmids, comprising two variants of bla_TEM-1, bla_CTX-M-55, and rmtB, which were flanked by four copies of IS26 (Fig. 5b). Despite the high conservation of the structure amongst our sequenced plasmids compared to the original pHN7A8 structure, the resistance gene content showed specific differences, such as the loss of fosA3 and the integration of bla_CTX-M-55 instead of bla_CTX-M-65, which could be favored by the abundant presence of IS26. These divergences were also present at SNP level comparing all the sequenced plasmids around the world with the same pHN7A8 origin of replication (46 plasmids), showing a high similarity among Barcelona plasmids, which were distantly related to all other plasmids. Almost all of them were isolated in Asia, especially China, where this plasmid type was predominantly confined and showed a common structure. Only two plasmids were located out of Asia, in Bolivia, but they were closely related to Asian plasmids with 9 SNPs of difference from the original pHN7A8, demonstrating a possible transcontinental dissemination²⁵. However, the pHN7A8 plasmids characterized in WWTPs of Barcelona presented between 90 and 103 SNPs against the Chinese plasmid (Supplementary Fig. 5a, https://microreact.org/project/yeLN26yJf).

Considering the ST224 isolates, they also harbored the rmtB gene in an IncFII plasmid type but similar to plasmid pC15-1a. This plasmid, originally described in Canada, was the first bla_CTX-M-15-carrying plasmid completely sequenced, with a size of 92,353 bp, and it is well-known for its broad geographical distribution and genetic content plasticity²⁶. pC15-1a-like plasmid identified in ST224 was 76,521 bp-sized and presented 67% nucleotide identity with the original pC15-1a. Plasmids with this same origin of replication were found in other STs identified in both wastewater, including all ST1196 isolates, and river E. coli (Fig. 5a). However, no 16S-RMTase gene was integrated in pC15-1a-like plasmids carried by these other STs, except for the ST615 isolate from WWTP water and one ST635 isolate from river water, which harbored the armA gene instead of rmtB and presented different plasmid sizes. The differences among pC15-1a-like plasmids were caused by the multiple mobile genetic elements integrated in the variable region of the plasmids, especially IS26, and were the cause of the numerous and diverse resistance gene combinations (Fig. 5b). Notably, and despite the different origins, the two armA-carrying pC15-1a-like plasmids shared a highly similar variable region, except for the absence of bla_TEM-1b in the ST635 isolate from river water, which could indicate plasmid transfer between these two niches, and the effect of the environment in the plasmid content shaping (Fig. 5b). According to the SNP-tree analysis, including worldwide plasmids with identical origins of replication (60 plasmids), the pC15-1a plasmid family diversified for a long period of time in multiple branches that disseminated in all continents, mainly associated to E. coli recovered from human samples, which explained the low nucleotide similarity between their members. pC15-1a-like plasmids identified in Barcelona waters were located in different cluster branches of the SNP-tree (Supplementary Fig. 5a). In the case of the sole ST131 isolate detected, the 16S-RMTase gene identified was armA, contrasting with the majority of wastewater isolates, and the gene was integrated in an IncR type plasmid, a different plasmid type to those responsible in the other STs and previously associated with armA dissemination²⁷.

In contrast to the E. coli from wastewater, all isolates from river environments, belonging to four different STs, harbored the armA gene encoded in a common IncHI2A type plasmid closely related to R478, except for the aforementioned ST635 isolate and one ST224 isolate that presented this gene integrated in an IncX1 type plasmid. R478 is a 274,762 bp resistance plasmid isolated in the USA in 1969 with a complex conjugative machinery²⁸. The size of the 5 R478-like IncHI2A plasmids identified in river isolates and sequenced by long-read WGS ranged between 321,526 and 334,149 bp, with ~50 kb more than the original plasmid, which resulted in a total nucleotide identity of 70.0–70.6%. However, the minimum pairwise nucleotide identity increased up to 95.1% considering only the backbone of the plasmids, and the nucleotide identity between all sequenced plasmids was 92.0% (Fig. 5c). SNP-tree analysis comprising IncHI2A plasmids isolated worldwide (37 plasmids) confirmed that all members were related, and mostly found in species belonging to genus Enterobacter such as the genetically closest one to Barcelona plasmids (47–50 SNPs), isolated in the United States of America (USA) (Supplementary Fig. 5b, https://microreact.org/project/yNTuhk1Jw). However, the R478-like plasmid identified in Barcelona has acquired a large gene content and incorporated it in the variable region of the plasmid, including the armA gene together with other multiple resistance genes, due to the abundance of mobile and integrative genetic elements located along the variable region, comprising different IS families and transposons (Fig. 5c). Remarkably, these IncHI2A plasmids also carried numerous resistance determinants to diverse metals, including the sets of genes ter, mer, and ars, which confer resistance to tellurium, mercury, and arsenic, respectively. These metal resistance mechanisms are commonly associated to IncHI2A plasmids²⁸ and they were directly linked with the metal metabolism functions specifically identified in river water isolates but not in wastewater ones, since most of the river E. coli harbored the IncHI2A plasmid and, therefore, their associated metal resistance determinants. Furthermore, the conjugative machinery of IncHI type plasmids is intimately regulated by htdA and trh genes, which are expressed in a temperature-dependent manner. At around 37 °C, htdA is expressed, and the resulting protein interacts with Trh factors, transcriptional factors of transfer-related genes, blocking their expression and repressing the conjugation process (Fig. 5c). At the optimal temperature for IncHI conjugation, between 22 and 28 °C, Trh factors are functional, thus the dissemination of this plasmid type has been related to natural water environments²⁹.

The three plasmid types carriers of 16S-RMTase genes in E. coli from aquatic environments of Barcelona (pHN7A8-like, pC15-1a-like, and R478-like plasmids) also harbored the key genetic elements that comprise a complete conjugative machinery: the relaxase or MOB gene (MOBF family for IncFII pHN7A8-like and pC15-1a-like plasmids and MOBH family for IncHI2A R478-like plasmids), the type 4 coupling protein (T4CP) gene and the type 4 secretion system (T4SS) ATPase VirB4 gene, together with 10–12 accessory genes belonging to the type F conjugative system (Figs. 5a, c). These findings confirmed the potential capacity of these plasmids to be transferred to other bacteria via conjugation, allowing the dissemination of their genetic resistance content, including the 16S-RMTase genes.

Worldwide emergence and dissemination of aminoglycoside-resistant E. coli STs are led by different evolutionary events

SNP-tree analysis of the chromosomes of the three main E. coli STs identified in Barcelona waters uncovered its relationship with other isolates belonging to these groups (40 genomes belonging to ST1196, 45 to ST224, and 11 to ST607) and its dissemination in the world (Supplementary Fig. 6, https://microreact.org/project/a9-1ZIze0). A total of 1322 SNP positions were extracted in non-recombinogenic regions of ST1196 E. coli chromosomes, of which between 0 and 5 SNPs were detected among Barcelona isolates and between 72 and 77 SNPs comparing these to the genetically closest isolate. Furthermore, the ST1196 genomic dataset exhibited a high temporal signal (correlation coefficient = 0.8102) and a constant mutational rate (Supplementary Fig. 7a), allowing the phylogeographic reconstruction of this ST. E. coli ST1196 was estimated to have emerged at the end of 1996, and the dissemination focus was located in Germany, from where this ST spread worldwide in a short period of time, especially to the USA at the beginning of 2000, which was the possible country of origin of the isolates identified in Barcelona WWTPs (Fig. 6) (See Supplementary Movie 1 and Supplementary Data 2 to visualize an interactive map in Google Earth website³⁰ by importing it as a project). We found no evidence of the spread of ST1196 from Barcelona to new locations, but new routes of dissemination from Germany, the USA, and Uganda to other countries were identified in the last two years. Most of the ST1196 E. coli have been sequenced from animal sources, mainly from poultry samples in Germany and Uganda, although the most closely related isolate to those from the Barcelona WWTPs was isolated from a plant in Florida (USA).

The ST224 E. coli chromosomes accumulated a total of 7246 SNPs along their non-recombinogenic regions. This divergence was mainly due to the existence of two clusters genetically quite distant (Supplementary Figure 6, https://microreact.org/project/a9-1ZIze0). Furthermore, the members belonging to the same cluster also presented numerous SNPs. In fact, Barcelona ST224 isolates from wastewater and river water were well-differentiated by 92–95 SNPs, indicating a possible divergent evolution in each environment from a common ancestor or an unrelated origin between them. However, WWTP isolates shared an identical chromosome (0 SNPs), whereas river isolates presented 6 SNPs between them. Focusing on the SNP analysis results of ST607 E. coli, the total SNPs detected in the chromosome dataset were 1713, but Barcelona river isolates were distributed in different clusters genetically separated by 177–181 SNPs, which highlighted the higher chromosomic variation exhibited by river isolates comparing to WWTP isolates (Supplementary Figure 6, https://microreact.org/project/a9-1ZIze0). Due to the heterogeneity among E. coli isolates, temporal signal and mutational rate of ST224 (correlation coefficient=0.1644, Supplementary Fig. 7b) and ST607 (correlation coefficient = 0.3548, Supplementary Fig. 7c) were not suitable to carry out posterior molecular clock analysis. Remarkably, no 16S-RMTase gene has been previously reported to be harbored by any of these three STs (ST1196, ST224, and ST607) (Supplementary Fig. 6, https://microreact.org/project/a9-1ZIze0).

Source: Ecology - nature.com