1. Microbiology and Infectious Disease

Mobilome-driven segregation of the resistome in biological wastewater treatment

  1. Laura de Nies
  2. Susheel Bhanu Busi
  3. Benoit Josef Kunath
  4. Patrick May
  5. Paul Wilmes  Is a corresponding author
  1. Systems Ecology Group, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Luxembourg
  2. Bioinformatics Core, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Luxembourg
  3. Department of Life Sciences and Medicine, Faculty of Science, Technology and Medicine, University of Luxembourg, Luxembourg
https://doi.org/10.7554/eLife.81196
Share this article
Cite this article
  1. Laura de Nies
  2. Susheel Bhanu Busi
  3. Benoit Josef Kunath
  4. Patrick May
  5. Paul Wilmes
(2022)
Mobilome-driven segregation of the resistome in biological wastewater treatment
eLife 11:e81196.
https://doi.org/10.7554/eLife.81196
  2. Download BibTeX
  3. Download .RIS

Abstract

Biological wastewater treatment plants (BWWTP) are considered to be hotspots for the evolution and subsequent spread of antimicrobial resistance (AMR). Mobile genetic elements (MGEs) promote the mobilization and dissemination of antimicrobial resistance genes (ARGs) and are thereby critical mediators of AMR within the BWWTP microbial community. At present, it is unclear whether specific AMR categories are differentially disseminated via bacteriophages (phages) or plasmids. To understand the segregation of AMR in relation to MGEs, we analyzed meta-omic (metagenomic, metatranscriptomic and metaproteomic) data systematically collected over 1.5 years from a BWWTP. Our results showed a core group of 15 AMR categories which were found across all timepoints. Some of these AMR categories were disseminated exclusively (bacitracin) or primarily (aminoglycoside, MLS and sulfonamide) via plasmids or phages (fosfomycin and peptide), whereas others were disseminated equally by both. Combined and timepoint-specific analyses of gene, transcript and protein abundances further demonstrated that aminoglycoside, bacitracin and sulfonamide resistance genes were expressed more by plasmids, in contrast to fosfomycin and peptide AMR expression by phages, thereby validating our genomic findings. In the analyzed communities, the dominant taxon Candidatus Microthrix parvicella was a major contributor to several AMR categories whereby its plasmids primarily mediated aminoglycoside resistance. Importantly, we also found AMR associated with ESKAPEE pathogens within the BWWTP, and here MGEs also contributed differentially to the dissemination of the corresponding ARGs. Collectively our findings pave the way toward understanding the segmentation of AMR within MGEs, thereby shedding new light on resistome populations and their mediators, essential elements that are of immediate relevance to human health.

Editor's evaluation

This paper reports important results regarding the presence and potential dissemination of antibiotic resistance genes in wastewaters by convincingly combining analysis of gene abundance, expression, and association with mobile genetic elements and bacterial taxa. Via systematic evaluation and implementation of multiple tools, the authors provide a valuable approach for monitoring antibiotic resistance genes in the environment and assessing their dispersal and possible risks to human health.

https://doi.org/10.7554/eLife.81196.sa0

Introduction

Throughout human history, bacterial infections have been a major cause of both disease and mortality (Bonilla and Muniz, 2009). The discovery as well as the subsequent development and medical use of antibiotics have provided effective treatment options which limited the development and spread of bacterial pathogens. However, the use of antibiotics has exacerbated the emergence of antimicrobial resistance (AMR) in both commensal and pathogenic bacteria (Wright, 2007). As a result, AMR, as the ‘silent pandemic’, has become a prevalent threat to human health (Brogan and Mossialos, 2016; Mahoney et al., 2021; O’Neill, 2014).

From a public health perspective, biological wastewater treatment plants (BWWTPs) are considered hotspots of AMR due to the convergence of antibiotics with resistant, potentially pathogenic microorganisms originating from both the general population as well as agriculture, healthcare services and industry (Alexander et al., 2020; Rodríguez-Molina et al., 2019). Additionally, the mobilization of antimicrobial resistance genes (ARGs) through rampant horizontal gene transfer (HGT) promotes the dissemination of AMR within the BWWTP microbial community (von Wintersdorff et al., 2016). Therefore, BWWTPs represent an environment exceptionally suited for the evolution and subsequent spread of AMR (Calero-Cáceres et al., 2014; Chen et al., 2013). To date, more than 32 studies have documented the role of BWWTPs as key reservoirs of AMR (Fouz et al., 2020). Furthermore, BWWTPs generally do not contain the necessary infrastructure to remove either ARGs or resistant bacteria, which are released into the receiving water via the effluent, promoting its spread in the environment at large (Alexander et al., 2020). Most often these are surface water bodies such as rivers, which contribute to the further dissemination of AMR and resistant bacteria among environmental microorganisms (Singer et al., 2016). Acquired resistance may in turn be carried over to humans and animals using these water resources. In fact, there is strong evidence suggesting that ARGs from environmental bacteria can be taken up by human-associated and pathogenic bacteria (Nadeem et al., 2020; Trinh et al., 2018). From an epidemiological and surveillance perspective, BWWTPs also provide samples representative of entire populations (Hendriksen et al., 2019). As such, BWWTPs have recently been crucial for the monitoring of SARS-CoV-2 within the human population (Herold et al., 2021). Overall, to increase our understanding of the dissemination of AMR and the underlying mechanisms as well as its general prevalence, it is necessary to map the resistome of various environments starting with biological BWWTPs because it is critical to unravel the extent to which they act as reservoirs for the dissemination of antimicrobial resistance genes (ARGs) to bacterial pathogens. Moreover, understanding the community-level overviews of the ARG potential and its expression, coupled with population-level linking, including to pathogens, may allow for efficient monitoring of pathogenic and AMR potential with broad impacts on human health.

The presence of resistance genes and mobile genetic elements (MGEs) along with sub-inhibitory antibiotic selection pressures may facilitate HGT of ARGs into new hosts through the mobilome (von Wintersdorff et al., 2016). Previous work has in particular shown that antibiotic selection pressures may alter HGT processes, thereby increasing the number of resistance elements which reside on mobile DNA (Datta and Hughes, 1983). Acquisition of ARGs via MGEs primarily occurs through two mechanisms: conjugation or transduction (MacLean and San Millan, 2019). In conjugation, plasmids carrying one or more resistance genes are transferred between microorganisms (Carattoli, 2013), while in transduction bacteriophages carrying ARGs infect bacteria and integrate their genome into those of the host thereby conferring resistance (Chiang et al., 2019). Of these mechanisms, conjugation is often thought to have the greatest influence on the dissemination of ARGs, while transduction is deemed less important (von Wintersdorff et al., 2016). In general terms, studies concerning AMR and its dissemination focus either on phage (Lood et al., 2017; Strange et al., 2021) or plasmids solely (Li et al., 2019). Alternatively, the two are treated collectively (Alexander et al., 2020; Che et al., 2019) without a comprehensive comparative analysis. This circumstance has created a knowledge gap whereby the contributions of plasmids and phages as independent entities to AMR transmission within complex communities, such as those found in biological BWWTPs, is largely unknown.

To shed light on the dissemination and potential segregation of AMR within MGEs in a WTTP microbial community, we leveraged longitudinal meta-omics data (metagenomics, metatranscriptomics, and metaproteomics). Samples collected for 51 consecutive weeks over a period of 1.5 years were used to characterize the resistome. We found that several bacterial orders such as Acidimicrobiales, Burkholderiales, and Pseudomonadales were associated with 29 AMR categories across all timepoints. Our longitudinal analysis suggests that MGEs are important drivers of AMR dissemination within BWWTPs. More importantly, we reveal that MGEs, that is plasmidomes and phageomes, contribute differentially to AMR dissemination. Furthermore, we observed this phenomenon in clinically-relevant taxa such as the ESKAPEE pathogens (Reza et al., 2019), for which plasmids and phages were exclusively associated with specific ARGs. Collectively, our data suggest that BWWTPs are critical reservoirs of AMR which show clear evidence for the segregation of distinct AMR genes within MGEs especially in complex microbial communities. In general, we believe that these findings may provide crucial insights into the segregation of the resistome via the mobilome in any and all reservoirs of AMR, including but not limited to animals, humans, and other environmental systems.

Results

Longitudinal assessment of the resistome within a BWWTP

To characterize the BWWTP resistome, we sampled a municipal BWWTP on a weekly basis over a 1.5-year period (ranging from 21-03-2011 to 03-05-2012) (Herold et al., 2020; Martínez Arbas et al., 2021). Metagenomic and metatranscriptomic reads were preprocessed, and both sets of reads co-assembled using the Integrated Meta-omic Pipeline as described previously (Narayanasamy et al., 2016). Subsequently, utilizing the PathoFact pipeline (de Nies et al., 2021) on the assembled contigs (Methods), we resolved the BWWTP resistome. This analysis revealed the presence of 29 different categories of AMR within the BWWTP. Subsequent longitudinal analyses highlighted enrichments in aminoglycoside, beta-lactam and multidrug resistance genes (Figure 1a). Concomitantly, we observed specific shifts in the AMR profiles over time. For example, a transient change at two timepoints (13-05-2011, 08-02-2012) highlighted a steep increase in resistance genes corresponding to glycopeptide resistance. Other AMR categories, such as diaminopyrimidine resistance, exhibited a less drastic but more fluid change in longitudinal abundance observable over multiple timepoints.

Figure 1 with 1 supplement see all
Download asset Open asset
Longitudinal metagenomic and metatranscriptomic assessment of AMR.

(a) ARG relative abundances over time within the BWWTP. (b) AMR categories at various timepoints categorized in four distinct groups based on presence/absence: Core (all timepoints), Prevalent (>75% of timepoints), Moderate (50–75% of timepoints), and Rare (<50% of all timepoints). (c) Relative abundance levels of expressed AMR categories over time within the BWWTP. Colors of all panels correspond to the AMR categories.

Additionally, AMR categories were found to persist over time within the BWWTP (Figure 1b). A core group of 15 AMR categories in total were identified and found to be present across the 1.5-year sampling period. These included aminoglycoside, beta-lactam, and multidrug resistance genes, which contributed the most to the pool of ARGs. A further six (aminocoumarin, aminoglycoside:aminocoumarin, elfamycin, nucleoside, triclosan, and unclassified) AMR categories were found to be prevalent (>75% of all timepoints), while another three AMR categories were moderately (50–75% of all timepoints) present over time (Figure 1b). Five other categories were rarely present within the BWWTP, with resistance corresponding to acridine dye only present at six of the timepoints. Altogether, this emphasized that the BWWTP resistome varies over time, substantiating the requirement for a longitudinal analysis to obtain an accurate overview of the community’s overall resistome.

Although the data thus far provided a clear overview of the BWWTP from a metagenomic perspective, it did not provide any information regarding AMR expression. We therefore utilized the corresponding metatranscriptomic dataset to investigate the expression of identified ARGs and monitor their changes, within the BWWTP, over time. In contrast to the metagenomic data, we observed a difference in AMR expression levels for several categories. Aminoglycoside, beta-lactam, and multidrug resistance identified at high levels in metagenomic information were also highly expressed within the BWWTP (Figure 1c). However, peptide resistance demonstrated the highest expression levels of all the AMR categories. We further investigated which ARG subtypes contributed to the identified peptide resistance category and found that ~90% of the expressed peptide resistance was directly contributed by a single resistance gene, YojI, which was found to be widely distributed among the major taxa comprising the BWWTP community such as the Comamonadaceae (Figure 1—figure supplement 1). YojI is typically associated with resistance to microcins by reducing the intracellular concentration of the toxic antibiotic peptide (Delgado et al., 2005). The high incidence of this gene indicates a broad adaptive strategy amongst the microbial populations in the BWWTP against these specific stressors.

Microbial community and co-occurrence patterns of AMR

Based on the previously identified microbial community (Herold et al., 2020), we hypothesized that the abundant and prevalent bacterial orders such as Acidimicrobiales were major contributors to the abundance in ARGs observable via metagenomics. To further investigate the contribution to AMR by the distinct microbial populations, we linked AMR genes to the contig-based taxonomic annotations of the assemblies (Methods). Herein, we identified a wide variety of taxonomic orders contributing to AMR, with multiple orders often contributing to the same resistance categories (Figure 2—figure supplement 1). Overall, taxa belonging to Acidimicrobiales, followed by Burkholderiales, were found to encode most of the ARGs (Figure 2a). Additionally, the abundance of ARGs linked to taxonomy varied over time. This was most noticeable during a five-week period (autumn: 02-11-2011 to 29-11-2011), where a decrease in abundance in ARGs linked to Acidimicrobiales and Bacteroidales was observed coinciding with an increase in ARG abundance in Pseudomonadales and Lactobacillales.

Figure 2 with 2 supplements see all
Download asset Open asset
Microbial population-linked AMR.

(a) Longitudinal ARG relative abundance levels linked to their corresponding microbial taxa (order level). Colors correspond to AMR categories. (b) Relative abundance of AMR categories linked to Candidatus Microthrix parvicella with colors corresponding to AMR categories. (c) Association network depicting co-occurrence patterns of individual antimicrobial resistance genes (ARGs) and microbial taxa on genus level. Nodes represent taxa or ARG with the node size representing the number of edges. The size of the edges represent the strength of interaction between the nodes.

Since the order Acidimicrobiales was found to be linked to the highest abundance in ARGs, we further resolved the taxonomic affiliation and identified the species Candidatus Microthrix parvicella (hereafter known as M. parvicella) to be the main contributor to AMR. M. parvicella was previously found to dominate this microbial community (Martínez Arbas et al., 2021) and is a well-characterized bacterium commonly occurring in the BWWTP (Calusinska et al., 2018). Overall, aminoglycoside, beta-lactam, multidrug, and peptide resistance were found to be abundant in this species (Figure 2b), with aminoglycoside resistance demonstrating the highest expression levels as confirmed through metatranscriptomic analysis (Figure 2—figure supplement 2). Although it was not surprising to find a high abundance of ARGs linked to this species, the longitudinal variation in the abundances of these ARGs was nevertheless surprising (Figure 2b). Furthermore, coupled to a decrease in the abundance of M. parvicella itself (Martínez Arbas et al., 2021), we observed an almost complete decrease in ARGs at two timepoints (23-11-2011 and 29-11-2011). However, the M. parvicella population recovered to levels resembling the earlier timepoints in conjunction with the abundances in ARGs toward the end of the sampling period (Figure 2a, Figure 2b), underlining their overall contribution to AMR within this BWWTP. Alternatively, it is plausible that the dominance of M. parvicella is attributable to the encoded ARGs, which in turn, may confer a fitness advantage.

In order to determine whether the abundances in ARGs may be directly associated with the community composition over time, co-occurrence patterns between ARG subtypes and taxa (genus level) were explored using the metagenomic data. Association network analyses (Figure 2c) demonstrated that ARGs, within or across ARG types and microbial taxa, showed clear and distinct co-occurrence patterns within the BWWTP. These patterns indicated a strong segregation of distinct, taxa-specific ARG subtypes within the BWWTP community over time. One clear example was that of M. parvicella which encoded different aminoglycoside resistance genes (Figure 2c). Thus, the abundance of this bacterium along with the aminoglycoside ARGs were highly correlated.

Monitoring pathogenic microorganisms within BWWTPs

In conjunction with the families observed within BWWTPs, we also found that certain ESKAPEE pathogens (Reza et al., 2019), such as Klebsiella spp. and Pseudomonas spp., demonstrated co-occurring patterns with ARGs (Figure 2c).

As previously mentioned, BWWTPs represent a collection of potentially pathogenic microorganisms originating from, among others, the human population. Moreover, evidence suggests that ARGs from environmental and commensal bacteria can spread to pathogenic bacteria through HGT (MacLean and San Millan, 2019). Therefore, we assessed the presence of AMR in the extended priority list of pathogens (Table 1), characterized as such by the WHO (Tacconelli et al., 2018), using both metagenomics and metatranscriptomics.

Table 1
WHO priority list for research and development of new antibiotics for antibiotics-resistant bacteria (Tacconelli et al., 2018).
BacteriaPriorityOrganism detectedResistance detected
Acinetobacter baumanniiCritical++
Pseudomonas aeruginosaCritical++
EnterobacteriaceaeCritical++
Enterococcus faeciumHigh++
Staphylococcus aureusHigh++
Helicobacter pyloriHigh++
Campylobacter sppHigh+-
Salmonella sppHigh++
Neisseria gonorrhoeaeHigh+-
Streptococcus pneumoniaeMedium++
Haemophilus influenzaeMedium+-
Shigella sppMedium++

Of the identified pathogens (Table 1), we found that Pseudomonas aeruginosa, both encoded and expressed the highest abundance of ARGs, followed by Acinetobacter baumannii, over time within the BWWTP (Figure 3). Moreover, an increase in ARG abundance and expression was observed in Pseudomonas aeruginosa during the time period, during which the otherwise dominant M. parvicella demonstrated reduced abundance (Figure 2b and Figure 3).

Figure 3
Download asset Open asset
Assessment of AMR associated with clinical pathogens.

ARG relative abundances encoded and expressed by clinical pathogens over time within the BWWTP, with colors corresponding to the identified pathogens.

Differential transmission of antimicrobial resistance via mobile genetic elements

As previously described (Beceiro et al., 2013; Wee et al., 2020), the mobilome is a major contributor to the dissemination of AMR within a microbial community. Consequently, to understand (i) the role of MGE-mediated AMR transfer within the BWWTP, and (ii) to identify differential contribution of the mobilome to the dissemination of AMR, we identified both plasmids and phages within the metagenome and linked these to the respective ARGs. While, as expected, the majority of ARGs was found to be encoded on the bacterial chromosome (Figure 4—figure supplement 1), we also found that plasmids contributed to an average of 10.8% of all ARGs, while phage contributed to an average of 6.8% of all resistance genes, in agreement with the general hypothesis that conjugation has the greatest influence on the dissemination of ARGs (adj.p<0.05, One-way ANOVA) (MacLean and San Millan, 2019). This phenomenon, however, varied across time within the BWWTP (Figure 4a).

Figure 4 with 2 supplements see all
Download asset Open asset
MGE-derived AMR within the BWWTP resistome.

(a) Overall relative abundance of MGEs encoding ARGs. Contribution of plasmids to AMR (average of 10.8% of all ARGs) was found significantly increased compared to phages (average of 6.8% of all ARGs) (adj.p <0.05, One-way ANOVA). Colors depict the different MGE predictions (phage, plasmid, ambiguous) (b) Boxplots depicting significant (adj.p < 0.05, Two-way ANOVA) differential abundances of ARGs encoded by plasmids (blue) vs phages (green). (c) Relative abundance of the six significantly different AMR categories encoded on phages over time, with colors corresponding to AMR categories. (d) Relative abundance of the six significantly different AMR categories encoded on plasmids over time, with colors corresponding to AMR categories.

When investigating the dissemination of AMR via MGEs, most reports typically focus on either phages or plasmids individually, or both as collective contributors to transmission (Slizovskiy et al., 2020). To date and to our knowledge, the respective contributions of phages and plasmids to AMR transmission have not been subjected to a comprehensive comparative analysis. To facilitate a systematic, comparative view of MGE-mediated AMR, we assessed the segregation of MGEs with respect to AMR and found that phages and plasmids contributed differentially to AMR (Figure 4—figure supplement 2). Specifically, we tested 28 AMR categories with respect to their association with MGEs and found a significant difference in six AMR categories when comparing ARGs encoded by phages and plasmids (adj.p <0.05, Two-way ANOVA)(Figure 4b). Aminoglycoside, bacitracin, MLS (i.e. macrolide, lincosamide, and streptogramin) and sulfonamide resistance were found to be primarily encoded by plasmids, whereas fosfomycin and peptide resistance were found to be associated with phages.

To further understand AMR in relation to the community dynamics, we investigated the abundance and segregation of the above-mentioned significant resistance categories at different timepoints within the BWWTP. We observed ARG abundances varied over time both in phages (Figure 4c) as well as plasmids (Figure 4d). For instance, the abundance in aminoglycoside and sulfonamide resistance, which was encoded primarily by plasmids (Figure 5—figure supplement 1a), fluctuated widely over time in both phages and plasmids (Figure 4c). Additionally, plasmid-mediated sulfonamide resistance was reduced at 23-11-2011, followed by its highest abundance a week later (29-11-2011), while subsequently again decreasing. Similarly, in line with the above observations, fosfomycin and peptide resistance genes, while segregating within phages, demonstrated significant fluctuations over time (Figure 4d). In addition to the metagenome, we also contextualized the localization of the expressed ARGs within MGEs based on the metatranscriptomic information. Specifically, we found that plasmids demonstrated a significantly increased expression of aminoglycoside along with bacitracin and sulfonamide resistance genes, while the expression of glycopeptide, mupirocin and peptide resistance genes were primarily enriched in phages (Figure 5a). These observations pertaining to plasmid-mediated AMR were in line with the metagenomic findings (Figure 4b). Only peptide resistance was observed to be expressed via phages in contrast to the differential enrichment of fosfomycin resistance observable in the metagenomic data.

Figure 5 with 2 supplements see all
Download asset Open asset
Taxonomic affiliations of MGE-derived resistance genes.

(a) Boxplot depicting significant differential abundance (n=51 per group, adj.p <0.05, Two-way ANOVA) of ARGs expressed in plasmids vs phages. (b) Tripartite network assessing the association of MGE-derived ARGs with the microbial taxa. Thickness of the lines representing potential niche-partitioning of the AMR category to one MGE over the other. Color of the line representing which MGE the AMR is linked to: green (phage), blue (plasmid) or black (both phage and plasmid). Asterisk denominates taxonomic orders which include known clinical pathogens. (c) Alluvial plot depicting the mean abundance (log10) of MGE-derived ARGs encoded (metagenome) and/or expressed (metagenome) by clinical pathogens. Colors of all panels correspond to the MGEs and AMR categories.

Taxonomic affiliations of MGE-derived resistance genes

When assessing the differential contributions of MGEs to AMR, we found congruency between plasmids and phages to the AMR categories and taxonomic affiliations (Figure 5b). For example, in the metagenomic data MGEs (phage and plasmid) were predominantly associated with the same AMR category and subsequently the same taxa. However, some exceptions were observed with specific taxa associated with AMR either through plasmids or phages. For instance, MLS resistance in Bacteroidales and Nostocales was mediated solely through plasmids, whereas the same resistance category was mediated by phage in Bifidobacteriales, indicating segregation of AMR between taxa and MGEs.

As most bacteria harbor MGEs, we queried whether the MGE-mediated AMR categories were linked to the abundance of some of the earlier reported taxa. Interestingly, we found that peptide resistance encoded by M. parvicella was solely associated with phages, while aminoglycoside resistance was primarily correlated with plasmids (Figure 5—figure supplement 1b). Other highly abundant taxa such as Pseudomonas and Comamonas (Figure 5—figure supplement 1c-d), on the other hand, were correlated with sulfonamide resistance in addition to aminoglycoside resistance encoded on plasmids (Figure 5b). This was further reflected within the metatranscriptome data where in taxa such as Acidimicrobiales the expression levels of aminoglycoside resistance were solely associated with plasmids (Figure 5—figure supplement 2a). Additionally, in the Burkholderiales family, peptide resistance was found to be expressed through phages (Figure 5—figure supplement 2b).

We also found a clear segregation of the mobilome with respect to individual pathogens in the metagenome. Interestingly, plasmids were exclusively associated with AMR in six out of the fourteen relevant taxa (Figure 5c). These included Streptococcus pneumoniae, Staphylococcus aureus, Shigella flexneri, Klebsiella pneumoniae, Enterobacter kobei, and Enterobacter hormaechei. Furthermore, the plasmids were also associated with conferring peptide, multidrug, MLS, beta-lactam, fluoroquinolone, bacitracin, aminoglycoside, aminoglycoside:aminocumarin and sulfonamide resistance. Phages were exclusively associated with glycopeptide and aminoglycoside resistance in Salmonella enterica. Overall, our results revealed for the first time the key segregation patterns of AMR via the mobilome in taxa that are of relevance to human health and disease. Moreover, substantiating the metagenomic data, the pathogenic bacteria S. pneumoniae, S. aureus, K. pneumoniae, E. kobei, and E. hormaeche were found to express ARGs solely associated with plasmids (Figure 5c). Collectively, these findings indicate an imminent threat to global health due to the potential dissemination of resistant pathogens across reservoirs.

Metaproteomic validation of AMR abundance and expression

In order to expand our findings with the expression (metatranscriptomic) analyses on the BWWTP, we further used the corresponding metaproteomic data to offer complementary information at the protein level. Similar to the metagenome data we found protein expression linked to aminoglycoside, beta-lactam and multidrug resistance, over time within the BWWTP (Figure 6—figure supplement 1). Proteins, especially those linked to multidrug resistance were found to increase over time.

To further unravel AMR expression and assess its stability across time, we estimated the normalized protein index (NPI) per gene, as discussed in the Methods, by integrating all of the multi-omic data. The estimated NPI demonstrated stable levels of aminoglycoside and multidrug resistance within the BWWTP (Figure 6a). Specifically, proteins conferring multidrug resistance were found to increase over time, which is in line with the gene- and expression-level observations. Furthermore, we contextualized the normalized proteins conferring AMR to their localization on MGEs. We identified five resistance categories, that is aminoglycoside, beta-lactam, sulfonamide, multidrug, and tetracycline resistance, to be expressed through MGEs (Figure 6b). Of these categories we found that aminoglycoside resistance, in concordance with the gene and expression levels, was more strongly associated with plasmids than phages (adj.p <0.05, Two-way ANOVA). We further found that the MGE-mediated AMR categories were associated with specific microbial taxa. with plasmid-mediated aminoglycoside resistance found to be strongly associated with the previously mentioned M. parvicella (Figure 6b). On the other hand, we did not identify any peptides associated with the ESKAPEE pathogens via metaproteomics.

Figure 6 with 1 supplement see all
Download asset Open asset
Integrative multi-omic assessment of AMR.

(a) Metagenomic and metatranscriptomic normalized protein levels linked to AMR within the WTTP over time. (b) Tripartite network assessing the normalized protein levels derived from MGEs and associated taxa. Boxplots depicting significant differential (n=51 per group, adj.p <0.05, Two-way ANOVA) abundance of aminoglycoside resistance in plasmid versus phage in Candidatus Microthrix parvicella as well as overall. Colors of all panels correspond to the MGEs and AMR categories.

Discussion

The surveillance of wastewaters for the identification of microbial molecular factors is a critical tool for identifying potential pathogens. This has been highlighted recently with the tracking of SARS-CoV-2 within wastewater treatment plants to assess viral prevalence and load within a given community (Westhaus et al., 2021). Such approaches have also been employed for screening for antimicrobial resistance at a population level (Kwak et al., 2015; Reinthaler et al., 2013). So far, several studies (Calero-Cáceres et al., 2014; Hendriksen et al., 2019; Parsley et al., 2010; Szczepanowski et al., 2009) have characterized the proliferation of ARGs and antibiotic resistant bacteria in BWWTP. Szczepanowski et al., 2009 identified 140 clinically relevant plasmid-derived ARGs in a BWWTP metagenome while Parsley et al., 2010 characterized ARGs from bacterial chromosomes, plasmids and in viral metagenomes found in a BWWTP. Further studies have shown that conventional BWWTP processes at best only partially remove ARGs from the effluent and may find their way into the urban water cycle (Hiller et al., 2019; Proia et al., 2018; Rodriguez-Mozaz et al., 2015). Wastewater treatment plants, therefore, are crucial reservoirs of AMR, whose monitoring may allow for early-detection of AMR within the human population feeding into the system. Here, we leveraged a systematic and longitudinal sampling scheme from a BWWTP to identify diverse AMR categories prevalent within the BWWTP microbial community. In line with the studies by Szczepanowski et al., 2009 and Parsley et al., 2010, we found up to 29 AMR categories with several ARGs within the BWWTP. More importantly, and unlike the previous studies, we linked the identified ARGs to clinically-relevant ESKAPEE pathogens, which represent a growing global threat to human health.

In our BWWTP samples, we identified a core group of 15 AMR categories that were ubiquitous at all timepoints. In line with the above-mentioned reports, the observed core resistance categories may reflect their abundance in the surrounding human population (Aarestrup and Woolhouse, 2020). This has previously been reported by Pärnänen et al., 2019, Su et al., 2017 and Hendriksen et al., 2019 where they showed that BWWTP AMR profiles correlate with clinical antibiotic usage as well as other socio-economic and environmental factors. Furthermore, bacteria are known to have innate defense mechanisms against inhibitory secondary metabolites from other taxa (Frost et al., 2018). Therefore, one must be cognizant of thprace phenomenon that the observed core group of AMR categories may also be a proxy for the abundance of specific resistant bacteria. Despite this observation, it is plausible that both anthropogenic and microbial sources for AMR play a role in the observed resistance categories within the BWWTP. For instance, we found that several AMR categories, including ancillary (prevalent, moderate, and rare) groups, were associated with M. parvicella within the BWWTP. While the mechanisms underlying the abundance of this taxon are numerous, the increased abundance of ARGs are likely to contribute to its fitness advantage in this metal- and antibiotic-rich environment (Baquero et al., 2008). Additionally, similar to the findings by Munck et al., 2015, we found a wide range of bacteria associated with AMR categories including Acidimicrobiales, Burkholderiales and Rhodocyclales. On the other hand, we report that taxa, including ESKAPEE pathogens, belonging to 25 bacterial orders were associated with 29 AMR categories, compared to the eight bacterial orders reported previously.

It is important to note that the mobilome plays a critical role in the dissemination of AMR within microbial communities. AMR from resistant bacteria within the BWWTP can quickly disseminate within the BWWTP (Che et al., 2019; Fouz et al., 2020), including transmission from pathogenic to commensal species (Blake et al., 2003; Brinkac et al., 2017). As a result, mediated through HGT, the BWWTP becomes a hotspot for resistant bacteria, which are then released back into the receiving environment (Turolla et al., 2018), and eventually the human population (Fouz et al., 2020; Newton and McClary, 2019). Therefore, to limit the dissemination of AMR, it is important to understand the role of MGEs within the BWWTP. Our comprehensive analyses identified the differential contributions of AMR transmission mediated via phage and plasmid (Figure 7). Specifically, we identified clear segregation of aminoglycoside, bacitracin, MLS and sulfonamide resistance categories with plasmids, while fosfomycin and peptide resistance were increasingly encoded and conferred via phages. While the association between these AMR categories and plasmids (Dubnau et al., 1981; Galimand et al., 2003; Han et al., 2015; Razavi et al., 2017) or phages (Torres-Barceló, 2018) are in line with previously reported results, differential analysis between MGEs has not been previously reported and has not been performed on multi-omic levels. In this study we report for the first time the systematic and extensive comparison of AMR encoded and expressed by phages versus plasmids. Our results indicating the segregation of ARGs within the ESKAPEE taxa via the MGEs further provide insights into potential modes of AMR transmission among pathogens. Although one cannot exclude the possibility of transmission of the above-mentioned ARGs via other MGEs, identifying potential segregation of MGEs in the transmission of ARGs brings us one step closer to identifying specific transmission paths and limiting the spread of AMR. For example, some studies have reported plasmid ‘curing’, the process by which plasmids are lost from bacterial populations, as a strategy against dissemination of AMR (Bouanchaud and Chabbert, 1971; Vrancianu et al., 2020). As described by Buckner et al., 2018 plasmid curing, as well as other anti-plasmid strategies, could both reduce AMR prevalence, and (re-)sensitize bacteria to antibiotics (Buckner et al., 2018).

Figure 7
Download asset Open asset
Separation of MGE-derived AMR within the BWWTP.

A graphical summary highlighting AMR categories found significantly increased in phage versus plasmid in all three omes.

Although several methods and tools exist for the identification of MGEs, the linkage to their respective bacterial hosts still remains a challenge. Our method presented here involves taxonomic annotations of reconstructed genomic information as a mean of linking bacteria to their corresponding MGEs using stringent filtering/identification criteria. An alternative approach may involve targeted sequencing of plasmids of interest and their respective hosts using the methodology established by Li et al., 2018. However, to use such targeted sequencing approaches one requires a priori information on which plasmids to focus on in contrast to our method which is agnostic to such prior information. Furthermore, Hi-C-based methods (Lieberman-Aiden et al., 2009; Press et al., 2017) capture inter-chromosomal junctions, such as plasmid-genome interactions. However, compared to metagenomic approaches, these methods require different and extensive processing of the samples prior to sequencing to induce covalent linkages among DNA molecules (Lieberman-Aiden et al., 2009; Press et al., 2017), therefore precluding the use of this method for large-scale characterization efforts. Combining these strategies in future studies with AMR categorization according to taxonomic affiliation as well as linkage to specific MGEs will provide novel strategies for characterizing and subsequently affecting MGE-mediated AMR.

By complementing the metagenomic analyses, metatranscriptomics conferred essential information regarding gene expression within the resistome. For instance, when comparing AMR expression levels of aminoglycoside, bacitracin, and sulfonamide mediated via MGEs, it is noticeable that expression levels in plasmids mirror the genomic content, that isthey exhibited higher levels of expression when compared to phage. On the other hand, glycopeptide and mupirocin resistance genes which were highly expressed in phages were not reflected within the metagenomic data. Additionally, we found the YojI resistance gene to be more highly expressed than any other ARGs. To facilitate resistance against the peptide antibiotic microcin J25, the outer membrane protein, TolC, in combination with YojI is required to export the antibiotic out of the cell (Delgado et al., 2005). Microcin J25 belongs to the group of ribosomally synthesized and post-translationally modified peptides (RiPPs) and has antimicrobial activity against pathogenic genera such as Salmonella spp. and Shigella spp. (Naimi et al., 2018). Interestingly, it has only recently been proposed as a treatment option against Salmonella enterica and has been discussed in recent years as a potential novel antibiotic (Ben Said et al., 2020). Based on these results, by considering that BWWTPs may reflect both the presence of AMR within the human population as well as be a hotspot of dissemination and generation of new AMR, surveillance of BWWTPs must be emphasized when developing new antibiotics. Our findings collectively suggest that the differential capacity of MGEs to disseminate AMR, coupled with longitudinal and expression-level analyses are crucial for monitoring human health conditions. More importantly, we report for the first time that BWWTP monitoring for AMR may allow for early detection of resistance mechanisms previously undescribed in BWWTP.

Finally, we applied an integrated multi-omic analysis approach to improve our knowledge on the functional potential of AMR and simultaneously validate the abundance and expression findings of the ARGs. By normalizing the metaproteomic results with the normalized expression of genes, we were able to assess the stability of expressed AMR across time. We find that our methodology allows for an unbiased assessment of overall expression accounting for gene copy abundance and expression. These findings support the notion that the AMR genes may serve as sentinels or indicators of the presence of particular antimicrobial agents. However, it is plausible that we are only identifying the most abundant proteins and/or proteins that are more stable over time, and do not capture the entirety of the proteome profiles. Factors such as protein decay rates (Cameron and Collins, 2014) among others, may additionally influence this assessment. Irrespective of these observations, we identified segregation of AMR categories with respect to plasmids and phages.

Our findings also highlight the potential for identifying segregation of AMR via specific MGEs with an aim toward possible therapeutic and mitigation strategies via for example plasmid curing. Furthermore, we demonstrate that longitudinal analyses are required to survey AMR within BWWTPs due to the variations in the resistome across time. These shifts may either be representative of a shift within the human population itself, which in turn could be associated with the concurrent use of antibiotics at a given time, or competition within the microbial community. In any case, an independent or static analysis of the various time points may show an incomplete view of the BWWTP resistome, thus underlining the importance of our longitudinal resistome analyses. Overall, our findings suggest that BWWTPs are critical reservoirs of AMR, potentially allowing for early detection and monitoring of pathogens. In addition, BWWTP monitoring may allow detection of resistance mechanisms linked to the introduction of new antimicrobials. Finally, BWWTPs may serve as a model for understanding the segregation of MGEs through AMR.

Methods

Sampling and biomolecular extraction

From within the anoxic tank of the Schifflange municipal biological wastewater treatment plant (located in Esch-sur-Alzette, Luxembourg; 49° 30′ 48.29″ N; 6° 1′ 4.53″ E) individual floating sludge islets were sampled according to previous described protocols (Herold et al., 2020). Sampling was performed starting on 21-03-21 till 03-05-2012 in approximately one-week intervals resulting in a total of 51 samples. DNA, RNA and proteins were extracted from the samples in a sequential co-isolation procedure as previously described (Roume et al., 2013).

Sequencing and data processing for metagenomics and metatranscriptomics

Paired-end libraries were generated for metagenomics with the AMPure XP/Size Select Buffer Protocol following a size selection step recommended by the standard protocol. Libraries for metatranscriptomics were prepared from RNA after washing stored extractions with ethanol and depletion of rRNAs with the Ribo-Zero Meta-Bacteria rRNA Removal Kit (Epicenter). Subsequently, the ScriptSeq v2 RNA-seq library preparation kit (Epicenter) was used for cDNA library preparation, followed by sequencing on an Illumina Genome Analyses IIx instrument with 100-bps paired-end protocol. Processing and iterative co-assembly of metagenomic and metatranscriptomic reads was done using the Integrated Meta-omic Pipeline (Narayanasamy et al., 2016) (IMP v1.3; available at https://r3lab.uni.lu/web/imp/). For read processing, Illumina Truseq2 adapters were trimmed, and reads of human origin were filtered out, followed by a de novo assembly with MEGAHIT (Li et al., 2015) v1.0.6. For the assembly, both metagenomic and metatranscriptomic reads were co-assembled to increase contiguity of the assemblies and to improve read usage. Additional information regarding the read coverage and depth for each sample is available at https://git-r3lab.uni.lu/laura.denies/lao_scripts.

Identification of antimicrobial resistance genes, mobile genetic elements and taxonomy

The assembled contigs from IMP were used as input for PathoFact (de Nies et al., 2021), for the prediction of antimicrobial resistance genes and MGEs. ARGs were further collapsed into their respective AMR categories, as identified by PathoFact in accordance with those provided by the Comprehensive Antibiotic Resistance Database (CARD) (Alcock et al., 2020). Furthermore, utilizing PathoFact, AMR genes were linked to predicted MGEs (i.e. plasmids and phages) to track transmission of AMR. By considering all different predictions of MGEs, a final classification was made by PathoFact based on the genomic contexts of the AMR genes encoded on plasmids, phages or the organismal chromosomes, including the classification of those that could not be resolved (ambiguous). The AMR genes that could not be assigned to either the MGEs or bacterial chromosomes were subsequently referred to as unclassified genomic elements.

Identified ARGs and their categories were further linked to associated microbial taxonomies using the taxonomic classification system Kraken2 (Wood et al., 2019). Kraken2 was run on the contigs using the maxikraken2_1903_140 GB (March 2019, 140GB) (https://lomanlab.github.io/mockcommunity/mc_databases.html) database (Wood et al., 2019). To ensure confidence in taxonomic classification, stringent cutoffs based on a minimum 70% identity of the contig with the database, across 90% of the contig length, was used (Wood et al., 2019). Since the same contig-based assembly file was used as input for PathoFact as well as for Kraken2, the prediction of antimicrobial resistance genes and mobile genetic elements as well as functional- and taxonomic- annotations were linked based on the contig they were encoded on (Figure 8).

Figure 8
Download asset Open asset
Identification of ARGs and contextualization of MGEs in relation to taxa.

A schematic diagram depicting the bioinformatic workflow to identify ARGs and the subsequent contextualization of MGEs in relation to microbial taxa. The green features represent the different tools and pipelines used while the red features highlight the data used and generated in this process.

Though several methods and tools exist for the identification of plasmids and/or phages, tools identifying the respective bacterial hosts for MGEs are rare. The hosts for phages can be determined based on classification and linking of CRISPR and spacer sequences between bacterial and phage sequences respectively (Bland et al., 2007; Martínez Arbas et al., 2021). Plasmid hosts, on the other hand, due to the competence of plasmids (including self-replication, etc.), are difficult to identify (Suzuki et al., 2010) or can only be classified against existing databases (Aytan-Aktug et al., 2022).

Therefore, to streamline this process, we used taxonomic associations as a means of linking bacteria with their corresponding MGEs using the stringent filtering/identification criteria. Based on the already established stringent cutoffs for plasmid/phage predictions, the original MGE predictions from PathoFact were used without any additional processing and the taxonomic classifications from Kraken2 were used to assign putative hosts. Finally, to obtain the gene copy number and transcriptome expression levels, both metagenomic and metatranscriptomic reads were independently mapped to the corresponding assembly files per sample. The raw read counts per contig as well as per individual ORF, as given by PathoFact, were determined with the featureCounts option. The relative abundance of the ARGs was calculated using the RNum_Gi method previously described by Hu et al., 2013.

Metaproteomics and data analyses

Raw mass spectrometry files were converted to MGF format using MSconvert (Chambers et al., 2012) with default parameters. The metaproteomic searches were performed using SearchGUI / PeptideShaker (Vaudel et al., 2015) for each time point. To generate the databases, each predicted protein sequence file was concatenated with the cRAP database of contaminants (common Repository of Adventitious Proteins, v 2012.01.01; The Global Proteome Machine) and with the human UniProtKB Reference Proteome (Consortium, 2021). In addition, inverted sequences of all protein entries were concatenated to the databases for the estimation of false discovery rates (FDRs). The search was performed using SearchGUI-3.3.20 (Barsnes and Vaudel, 2018) with the X!Tandem (Langella et al., 2017), MS-GF+ (Kim and Pevzner, 2014) and Comet (Eng et al., 2013) search engines using the following parameters: Trypsin was used as the digestion enzyme and a maximum of two missed cleavage sites was allowed. The tolerance levels for identification were 10 ppm for MS1 and 15ppm for MS2. Carbamidomethylation of cysteine residues was set as a fixed modification and oxidation of methionines was allowed as variable modification. Peptides with a length between 7 and 60 amino acids and with a charge state composed between +2 and+4 were considered for identification. The results from SearchGUI were merged using PeptideShaker-1.16.45 (Vaudel et al., 2015) and all identifications were filtered in order to achieve a peptide and protein FDR of 1%.

Each predicted protein sequence corresponded to the predicted ORFs generated by the Prodigal (version 2.6.3) predictions included in PathoFact. As such predicted protein sequences matched the ARG annotation of the ORFs as provided by PathoFact.

Multi-omic data integration

To further improve upon the understanding of the AMR expression and assess its stability across time, we estimated the normalized protein index (NPI) per gene, by integrating the multi-omic data. To estimate the NPI, we first normalized the metaT abundance based on per gene copy numbers obtained via the metagenomic abundance:

NPI =NmetaproteomeNmetatranscriptome / Nmetagenome

This, the normalized expression of genes, yields the per copy expression of ARGs within each AMR category. Subsequently, the normalized expression was used to standardize the metaP abundances for those genes where the necessary data was available.

MGE partition assessment

To assess the segregation of MGEs through AMR we determined niche regions and overlap using the nicheROVER R package (Swanson et al., 2015). nicheROVER uses Bayesian methods to calculate niche regions and pairwise niche overlap using multidimensional niche indicator data (i.e. stable isotopes, environmental variables). As such, using AMR as the indicator data, we extended the application of nicheRover to calculate the probability for the size of the niche area of one MGE inside that of the other, and vice versa. We calculated the segregation size estimate for each MGE and additionally generated the posterior distributions of μ (population mean) for each AMR category in all omics. We further computed the niche overlap estimates between MGEs with a 95% confidence interval over 10,000 iterations.

Data analysis

Figures for the study including visualizations derived from the taxonomic and functional analyses were created using version 3.6 of the R statistical software package (R Development Core Team, 2013). A paired two-way ANOVA (Analysis of Variance, accounting for MGEs and AMR, or taxa and AMR) within the nlme package was used for identifying statistically significant differences. Tripartite and the association networks were generated using the SpiecEasi (Kurtz et al., 2015) R package where a weighted adjacency matrix was generated using the Meinhausen and Buhlmann (mb) algorithm, with a nlambda of 40, and lambda minimum ratio at 0.001. The analyses were bootsrapped with n=999 to avoid overfitting, autocorrelations and false network associations. The network was further refined, selecting for positive edges, with a degree greater than the mean-degree of the initial network. Additionally, edges with a correlation value greater than 0.6 were retained where the associated significance was below an adjusted p-value <0.05. The association networks were based on two types of edges, that is, co-occurences between taxa based on ARG abundances and ARG-taxon edges based on ARGs encoded on specific taxa. The igraph (Csardi and Nepusz, 2006) package was used in R to render the graphics for the network. All code for visualization and analysis is available at: https://git-r3lab.uni.lu/laura.denies/lao_scripts.

Data availability

The genomic FASTQ files used in this work (previously published) are publicly available at NCBI BioProject PRJNA230567. Metaproteomic data (previously published) are publicly available at the PRIDE database under accession number PXD013655. The open-source tools and algorithms used for the data analyses are reported in the Methods section, including relevant flags used for the various tools. Additionally, custom code for further analysis and generation of the figures can be found at: https://git-r3lab.uni.lu/laura.denies/lao_scripts.

The following previously published data sets were used

References

  11. Book
    1. Bonilla AR
    2. Muniz KP.
    (2009)
    Antibiotic Resistance: Causes and Risk Factors, Mechanisms and Alternatives
    Nova Science Publishers.
    1. Csardi G
    2. Nepusz T
    (2006)
    The igraph software package for complex network research
    InterJournal, Complex Systems 1695:1–9.
  29. Book
    1. Dubnau D
    2. Grandi G
    3. Grandi R
    4. Gryczan TJ
    5. Hahn J
    6. Kozloff Y
    7. Shivakumar AG.
    (1981)
    Regulation of Plasmid Specified MLS-Resistance in Bacillus subtilis by Conformational Alteration of RNA Structure In
    In: Levy SB, Clowes RC, Koenig EL, editors. Molecular Biology, Pathogenicity, and Ecology of Bacterial Plasmids. Boston, MA: Springer US. pp. 157–167.
  57. Book
    1. O’Neill J
    (2014)
    The Review on Antimicrobial Resistance
    London, United Kingdom: wellcome trust.
  63. Software
    1. R Development Core Team
    (2013) R: a language and environment for statistical computing
    R Foundation for Statistical Computing, Vienna, Austria.
    http://www.r-project.org
  80. Conference
    1. Varrette S
    2. Bouvry P
    3. Cartiaux H
    4. Georgatos F
    (2014) Management of an academic HPC cluster
    The UL experience2014 International Conference on High Performance Computing Simulation (HPCS. pp. 959–967.
    https://doi.org/10.1109/HPCSim.2014.6903792

Decision letter

  1. María Mercedes Zambrano
    Reviewing Editor; CorpoGen, Colombia
  2. Bavesh D Kana
    Senior Editor; University of the Witwatersrand, South Africa
  3. María Mercedes Zambrano
    Reviewer; CorpoGen, Colombia

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Decision letter after peer review:

[Editors’ note: the authors submitted for reconsideration following the decision after peer review. What follows is the decision letter after the first round of review.]

Thank you for submitting the paper "Mobilome-driven segregation of the resistome in biological wastewater treatment" for consideration by eLife. Your article has been reviewed by 2 peer reviewers, including María Mercedes Zambrano as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by a Senior Editor.

Comments to the Authors:

We are sorry to say that, after consultation with the reviewers, we have decided that this work will not be considered further for publication by eLife.

Specifically, more methodological details and rigorous statistical analyses are needed in several places, as well as more thorough and clearer information to fully assess the validity of the results and the subsequent claims.

Reviewer #1 (Recommendations for the authors):

This work aims to provide a broad view of antimicrobial resistance determinants in a wastewater treatment plant by conducting a multi-omics approach on samples taken over 1.5 years. The use of several omics approaches and the distinction of markers associated with either phages or plasmids in the system represent strengths of the study. However, the lack of clarity and methodological details weaken the interpretation and impact of the results.

1. While the multiple omics approach provides valuable insight, it is not entirely clear how expression data adds to what is already known regarding ARGs in the environment and how this expression information sheds new light on the resistome or can be "crucial for monitoring human health conditions" (L412). In this respect it would have been good to see if their findings on ARGs segregating between phages and plasmids correlate with what is known about resistances, at least in the relevant taxa (eg, L 288-9).

2. Since the results are based on assembled data (L 158), the authors should provide additional and sufficient information regarding these assemblies and how they can be confident of taxa identification to the species level (such as the ESKAPEE group) and the linkage of ARGs to particular microorganisms. This is important, given that much of the analysis and ensuing conclusions are based on linkage of ARGs and MGEs to particular taxa. It would also be important to know how they managed to precisely map expression data to these MGEs and not to chromosomal loci and if perhaps this could also explain some of the discrepancies observed (eg, L397-8).

3. It seems that about 18% of all ARGs are associated with MGEs, yet no additional mention is made of the remaining genes, which presumably are chromosomal. Please clarify this in the document and also indicate how this could affect the transcriptomic and proteomic analysis.

4. It is intriguing that they observe high levels of peptide resistance in their metatranscriptome (L146) associated with resistance in E. coli, which is not a major contributor in the community and the gene was apparently not abundant in their metagenomic analysis. This discrepancy is not fully addressed and is relevant because it can also reflect biases in the approach or conditions of the sites being analyzed.

5. It seems that the majority of markers are associated with abundant taxa, which is to be expected, but it is not clear if this is the reason why the authors hypothesize that the majority of ARGs would be in the major taxon M. parvicella (L169 and L175) or why these ARGs would confer fitness advantage (L184), unless of course this correlates with antibiotics present, which was not determined.

6. The document can also benefit from revision of the text, particularly the discussion, to clarify some claims and make the document easier to follow. There are some statements that do not accurately reflect either the scope of the work or the results obtained. Some of these include:

L79 – why would ARGs and sub-inhibitory antibiotic concentrations (which should not be lethal) promote HGT?

L 96 – They do not really shed light on "the evolution" and at best suggest dissemination differences.

L102 – The analysis suggests but does not necessarily "demonstrate" that MGEs are important drivers of AMR dissemination. In addition, it does not necessarily follow that "assessing the activity of the ARGs" can lead to understanding the underlying mechanisms.

L186 – why population sizes? Did you actually look at association with the size?

L192 – Is this figure 3a?

On L202 it says that they "assessed the acquisition and dissemination of AMR". The data shows only associations, which is not necessarily acquisition and at best hints at modes of dissemination.

L292-94- the findings indicate but are not themselves a threat or have potential for dissemination. Please correct.

L296 – it does not seem that the proteomics validated transcriptomics data, but rather are more of a complementary analysis (L298). Likewise, on L416- It does not seem that proteomics data validated the abundance and expression findings.

L413 and 438 – Monitoring based on sequence analysis is based on comparisons using known data and therefore it would be difficult to detect undescribed or novel mechanisms.

Reviewer #2 (Recommendations for the authors):

In this paper, de Nies and colleagues reanalyse their previously released, ground-breaking multi-omics time series dataset collected from a wastewater treatment plant microbial community over a 14 month period (with 51 sampling days), with the aim of exploring the type and diversity of AMR-related genes present in this microbial community, the type of mobile genetic element they are associated with (plasmid or phage), their association with microbial groups and the nature of their temporal dynamics.

This work holds strong interest and utility for workers in the areas of antimicrobial resistance and One Health, and, more broadly, for our understanding the nature of the interconnections between human and natural ecosystems.

The starting point is the metagenome assemblies constructed from the collected nucleic acid sequence datasets, and from which contigs form the units of analysis. The analysis is organised around classifying AMR related genes within contig sequences, with particular reference to distinguishing between the type of mobile genetic element (MGE) that harbour them, namely either plasmids or bacteriophage/prophages (assumed to be incorporated into a host genome via transduction). While there are definitive strengths of using the contig-level analysis, I feel far more detail could be provided in the methods description to convince the reader of the validity of the analysis, particularly in relation to drawing associations between MGE sequences and host microbial groups. Exactly how this is done, and the relative strengths and weaknesses, are not fully described and this limits the extent to which the reader can evaluate the overall approach and findings.

At certain places in the manuscript, I feel there is an over-reliance of visual data analysis, and more rigorous statistical analysis should be used to test the hypotheses in question; this includes (1) the analysis summarised on line 222 concerning the proportion of ARG/s attributable to plasmids and phages, and (2) the analyses presented in Figure 2C. I note there is also an analysis examining ESKAPEE pathogens, however, there is no direct evidence that the genome "backbone" of these organisms are actually present.

Materials and methods

For an analysis of considerable complexity, the methods section needs more detail and clarity, as it is sometimes hard to understand how a specific result as been obtained.

Firstly, the metagenome and metatranscriptome reads from all n=51 samples are co-assembled (line 462-463 ), but it is not clearly stated whether the both sets of data are assembled together?, This appears to be the case, as the reader is given the impression there is a single set of contigs in play? (line 466).

Secondly, a critical part of the methodology relates to the annotation of MGE sequences, and the inference of association between MGE sequences and bacterial groups; however. My understanding is that PathoFact will predict whether a contig arises, if classifiable, from either a plasmid or a phage, or whether it derives from a chromosomal genome. The taxonomic annotations for all contig are then made with Kraken. One possibility is that a contig arises from a chromosomal genome and contains an AMR sequence, which seems straightforward, as cognate AMR sequence will be unambiguously associated with the Kraken taxonomic annotation. Whereas another possibility is that the contig arises from a non-chromosomal MGE context; but in that case, it is a little unclear how the taxonomic assignments are actually made (i.e. over lines 477-489); particularly the logic that supports the statement "…to track transmission of AMR between taxa". It is very unclear how these decisions are actually made, but it comes across as rather subjective e.g. "…considering all different predictions of MGEs…". Needless to say, these inferences are critical for the rest of the paper, so can the authors please describe the approach in more detail, including backfilling the underlying scientific logic? A schematic diagram would be useful showing how different sources of data are linked.

Results

Figure 1. What is the level of read coverage that underpins the relative abundance measures used? Can some short statistical summary to cover that be included so the reader can calibrate?

Line 123. Changes evident at dates 13-05-2011 and 08-02-2012: perhaps these are better described as "transient" [changes]?

Line 127 (paragraph). To what extent are these changes simply secondary to abundance?

Line 158 "assemblies". This is the first place in the manuscript that reference is made to the fact that sequence data are assembled: please provide a sentence to cover this?

Line 160. Supp. Figure 2: what do the colours of the circles atop the bars denote? In general figure legends for Supp Figuresneed more detail.

Line 167+. AMR dynamics within Ca. M. parvicella. Line 175 "..confirmed

175 through metatranscriptomic analysis (Figure 2c): I cannot see how Figure 2C supports this statement. Line 179 "…29-22-2011…": I presume you mean 29-11-2011 [the second of the two adjacent points in Figure 2b that show near-zero signal levels?].

Line 185+ Figure 2C. In places, this figure is incomprehensible: can you offset the labels for improved clarity? "…clear and distinct co-occurrence patterns": there is no statistical support for this conclusion provided.

Line 194+. Where are the data that support the genomic "backbone" of these pathogens are detectable in this dataset? Also relevant for Figure 5C. The lack of protein level expression signals (line 316) from these organisms may also be a relevant observation in relation to this point.

Line 222. "…we found that plasmids contributed to an average of 10.8% of all ARGs, while phage contributed to an average of 6.8% of all resistance genes, confirming the general hypothesis that conjugation has the greatest influence on the dissemination of ARGs". Where is the statistical support for this conclusion? Do you think the effect size is scientifically significant?

Line 235. Figure 4B. What multiple correction procedure used and what were the total number of tests corrected for? (see also Figure 5A).

[Editors’ note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Mobilome-driven segregation of the resistome in biological wastewater treatment" for further consideration by eLife. Your revised article has been evaluated by Bavesh Kana (Senior Editor) and a Reviewing Editor.

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

Please revise the text to address remaining issues regarding the methodology and statistics, figures, legends and biological aspects to make the document more precise and clear. You will find further details in the reviewers comments below.

Reviewer #1 (Recommendations for the authors):

Thank your revision of the paper and responses to my comments, with have addressed many of my concerns. Many parts of the paper are now much clearer. For comments for which I still have concerns, I have listed below, as numbered:

Comment 2.7

In your response, you state "Numbers of taxa are highlighted at the top of each bar, with colors corresponding to the AMR categories on the x-axis.". This will only make sense if these colours are used consistently across figures (which they appear to be, but can you confirm or otherwise?). If so, some reference to that fact should be made.

Comment 2.9

Figure 2C remains challenging. The resolution appears to have increased but the nodes with long annotations are extremely difficult to read, as are the overlapping labels. I don't think it is unreasonable to say that this figure in its present form is not suitable for publication. Have you considered a heatmap (of the adjacency matrix) as an alternative?

Comments 2.12

Firstly I do not see how these tests can be described a two-way ANOVA, given there appears to be only one factor that is tested e.g. phage vs plasmid. Can you confirm that multiple testing has been performed using the total number of tests: that is, the number of ARG categories, encoded (Figure 4) or expressed (Figure 5)? In terms of visualisation, some of the horizontal "significance bars" are solid at both ends, while some show some kind of directionality. Can you clarify what these each mean?

Reviewer #2 (Recommendations for the authors):

In this revised version the authors have done a good job of addressing the reviewers' various concerns, such as rephrasing some claims and including additional needed methodological information and statistical support. As consequence, the manuscript and their interpretation of the results are clearer. There are still some small points to be addressed for further clarification of the text.

1. Despite making it clear in the rebuttal that sub-inhibitory antibiotic concentrations can in fact induce HGT, the same does not apply (as far as I know and despite the references shown) for ARGs. Resistant genes can be selected as a consequence of a particular pressure, but I see no indication that these genes actually "facilitate HGT of ARGs into new hosts" as is suggested in the text (L79). Please revise.

2. Is there a possible explanation for the fact that, for example, there were no markers associated with plasmids (or very few) in Salmonella, even though plasmids are in fact associated with ARGs? It might also be relevant to mention possible limitations of this approach (as in any approach) since this might help researchers in the future.

3. Lines 496-8 mention the lack of tools to link MGEs and hosts, but the authors should look at the following reference: https://doi.org/10.1101/198713.

4. L422 – In the response to a previous comment, the authors state that the approach can detect "mechanisms novel to a specific microbial reservoir". Since this phrase was not included in the text itself, the sentence is still misleading. They should make clear that when they refer to novel mechanisms, they refer to a specific environment being studied (particularly in L422; perhaps also on L 445).

https://doi.org/10.7554/eLife.81196.sa1

Author response

[Editors’ note: the authors resubmitted a revised version of the paper for consideration. What follows is the authors’ response to the first round of review.]

Reviewer #1 (Recommendations for the authors):

This work aims to provide a broad view of antimicrobial resistance determinants in a wastewater treatment plant by conducting a multi-omics approach on samples taken over 1.5 years. The use of several omics approaches and the distinction of markers associated with either phages or plasmids in the system represent strengths of the study. However, the lack of clarity and methodological details weaken the interpretation and impact of the results.

We recognize the reviewer’s appreciation of the strength of the study as well as their concerns regarding the lack of clarity and methodological details. To address these concerns, we have included additional information on the methodology. We have also thoroughly revised the manuscript to improve clarity as described in detail below.

1. While the multiple omics approach provides valuable insight, it is not entirely clear how expression data adds to what is already known regarding ARGs in the environment and how this expression information sheds new light on the resistome or can be "crucial for monitoring human health conditions" (L412). In this respect it would have been good to see if their findings on ARGs segregating between phages and plasmids correlate with what is known about resistances, at least in the relevant taxa (eg, L 288-9).

We thank the reviewer for this comment. Although metagenomics is useful to characterize the resistome, it does not provide information regarding the actual expression of antibiotic resistance genes (ARGs). In contrast, metatranscriptomics and metaproteomics identifies the ARGs encoded in the metagenome that are actively expressed within the community. For instance, expression indicates that the relevant ARGs are actively expressed as the organisms are exposed to the corresponding antibiotics. As such, linked with metagenomic data, metatranscriptomic data provides a taxonomic, genomic and functional overview of the BWWTP resistome. Thereby potentially reflecting the presence of ARGs in the general human population, as well as the further dissemination in the receiving environment. Moreover, metatranscriptomic data, in addition to metagenomic data, provides essential additional perspectives not only of the short-term adaptation to existing environmental conditions but also regarding selective pressures shaping the resistome.

Regarding the segregation of ARGs according to MGE, while to our knowledge no other study has assessed the differential contribution of MGEs to the dissemination of AMR, various studies have investigated the dissemination of AMR by MGEs in clinically relevant taxa. For instance, correlating with our findings a study by Gupta et al. (PMID: 23734150) found that phages are associated with aminoglycoside resistance, among others, in Salmonella enterica. Additionally, analogous to our findings, they found that plasmid encoded resistance is heavily associated with pathogenic bacteria such as S. pneumoniae, S. aureus, K. pneumoniae, E. kobie and E. hormaeche. This has been further highlighted in lines 381-384 in the revised manuscript.

2. Since the results are based on assembled data (L 158), the authors should provide additional and sufficient information regarding these assemblies and how they can be confident of taxa identification to the species level (such as the ESKAPEE group) and the linkage of ARGs to particular microorganisms. This is important, given that much of the analysis and ensuing conclusions are based on linkage of ARGs and MGEs to particular taxa. It would also be important to know how they managed to precisely map expression data to these MGEs and not to chromosomal loci and if perhaps this could also explain some of the discrepancies observed (eg, L397-8).

We thank the reviewer for highlighting this point. Based on this comment, we have restructured and included additional information regarding the methodology in the revised manuscript (lines 457-514). Additionally, we have included an additional figure representing a schematic overview of the methods used (Figure 8).

With respect to the reviewer’s first point relating to this comment, both metagenomic and metatranscriptomic reads were preprocessed and assembled using the Integrated Meta-omic Pipeline (IMP v1.3; available at https://r3lab.uni.lu/web/imp/, PMID: 27986083). Specifically, during processing Illumina Truseq2 adapters were trimmed, and reads of human origin were filtered out, which was followed by a de novo assembly using MEGAHIT (v1.2). For the assembly, both metagenomic and metatranscriptomic reads were co-assembled to increase the contiguity of the assemblies (lines 465-472, in the revised manuscript), in an iterative process to optimize read usage (PMID: 27986083).

The contig-based assembly file was then used as input for the taxonomic classification system Kraken2 which subsequently assigned the bacterial taxonomy, including where relevant the ESKAPEE group, to each contig. To ensure high levels of confidence in taxonomic classification, stringent cutoffs such as a minimum 70% identity of the contig with the database, across 90% of the contig length was used. Since the same contig-based assembly file was used as input for PathoFact, and as such for the prediction of antimicrobial resistance genes and mobile genetic elements, functional- and taxonomic- annotations were linked based on the contig they originated from (lines 487-496, in the revised manuscript).

Finally, to calculate the gene copy number and transcriptome expression, both the metagenomic and metatranscriptomic reads were independently mapped to the assembly files. The tool FeatureCount was then used to extract the raw read counts per contig as well as per individual ORF, as provided by PathoFact (lines 511-514). By using the above-described stringent methods and based on an identical approach in our previous (peer-reviewed) publications (PMIDs: 33597026, 35484157, DOI: 10.1038), we were able to precisely map the expression data to MGEs, including those distinct from the chromosomal loci, with high levels of confidence.

3. It seems that about 18% of all ARGs are associated with MGEs, yet no additional mention is made of the remaining genes, which presumably are chromosomal. Please clarify this in the document and also indicate how this could affect the transcriptomic and proteomic analysis.

We thank the reviewer for highlighting that the remaining resistance genes are found encoded on the bacterial chromosome and have expanded on this in the revised manuscript (lines 226-228).

Regarding the transcriptomic and proteomic analysis, these are not affected at all given that all expressed genes and proteins are taken into account (i.e. both chromosomal and MGE-derived genes).

4. It is intriguing that they observe high levels of peptide resistance in their metatranscriptome (L146) associated with resistance in E. coli, which is not a major contributor in the community and the gene was apparently not abundant in their metagenomic analysis. This discrepancy is not fully addressed and is relevant because it can also reflect biases in the approach or conditions of the sites being analyzed.

We thank the reviewer for picking up on this point. While resistance against the peptide antibiotic microcin J25 mediated by YojI has first been described in E. coli, we found this resistance gene to be widely distributed among Comamonadaceae members, among others, which are major members in the community. We have now adjusted the manuscript to reflect this fact in lines 149153.

5. It seems that the majority of markers are associated with abundant taxa, which is to be expected, but it is not clear if this is the reason why the authors hypothesize that the majority of ARGs would be in the major taxon M. parvicella (L169 and L175) or why these ARGs would confer fitness advantage (L184), unless of course this correlates with antibiotics present, which was not determined.

We thank the reviewer for highlighting this point. M. parvicella generally, as confirmed in previous studies, is found to be the major taxon of the studied wastewater microbial community (PMID: 25424998, 33077707, 28721231). While the mechanisms underlying the abundance of this taxon are numerous, we hypothesized that their success within the WWTP microbial community may also be driven by an increased abundance of antimicrobial resistance genes. Owing to the nature of the sample, i.e. wastewater derived from commercial and anthropogenic influences, we posited that the relatively higher abundance of M. parvicella may be driven by the presence of several ARGs, which are likely to be introduced into the WWTP via the intake system. This is further underlined by the expression levels of ARGs encoded by M. parvicella indicating a fitness advantage and, thus, allowing it to thrive in the general antibiotic rich environment of a wastewater treatment plant.

6. The document can also benefit from revision of the text, particularly the discussion, to clarify some claims and make the document easier to follow. There are some statements that do not accurately reflect either the scope of the work or the results obtained. Some of these include:

L79 – why would ARGs and sub-inhibitory antibiotic concentrations (which should not be lethal) promote HGT?

Various studies have reported both the development of antimicrobial resistance mechanisms as well as the selection of resistant strains (as cited in the manuscript: PMID: 27029295, 20159551, 14506034 and 7726525). Furthermore, there has been further evidence suggesting that subinhibitory concentrations may significantly increase the frequency of dissemination of many types of MGEs (PMID: 21845185, 3710955, 3343220 and 18359835).

L 96 – They do not really shed light on "the evolution" and at best suggest dissemination differences.

The passage has been rephrased as suggested by the reviewer (line 96).

L102 – The analysis suggests but does not necessarily "demonstrate" that MGEs are important drivers of AMR dissemination. In addition, it does not necessarily follow that "assessing the activity of the ARGs" can lead to understanding the underlying mechanisms.

The passage has been corrected as suggested by the reviewer (lines 101-103).

L186 – why population sizes? Did you actually look at association with the size?

We thank the reviewer for highlighting this point and have corrected the passage in the revised manuscript (lines 189-191).

L192 – Is this figure 3a?

This is indeed figure 3a. We thank the reviewer for pointing this out and have adjusted this in the revised manuscript.

On L202 it says that they "assessed the acquisition and dissemination of AMR". The data shows only associations, which is not necessarily acquisition and at best hints at modes of dissemination.

As suggested by the reviewer we have revised the paragraph to highlight the modes of dissemination instead of acquisition.

L292-94- the findings indicate but are not themselves a threat or have potential for dissemination. Please correct.

The paragraph has been corrected as suggested by the reviewer.

L296 – it does not seem that the proteomics validated transcriptomics data, but rather are more of a complementary analysis (L298). Likewise, on L416- It does not seem that proteomics data validated the abundance and expression findings.

The metaproteomic data confirmed the identified resistance mechanisms in the metagenomic data as functional active. However, we agree with the reviewer that the metaproteomic analysis should be seen more as a complementary analysis than validation and have revised the section as suggested.

L413 and 438 – Monitoring based on sequence analysis is based on comparisons using known data and therefore it would be difficult to detect undescribed or novel mechanisms.

As correctly noted by the reviewer, sequence analysis is based on known data and the corresponding databases. Although using machine learning based searches allows us to identify homologs of known antimicrobial resistance genes and mechanisms with current methods, it is still difficult to detect completely novel mechanisms. However, these methods allow for the early detection of mechanisms novel to a specific microbial reservoir. For example, the identified resistance gene Yoji, was already reported as a resistance mechanism against microcin J25 in 2005. Nevertheless, resistance was expected to be rare and microcin J25 has even recently been investigated as a potential new commercial antibiotic. Yet, within our BWWTP we not only identified the presence of Yoji, heretofore unknown in BWTTPs, we also identified high expression levels of this resistance mechanism.

Reviewer #2 (Recommendations for the authors):

In this paper, de Nies and colleagues reanalyse their previously released, ground-breaking multi-omics time series dataset collected from a wastewater treatment plant microbial community over a 14 month period (with 51 sampling days), with the aim of exploring the type and diversity of AMR-related genes present in this microbial community, the type of mobile genetic element they are associated with (plasmid or phage), their association with microbial groups and the nature of their temporal dynamics.

This work holds strong interest and utility for workers in the areas of antimicrobial resistance and One Health, and, more broadly, for our understanding the nature of the interconnections between human and natural ecosystems.

The starting point is the metagenome assemblies constructed from the collected nucleic acid sequence datasets, and from which contigs form the units of analysis. The analysis is organised around classifying AMR related genes within contig sequences, with particular reference to distinguishing between the type of mobile genetic element (MGE) that harbour them, namely either plasmids or bacteriophage/prophages (assumed to be incorporated into a host genome via transduction). While there are definitive strengths of using the contig-level analysis, I feel far more detail could be provided in the methods description to convince the reader of the validity of the analysis, particularly in relation to drawing associations between MGE sequences and host microbial groups. Exactly how this is done, and the relative strengths and weaknesses, are not fully described and this limits the extent to which the reader can evaluate the overall approach and findings.

At certain places in the manuscript, I feel there is an over-reliance of visual data analysis, and more rigorous statistical analysis should be used to test the hypotheses in question; this includes (1) the analysis summarised on line 222 concerning the proportion of ARG/s attributable to plasmids and phages, and (2) the analyses presented in Figure 2C. I note there is also an analysis examining ESKAPEE pathogens, however, there is no direct evidence that the genome "backbone" of these organisms are actually present.

Materials and methods

For an analysis of considerable complexity, the methods section needs more detail and clarity, as it is sometimes hard to understand how a specific result as been obtained.

Firstly, the metagenome and metatranscriptome reads from all n=51 samples are co-assembled (line 462-463 ), but it is not clearly stated whether the both sets of data are assembled together?, This appears to be the case, as the reader is given the impression there is a single set of contigs in play? (line 466).

We thank the reviewer for this comment. For the assembly, both metagenomic and metatranscriptomic reads were indeed co-assembled, i.e. assembled together, to increase the contiguity of the assemblies and enhance read usage (PMID: 27986083). Specifically, both metagenomic and metatranscriptomic reads were processed and assembled together in a single set of contigs using the Integrated Meta-omic Pipeline (IMP v1.3; available at https://r3lab.uni.lu/web/imp/). This has been extensively described in the article describing the IMP workflow (PMID: 27986083). However, to clarify this further as pointed out by the reviewer, we have expanded lines 468-473, in the revised manuscript.

Secondly, a critical part of the methodology relates to the annotation of MGE sequences, and the inference of association between MGE sequences and bacterial groups; however. My understanding is that PathoFact will predict whether a contig arises, if classifiable, from either a plasmid or a phage, or whether it derives from a chromosomal genome. The taxonomic annotations for all contig are then made with Kraken. One possibility is that a contig arises from a chromosomal genome and contains an AMR sequence, which seems straightforward, as cognate AMR sequence will be unambiguously associated with the Kraken taxonomic annotation. Whereas another possibility is that the contig arises from a non-chromosomal MGE context; but in that case, it is a little unclear how the taxonomic assignments are actually made (i.e. over lines 477-489); particularly the logic that supports the statement "…to track transmission of AMR between taxa". It is very unclear how these decisions are actually made, but it comes across as rather subjective e.g. "…considering all different predictions of MGEs…". Needless to say, these inferences are critical for the rest of the paper, so can the authors please describe the approach in more detail, including backfilling the underlying scientific logic? A schematic diagram would be useful showing how different sources of data are linked.

We thank the reviewer for this comment. As detailed in our response to comment 1.X. by reviewer 1, we have now provided much more methodological details in the revised manuscript. As indicated by the reviewer, PathoFact does indeed predict whether a contig arises, if classifiable, from either a plasmid or a phage through application of three different prediction methods (PlasFlow, DeepVirFinder and VirSorter). Subsequently, the bacterial taxonomy, such as the ESKAPEE pathogens, are assigned, as highlighted by the reviewer using the taxonomic classification system Kraken2. To ensure confidence in taxonomic classification, stringent cutoffs such as a minimum 70% identity of the contig with the database, across 90% of the contig length was used (lines 487-496). As per the reviewer’s suggestion, and as described in Comment 1.3, to clarify the linkage of different sources of data, a schematic diagram (Supp. Figure 9) has been added to the revised manuscript.

We would further like to clarify that though several tools exist for the identification of plasmids and/or phages, for example: metaplasmidspades, PlasFlow, DeepVirFinder and VIBRANT; tools for comprehensively identifying the respective bacterial hosts for MGEs do not exist. Simultaneously, the hosts for phages can be determined based on classification and linking of CRISPR and Spacer sequences between bacterial and phage sequences, respectively. Since most tools, i.e. SpacePHARER and PHIST, limit this analysis to high-quality viral sequences thus precluding a large fraction of putative identified phages; we did not pursue this approach. The hosts of plasmids, on the other hand, due to the competence of plasmids (e.g. self-replication, broad host range) the precise nature of plasmid ranges are difficult to identify. Therefore, to streamline this process, we used taxonomic associations as a means of linking bacteria with MGEs, given the stringent filtering/identification criteria used. Since we already used stringent cutoffs for plasmid/phage predictions, the original MGE predictions from PathoFact were used asis. This was followed by the taxonomic classification using Kraken2. We have now further clarified these points in the revised manuscripts in lines 493-514.

Results

Figure 1. What is the level of read coverage that underpins the relative abundance measures used? Can some short statistical summary to cover that be included so the reader can calibrate?

As suggested by the reviewer we have included additional information on the read coverage and depth in the revised manuscript (lines 472-473). We have further included more detailed information including tables or read counts and coverage in our git repository: https://git-r3lab.uni.lu/laura.denies/lao_scripts.

Line 123. Changes evident at dates 13-05-2011 and 08-02-2012: perhaps these are better described as "transient" [changes]?

We thank the reviewer for highlighting this and have revised the description as suggested.

Line 127 (paragraph). To what extent are these changes simply secondary to abundance?

As suggested by the reviewer, the changes in relation to the presence and absence of specific AMR categories may indeed to an extent be secondary to abundance. For instance, as would be expected, highly abundant AMR categories such as genes conferring resistance to aminoglycoside, β-lactam as well as multidrug resistance were found across all timepoints. Meanwhile, lowly abundant AMR categories such as bicyclomycin resistance were found only sporadically within the BWWTP, indicating an association between AMR abundance and prevalence throughout time. However, one should also note the identified exceptions such as fosfomycin resistance, which, although less abundant, was found to be present across all timepoints, demonstrating that even lowly abundant AMR categories may still be prevalent over time.

Line 158 "assemblies". This is the first place in the manuscript that reference is made to the fact that sequence data are assembled: please provide a sentence to cover this?

The paragraph has been revised as suggested (lines 117-121).

Line 160. Supp. Figure 2: what do the colours of the circles atop the bars denote? In general figure legends for Supp Figuresneed more detail.

The different color circles on top of the bars in Supp. Figure 2. denote the different categories of AMR. As suggested by the reviewer, we have adjusted the figure legends in the revised manuscript accordingly.

Line 167+. AMR dynamics within Ca. M. parvicella. Line 175 "...confirmed 175 through metatranscriptomic analysis (Figure 2c): I cannot see how Figure 2C supports this statement. Line 179 "…29-22-2011…": I presume you mean 29-11-2011 [the second of the two adjacent points in Figure 2b that show near-zero signal levels?].

We thank the reviewer for highlighting these points. In agreement with the metagenomic data, ARGs (e.g. aminoglycoside, β-lactam and multi-drug resistance) were found to be highly expressed by M. parvicella, likely driven by selective pressure within the antibiotic-rich environment of the BWWTP. We have included an additional supplementary figure to highlight expressions levels of ARGs encoded by M. parvicella (Supp. Figure 3 in the revised manuscript). Additionally, we have corrected line 179 (in the revised manuscript) as suggested by the reviewer.

Line 185+ Figure 2C. In places, this figure is incomprehensible: can you offset the labels for improved clarity? "…clear and distinct co-occurrence patterns": there is no statistical support for this conclusion provided.

We thank the reviewer for highlighting this and have adjusted the figure labels. Additionally, we used edges with correlation values of greater than 0.6 correcting for multiple testing where the adjusted p-value <0.05. We have now revised the methods to include this information (lines 565576).

Line 194+. Where are the data that support the genomic "backbone" of these pathogens are detectable in this dataset? Also relevant for Figure 5C. The lack of protein level expression signals (line 316) from these organisms may also be a relevant observation in relation to this point.

We thank the reviewer for the comment. As highlighted in our response to comments 1.2 and 2.2 above, the genomic backbone of the pathogens was based on taxonomic classifications of the contigs using stringent cutoffs. We have now included said data, i.e. taxonomic classification of pathogens and the resistance genes they encode, within our publicly available git repository: https://git-r3lab.uni.lu/laura.denies/lao_scripts.

With respect to the protein level expression, we agree that an overall lack of protein expression is relevant, and therefore now discuss this in detail in the Discussion section in lines 430-433. In this context, there are many factors which affect the detection of proteins, not least the presence of proteases in the WWTP environment. In some cases, the proteins secreted may act as lytic molecules or potentially even degraded. Dreyer et al. (https://doi.org/10.1038/s41598-019-478439) reported the proteolytic degradation of bacteriocins, positing the possibility that some of the proteins may be undetectable due to technical reasons. The revised discussion now includes details reflecting these circumstances.

Line 222. "…we found that plasmids contributed to an average of 10.8% of all ARGs, while phage contributed to an average of 6.8% of all resistance genes, confirming the general hypothesis that conjugation has the greatest influence on the dissemination of ARGs". Where is the statistical support for this conclusion? Do you think the effect size is scientifically significant?

We thank the reviewer for highlighting this point. We found the contribution of plasmids contributing to AMR significantly increased compared to phages (p-value <0.05, Two-way ANOVA) and have updated the revised manuscript with the required statistical support (lines 228233).

Line 235. Figure 4B. What multiple correction procedure used and what were the total number of tests corrected for? (see also Figure 5A).

To determine the statistical difference of the results visualized in Figure 4b and Figure 5a, a two-way ANOVA has now been performed, after which the multiple comparison procedure Tukey HSD was performed using the fitted ANOVA as an argument.

[Editors’ note: what follows is the authors’ response to the second round of review.]

The manuscript has been improved but there are some remaining issues that need to be addressed, as outlined below:

Please revise the text to address remaining issues regarding the methodology and statistics, figures, legends and biological aspects to make the document more precise and clear. You will find further details in the reviewers comments below.

Reviewer #1 (Recommendations for the authors):

Thank your revision of the paper and responses to my comments, with have addressed many of my concerns. Many parts of the paper are now much clearer. For comments for which I still have concerns, I have listed below, as numbered:

We thank the reviewer for the recognition of our work and have addressed the remaining concerns in the comments below.

Comment 2.7

In your response, you state "Numbers of taxa are highlighted at the top of each bar, with colors corresponding to the AMR categories on the x-axis.". This will only make sense if these colours are used consistently across figures (which they appear to be, but can you confirm or otherwise?). If so, some reference to that fact should be made.

As the reviewer noted the colors corresponding to the AMR categories are indeed used consistently across all figures. As suggested, we have included a reference to this in the figure legends.

Comment 2.9

Figure 2C remains challenging. The resolution appears to have increased but the nodes with long annotations are extremely difficult to read, as are the overlapping labels. I don't think it is unreasonable to say that this figure in its present form is not suitable for publication. Have you considered a heatmap (of the adjacency matrix) as an alternative?

We appreciate the reviewer’s concerns regarding Figure 2c. We have considered a heatmap of the adjacency matrix as an alternative, but this resulted in the loss of information provided by the clusters as observed in the network. However, we have revised the original network and adjusted the long annotations as well as the overlapping labels to improve readability. ARGs have been color coded by the corresponding AMR category, with further information regarding specific ARGs provided in a supplementary adjacency matrix (https://git-r3lab.uni.lu/laura.denies/lao_scripts//tree/master/Data). Additionally, the fonts and size of the taxa annotations have been adjusted to improve readability. These adjustments should address the reviewer’s concerns regarding this figure.

Comments 2.12

Firstly I do not see how these tests can be described a two-way ANOVA, given there appears to be only one factor that is tested e.g. phage vs plasmid. Can you confirm that multiple testing has been performed using the total number of tests: that is, the number of ARG categories, encoded (Figure 4) or expressed (Figure 5)? In terms of visualisation, some of the horizontal "significance bars" are solid at both ends, while some show some kind of directionality. Can you clarify what these each mean?

We thank the reviewer for highlighting these points. Regarding the two-way ANOVA, the different ARG categories were taken into account as an additional factor. More specifically, in the metagenomic data, a total of 28 ARG categories were tested, of which 6 (aminoglycoside, bacitracin, fosfomycin, MLS, peptide and sulfonamide resistance) were significantly differentially encoded on phages in comparison to plasmids. In the metatranscriptomic data, the same 28 ARG categories were tested of which 6 (aminoglycoside, bacitracin, glycopeptide, mupirocin, peptide and sulfonamide resistance) were significantly differentially expressed on phages versus plasmids. We have clarified this further in the text by further elaborating on the description of the statistical analyses (lines 242-244 and lines 584-586).

Regarding the visualizations, the bars indicating significance should be straight and solid at both ends. We thank the reviewer for pointing out this artifact. We have fixed it in the revised manuscript.

Reviewer #2 (Recommendations for the authors):

In this revised version the authors have done a good job of addressing the reviewers' various concerns, such as rephrasing some claims and including additional needed methodological information and statistical support. As consequence, the manuscript and their interpretation of the results are clearer. There are still some small points to be addressed for further clarification of the text.

We thank the reviewer for their appreciation of the revised manuscript. Please see below our responses to the remaining comments.

1. Despite making it clear in the rebuttal that sub-inhibitory antibiotic concentrations can in fact induce HGT, the same does not apply (as far as I know and despite the references shown) for ARGs. Resistant genes can be selected as a consequence of a particular pressure, but I see no indication that these genes actually "facilitate HGT of ARGs into new hosts" as is suggested in the text (L79). Please revise.

We thank the reviewer for highlighting this point and appreciate their concern. Sub-inhibitory concentrations of antibiotics can induce HGT as well as select for genes conferring resistance and/or induce resistant mutations (PMID: 26925045). Additionally, as shown by Datta and Hughes (PMID: 6316165), selection pressures may alter HGT processes, increasing the number of resistance elements which reside on mobile DNA, thereby indirectly facilitating the transmission of said resistance elements to new hosts. However, we appreciate the reviewer’s concerns and have revised the manuscript to clarify and elaborate on these aspects (lines 78-83).

2. Is there a possible explanation for the fact that, for example, there were no markers associated with plasmids (or very few) in Salmonella, even though plasmids are in fact associated with ARGs? It might also be relevant to mention possible limitations of this approach (as in any approach) since this might help researchers in the future.

We thank the reviewer for bringing up this point. While Salmonella was identified throughout the time series in the BWWTP, the overall levels of detection for this genus were quite low (relative abundance of 0.01-0.02%). In addition, resistance genes encoded on the Salmonella genome were only observed at 5 timepoints, which may be related to the low abundances within this BWWTP. This may also be the reason that, any resistance-encoding plasmids associated with Salmonella were undetected.

As highlighted by the reviewer in comment #3 (below), there are Hi-C based methods which provide the capacity to link MGEs and hosts during sample preparations and subsequent sequencing (PMID: 19815776 and https://doi.org/10.1101/198713). We have now included this aspect in the Discussion in lines 410-418. Other methods include targeted sequencing of plasmids of interest using the methodology established by Li et al. (PMID: 29325009), which has been highlighted in the revised Discussion section in lines 402-409.

3. Lines 496-8 mention the lack of tools to link MGEs and hosts, but the authors should look at the following reference: https://doi.org/10.1101/198713.

We thank the reviewer for highlighting this point. While Hi-C-based methods allow capturing of inter-chromosomal junctions, such as plasmid-genome interactions, it requires a different processing of the samples prior to sequencing to induce covalent linkage among DNA molecules. As such we did not include said methods in our manuscript. However, we agree with the reviewer that Hi-C-based methods provide an important method to link MGEs and hosts and have rephrased the corresponding section accordingly (lines 410-418).

4. L422 – In the response to a previous comment, the authors state that the approach can detect "mechanisms novel to a specific microbial reservoir". Since this phrase was not included in the text itself, the sentence is still misleading. They should make clear that when they refer to novel mechanisms, they refer to a specific environment being studied (particularly in L422; perhaps also on L 445).

We thank the reviewer for highlighting this point and have now revised the paragraph as suggested (lines 439-441 and lines 465-468).

https://doi.org/10.7554/eLife.81196.sa2

Article and author information

Author details

  1. Laura de Nies

    Systems Ecology Group, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6483-7489

  2. Susheel Bhanu Busi

    Systems Ecology Group, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg
    Contribution
    Resources, Data curation, Formal analysis, Investigation, Visualization, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7559-3400

  3. Benoit Josef Kunath

    Systems Ecology Group, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg
    Contribution
    Formal analysis
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-3356-8562

  4. Patrick May

    Bioinformatics Core, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg
    Contribution
    Data curation, Software, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8698-3770

  5. Paul Wilmes

    1. Systems Ecology Group, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg
    2. Department of Life Sciences and Medicine, Faculty of Science, Technology and Medicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg
    Contribution
    Conceptualization, Resources, Software, Supervision, Funding acquisition, Validation, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    paul.wilmes@uni.lu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6478-2924

Funding

Fonds National de la Recherche Luxembourg (CORE/13684739)

  • Paul Wilmes

European Research Council (ERC-CoG 863664)

  • Paul Wilmes

Fonds National de la Recherche Luxembourg (PRIDE/11823097)

  • Paul Wilmes
  • Laura de Nies

Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (CRSII5_180241)

  • Susheel Bhanu Busi

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We are thankful for the assistance of Audrey Frachet Bour, Lea Grandmougin, Janine Habier, Laura Lebrun (LCSB) for laboratory support. We acknowledge the valuable input from Rashi Halder at the LCSB Sequencing Platform with respect to library preparation. The computational analyses were performed at the HPC facilities at the University of Luxembourg (Varrette et al., 2014).

Senior Editor

  1. Bavesh D Kana, University of the Witwatersrand, South Africa

Reviewing Editor

  1. María Mercedes Zambrano, CorpoGen, Colombia

Reviewer

  1. María Mercedes Zambrano, CorpoGen, Colombia

Version history

  1. Preprint posted: November 15, 2021 (view preprint)
  2. Received: June 19, 2022
  3. Accepted: September 15, 2022
  4. Accepted Manuscript published: September 16, 2022 (version 1)
  5. Version of Record published: November 8, 2022 (version 2)

Copyright

© 2022, de Nies et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,301
    Page views
  • 299
    Downloads
  • 3
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Laura de Nies
  2. Susheel Bhanu Busi
  3. Benoit Josef Kunath
  4. Patrick May
  5. Paul Wilmes
(2022)
Mobilome-driven segregation of the resistome in biological wastewater treatment
eLife 11:e81196.
https://doi.org/10.7554/eLife.81196

Categories and tags

Further reading

    1. Microbiology and Infectious Disease

    Antimicrobial Resistance: Adapt or perish

    Rohan BH Williams
    Insight

    Microbial communities in wastewater treatment plants provide insights into the development and mechanisms of antimicrobial resistance.

    1. Genetics and Genomics
    2. Microbiology and Infectious Disease

    Enterobacterales plasmid sharing amongst human bloodstream infections, livestock, wastewater, and waterway niches in Oxfordshire, UK

    William Matlock, Samuel Lipworth ... REHAB Consortium
    Research Article Updated

    Plasmids enable the dissemination of antimicrobial resistance (AMR) in common Enterobacterales pathogens, representing a major public health challenge. However, the extent of plasmid sharing and evolution between Enterobacterales causing human infections and other niches remains unclear, including the emergence of resistance plasmids. Dense, unselected sampling is essential to developing our understanding of plasmid epidemiology and designing appropriate interventions to limit the emergence and dissemination of plasmid-associated AMR. We established a geographically and temporally restricted collection of human bloodstream infection (BSI)-associated, livestock-associated (cattle, pig, poultry, and sheep faeces, farm soils) and wastewater treatment work (WwTW)-associated (influent, effluent, waterways upstream/downstream of effluent outlets) Enterobacterales. Isolates were collected between 2008 and 2020 from sites <60 km apart in Oxfordshire, UK. Pangenome analysis of plasmid clusters revealed shared ‘backbones’, with phylogenies suggesting an intertwined ecology where well-conserved plasmid backbones carry diverse accessory functions, including AMR genes. Many plasmid ‘backbones’ were seen across species and niches, raising the possibility that plasmid movement between these followed by rapid accessory gene change could be relatively common. Overall, the signature of identical plasmid sharing is likely to be a highly transient one, implying that plasmid movement might be occurring at greater rates than previously estimated, raising a challenge for future genomic One Health studies.

    1. Computational and Systems Biology
    2. Microbiology and Infectious Disease

    Candida albicans exhibits heterogeneous and adaptive cytoprotective responses to antifungal compounds

    Vanessa Dumeaux, Samira Massahi ... Michael T Hallett
    Research Article Updated

    Candida albicans, an opportunistic human pathogen, poses a significant threat to human health and is associated with significant socio-economic burden. Current antifungal treatments fail, at least in part, because C. albicans can initiate a strong drug tolerance response that allows some cells to grow at drug concentrations above their minimal inhibitory concentration. To better characterize this cytoprotective tolerance program at the molecular single-cell level, we used a nanoliter droplet-based transcriptomics platform to profile thousands of individual fungal cells and establish their subpopulation characteristics in the absence and presence of antifungal drugs. Profiles of untreated cells exhibit heterogeneous expression that correlates with cell cycle stage with distinct metabolic and stress responses. At 2 days post-fluconazole exposure (a time when tolerance is measurable), surviving cells bifurcate into two major subpopulations: one characterized by the upregulation of genes encoding ribosomal proteins, rRNA processing machinery, and mitochondrial cellular respiration capacity, termed the Ribo-dominant (Rd) state; and the other enriched for genes encoding stress responses and related processes, termed the Stress-dominant (Sd) state. This bifurcation persists at 3 and 6 days post-treatment. We provide evidence that the ribosome assembly stress response (RASTR) is activated in these subpopulations and may facilitate cell survival.