Skip to main content

Disentangling abiotic and biotic effects of treated wastewater on stream biofilm resistomes enables the discovery of a new planctomycete beta-lactamase

Abstract

Background

Environmental reservoirs of antibiotic resistance pose a threat to human and animal health. Aquatic biofilms impacted by wastewater effluent (WW) are known environmental reservoirs for antibiotic resistance; however, the relative importance of biotic factors and abiotic factors from WW on the abundance of antibiotic resistance genes (ARGs) within aquatic biofilms remains unclear. Additionally, experimental evidence is limited within complex aquatic microbial communities as to whether genes bearing low sequence similarity to validated reference ARGs are functional as ARGs.

Results

To disentangle the effects of abiotic and biotic factors on ARG abundances, natural biofilms were previously grown in flume systems with different proportions of stream water and either ultrafiltered or non-ultrafiltered WW. In this study, we conducted deep shotgun metagenomic sequencing of 75 biofilm, stream, and WW samples from these flume systems and compared the taxonomic and functional microbiome and resistome composition. Statistical analysis revealed an alignment of the resistome and microbiome composition and a significant association with experimental treatment. Several ARG classes exhibited an increase in normalized metagenomic abundances in biofilms grown with increasing percentages of non-ultrafiltered WW. In contrast, sulfonamide and extended-spectrum beta-lactamase ARGs showed greater abundances in biofilms grown in ultrafiltered WW compared to non-ultrafiltered WW. Overall, our results pointed toward the dominance of biotic factors over abiotic factors in determining ARG abundances in WW-impacted stream biofilms and suggested gene family-specific mechanisms for ARGs that exhibited divergent abundance patterns. To investigate one of these specific ARG families experimentally, we biochemically characterized a new beta-lactamase from the Planctomycetota (Phycisphaeraceae). This beta-lactamase displayed activity in the cleavage of cephalosporin analog despite sharing a low sequence identity with known ARGs.

Conclusions

This discovery of a functional planctomycete beta-lactamase ARG is noteworthy, not only because it was the first beta-lactamase to be biochemically characterized from this phylum, but also because it was not detected by standard homology-based ARG tools. In summary, this study conducted a metagenomic analysis of the relative importance of biotic and abiotic factors in the context of WW discharge and their impact on both known and new ARGs in aquatic biofilms.

Video Abstract

Background

Antibiotic resistance is a major threat to public health. Among the numerous factors contributing to antibiotic resistance, anthropogenic pollution plays a role in altering the abundance and prevalence of antibiotic-resistance genes (ARGs) in environmental reservoirs [1,2,3,4]. The rapid increase in the consumption of antibiotics and the release of micropollutants (MPs) into the environment stimulates the spread of antibiotic-resistance genes (ARGs), posing a threat to modern medical procedures and treatments that rely on effective antibiotics [5].

Aquatic ecosystems that receive effluent from wastewater treatment plants (WWTPs) are influenced by biotic factors, such as antibiotic-resistant bacteria released from WWTPs [6], and abiotic factors [7, 8] including nutrients and MPs that are not completely eliminated by conventional WWTPs [9]. Several recent studies have demonstrated that exposure to WWTP effluent alters the taxonomic and ARG compositions of microbial communities in downstream rivers [10,11,12,13,14]. However, disentangling the effects of abiotic and biotic factors on ARG abundance is challenging in the environment. Here, to systematically distinguish abiotic and biotic factors, we built on a previous research study in which controlled flume systems were used to grow biofilms in ultrafiltered and non-ultrafiltered WWTP effluent diluted with different proportions of stream water [10, 11]. Ultrafiltration (UF) treatment (pore size 0.4 µm) removed most of the microorganisms and particles (biotic factors) while maintaining consistent nutrient and MP levels, including antibiotics. Therefore, separating biotic factors from abiotic factors in UF treatment enabled the investigation of the factors shaping the “resistome”, i.e., the total composition of ARGs in stream biofilms influenced by WWTP discharge.

Many studies on stream biofilms have relied on amplicon sequencing methods for taxonomic and/or ARG profiling [10, 15,16,17]. In contrast to amplicon sequencing, short-read and long-read shotgun metagenomic sequencing studies of streams and stream biofilms [13, 18, 19] do not introduce polymerase chain reaction (PCR) amplification bias. As we demonstrated in this study, PCR-free metagenomic sequencing also enables the identification of sequence-divergent ARGs not detected by conventional primer pairs. Amplicon sequencing is more sensitive because it can aid in the detection of low-abundance species, but it is also biased toward the amplification of specific taxa during PCR [20, 21]. Specifically, PCR bias limits the detection of some bacterial phyla such as Planctomycetota, for which the commonly used 16S rRNA PCR primers do not properly bind [22, 23]. The Planctomycetota are undersampled by classical amplicon-based methods [23] and are a focus of this study because of their relevance in WWTP [24], lacustrine [25], and riverine biofilm (periphyton) communities [26]. We conducted deep shotgun metagenomic sequencing of stream biofilms and water samples in this study, which enabled us to sample less well-studied microorganisms and their functional genes including ARGs.

Here, we demonstrated a combined computational and biochemical workflow to investigate the effects of wastewater effluent (WW) on the microbiome taxonomic composition, functional enzyme diversity, and ARG abundance in flume biofilms grown in natural stream water mixed with different proportions of ultrafiltered and non-ultrafiltered WW [10, 11]. Specifically, we leveraged the functional genes identified through shotgun metagenomics to heterologously express and purify an uncharacterized ARG from the phylum Planctomycetota. Despite sharing low sequence identity with computationally-identified ARGs, the planctomycete-derived enzyme displayed beta-lactamase activity. Overall, this study enabled the dissection of biotic and abiotic factors that affect stream biofilm resistors and led to the characterization of both known and new ARGs.

Methods

Biofilm cultivation and sampling

The flume systems and biofilm cultivation experiments which generated the samples sequenced in this study were conducted and described previously [10, 11]. Briefly, the system consisted of an indoor experimental flow-through channel system utilizing two pumps that controlled the flow rate to a maximum of 10 m3/h and fed two 1.2 m3 buffer tanks. These channels were fed with different mixtures of two different water sources: (1) a peri-urban stream (Chriesbach; Dübendorf, Switzerland) and (2) wastewater effluent from a pilot-size wastewater treatment plant (100-person equivalents). Stream water was mixed with wastewater effluent  at the following percentages: 0% (control), 30% and 80% WW non-ultrafiltered as well as 30% and 80% ultrafiltered WW (UF) using a filter with a pore size of 0.4 µm to remove particles, e.g., suspended solids and microorganisms. Twenty flow-through channels, each containing five glass slides were assigned using a randomized design to one of the five treatments (0%, 30% WW, 80% WW, 30% UF WW, 80% UF WW), resulting in four biological flume replicates per treatment. Biofilms that naturally formed on the slides under each of the different treatment conditions were allowed to develop for a period of 4 weeks. Automated online monitoring systems were used for parameters such as flow rate, water conductivity, and temperature. Stream and WW mixing ratios were adjusted daily based on conductivity measurements [10, 11]. An extended panel of 20 different physicochemical parameters (e.g., micropollutant concentrations, nutrients) was monitored weekly through grab samples as well as two passive samplers [10, 11].

DNA extraction and metagenomics assessment

Biofilms were scraped off of each of the five glass slides per channel and suspended in artificial river water, of which 2 mL were used for DNA extraction as described previously [11, 27]. Filters from 100 mL of stream water, WW, and UF WW samples were also included for sequencing. Total genomic DNA was extracted with a DNeasy PowerBiofilm Kit (QIAGEN) for a total of 75 samples (Table S1). The quality and quantity of the DNA extracts were measured using a NanoDrop One spectrophotometer. Extracted DNA samples were stored at − 20 °C and shotgun metagenomic sequencing was performed using the Illumina platform NovaSeq with a paired-end (2 × 150 bp; an average of 20 Gbp raw data per sample) strategy. Library construction and sequencing were performed by Novogene (Cambridge, UK). The reads obtained from shotgun metagenomics and those containing adapters and low-quality reads (N > 10% and quality score ≤ 5) were removed by Novogene, and the read qualities were double-checked by the authors using Fastp v.0.20.0 (https://github.com/OpenGene/fastp). De novo assembly was performed using MEGAHIT v1.2.9 [28]. Open reading frames (ORFs) were predicted from the assembled contigs using Prodigal v2.6.3 [29].

Taxonomic and functional profiling

The taxonomic composition was profiled using the degenerate consensus marker gene profiling tool, mTAGs v1.0.4 with 1000 maximum accepts and 1000 maximum rejects [30] using the cleaned, paired-end short reads. The mTAGs profiles for all samples were merged using the mTAGs merge function. Genus- and phylum-level per-sample count tables were used for subsequent analysis. LEfSe (linear discriminant analysis effect size) analysis [31] was used to identify taxa that significantly differed between treatment groups. For analysis of functional genes, counts corresponding to Enzyme Commission (EC) number classes were compiled as described previously [32]. Briefly, fourth-level EC number counts were calculated using the EC-annotated prokaryotic fraction of the UniProt database using DIAMOND blastx v2.0.6 [33] and a minimum bit score cutoff of 50 (all other parameters set to default).

Reads-based ARG detection and ordination analyses

The ARG resistome was profiled using the DeepARG v1.0.2 short reads pipeline with default settings [34]. ARG counts normalized by the 16S-rRNA counts within the output files for each sample were merged into an abundance table. All data were analyzed and visualized in R v4.3.0 with the tidyverse package [35]. Principal coordinate analysis (PCoA) of the microbiome composition was performed using the genus abundance table. Unaligned and unassigned reads were removed. Genus counts were normalized using the trimmed mean of M-values (TMM) method within the EdgeR package [36,37,38]. Normalized genus counts were square-root-transformed. PCoA of the resistome composition was calculated from 16S-normalized DeepARG abundances [34]. PCoA analysis of functional enzyme composition was calculated from normalized and Hellinger-transformed EC counts. For all the data sets, pairwise Bray–Curtis dissimilarity between samples were computed using the “vegdist” function from the vegan package [39]. PCoA was performed using the “pcoa” function from the ape package [40]. Procrustes analysis comparing the PCoA vectors was conducted using the ‘procrust’ and ‘protest’ functions from the vegan package [39]. The Kruskal–Wallis test was applied to test for significant correlations between the normalized average coverage of different ARGs and wastewater effluent proportions (WW%), and p values were corrected for multiple hypothesis testing using the Benjamini–Hochberg method [41].

ARG identification and coverage analysis from metagenome assemblies

For a complementary and separate approach to the metagenomic short reads pipeline for ARG analysis, predicted open reading frames (ORFs) from assembled metagenomes were identified by protein BLAST v2.9.0 [42] search against the Structured Antibiotic Resistance Genes (SARG) protein database v2.2, with stringent e-values and query coverage cutoffs of 0.1 and 20, respectively [43]. Using thresholds established previously in similar samples [44] for the detection of ARGs in riverine metagenomes, BLAST hits with more than 85% amino acid identity and an alignment length of 100 amino acids were retained. BLAST hits were then clustered using CD-HIT v4.6.8 [45, 46], with parameters –n 5 and –c 0.7. The 102 non-singleton cluster representatives of SARG hits were selected for further analysis. Metagenomic coverage of selected ARG cluster representatives was calculated across all metagenomes based on methods previously described [44]. Briefly, raw reads were mapped back to BLAST hits using Bowtie2 v2.4.4 [47] with the parameters shown in Table S2. Coverage depth was calculated using samtools v1.12 and average coverage was calculated as described previously [48]. The average coverage of the resulting 102 ORFs was normalized to reads per million as previously described [44].

Bioinformatic analyses of beta-lactamases

Metallo-beta lactamase fold proteins were identified by querying assembled proteins from the 75 metagenomes and additional publicly available metagenomes from MP enrichment studies (NCBI accession: PRJNA725625) using BLAST v.2.9.0 [49] with relaxed parameters (e value = 0.1, qcov_hsp_perc = 20, all other parameters set to default) to detect homologs with low sequence identity to known ARGs. Taxonomic diversity and trends in abundances of metallo-beta lactamase fold homologs across different WW and UF WW conditions were examined leading to the prioritization of two metallo-beta-lactamase homologs (NCBI accession numbers MBX3358097.1 and MBX3359331.1) for experimental characterization. Structural modeling of MBX3358097.1 was conducted using ColabFold AlphaFold 2.0 [50] with default parameters. AlphaFold models were assessed for quality based on predicted local distance difference test (pLDDT) scores and aligned pairwise with the reference structure in the Protein Data Bank (PDB ID: 3NKS) using the PyMOL “super” command with default parameters.

Cloning and protein expression

For functional characterization, two full-length Phycisphaeraceae (Planctomycetota) proteins (NCBI accession numbers MBX3358097.1 and MBX3359331.1) were codon-optimized for E. coli codon usage using the Integrated DNA Technologies Codon Optimization Tool (https://eu.idtdna.com/pages/tools/). Signal peptides were detected using SignalP v6.0 [51] and removed where present (Table S3). The genes were synthesized as gBlocks by Integrated DNA Technologies. The genes were cloned individually by Gibson Assembly into the multiple cloning site 1 of the pETDuet-1 vector (ampR resistance marker) with N-terminal hexahistidine tags (Table S3). The constructs were transformed into E. coli DH5α cells and verified by Sanger sequencing (Microsynth, Balgach, Switzerland). The sequence-verified constructs were subsequently transformed into E. coli BL21 (DE3) cells (New England Biolabs Frankfurt, Germany). Starter cultures (5 mL) were grown in LB with 100 μg/ml ampicillin overnight at 37 °C. From the starter cultures, 2.5 mL was added to 400 mL LB with 100 μg/mL ampicillin in a 1-L baffled Erlenmeyer flask and grown to an optical density of 0.5–0.8 at 37 °C. Protein expression was induced by adding 0.1 mM isopropyl β-D-1-thiogalactopyranoside (IPTG) and cultures were incubated for 18 h at 16 °C.

Protein purification

Expression cultures were poured into 50 mL Falcon tubes and centrifuged at 4000 RCF for 15 min. Supernatants were decanted, and cell pellets were resuspended in 4 mL buffer A (20 mM Tris, 500 mM NaCl, 20 mM imidazole, 10% glycerol, pH = 8) per gram of pellet weight. Cells were lysed by 2 cycles of sonication using a 6-mm needle, 20% amplitude and 10 × 10 s pulses with a rest time of 10 s in between (total sonication time per cycle = 3 min). Lysates were centrifuged at 20,000 RCF at 4 °C for 1 h. Cleared lysate was filtered through a 0.2-μm syringe filter, loaded into a ÄKTA Pure system (Cytiva Europe GmbH, Glattbrugg, Switzerland), and injected onto a HisTrap FF 5 mL column (Cytiva, Europe GmbH, Glattbrugg, Switzerland). The column was washed with 20 column volumes of buffer A. His-tagged proteins were eluted in 3 column volumes of buffer B (same as buffer A but with an imidazole concentration of 500 mM), and 8 × 2 mL fractions were collected. The fractions with the highest protein concentrations were determined qualitatively by Bradford assay [52]. Imidazole was removed by running the sample over a PD-10 desalting column (Cytiva Europe GmbH, Glattbrugg, Switzerland) and eluting in buffer (5 mM Tris, 30 mM NaCl, 10% glycerol, 0.1 mM ZnCl2, pH = 8). ZnCl2 was added because the enzyme of interest is suspected to be Zn2+ dependent. The purified protein was flash-frozen using liquid nitrogen in 200 μl aliquots and subsequently stored at − 80 °C.

β-Lactamase activity assay

The purified, soluble protein (MBX3358097.1) was tested with standard β-lactamase activity detection assays using the chromogenic cephalosporin analog, nitrocefin (Sigma Aldrich) [53]. Nitrocefin solutions were prepared in a 40-mM stock solution in DMSO and then diluted in 20 mM Tris buffer (pH = 7.0) immediately before use in a 1 mM working solution. Briefly, standard nitrocefin assays were performed in a 96-well plate by combining 180 μl 20 mM Tris buffer with 0.1 mM ZnCl2 (pH = 7.0) with 10 μL purified enzyme and 10 μL nitrocefin working solution. The plate was incubated in a plate reader (BioTek Synergy H1, Agilent) at 25 °C, and the absorbance at 482 nm was measured every minute for 2 h. As negative controls, boiled enzyme, and buffer-only conditions were included in every assay. To obtain the boiled enzyme control, aliquots of purified protein were incubated at 100 °C for 15 min. Reaction rates were quantified from continuous measurements as described previously by Ryu et al. using the molar extinction coefficient of nitrocefin (12,124 M−1 cm−1) in 20 mM Tris [54].

Results and discussion

Biotic factors are more dominant than abiotic factors for altering the taxonomic and functional composition of biofilms grown in WWTP effluent

We deeply sequenced metagenomes from a total of 75 biofilm and water samples to an average sequencing depth of 138 million reads per sample (max. 204 million reads). After metagenomic assembly, we analyzed the taxonomic composition of stream biofilm communities to assess whether community composition was correlated with different experimental treatments: 30% WW, 80% WW, 30% UF WW, 80% UF WW, or stream water (0% WW). The experimental treatment was correlated with taxonomic composition for both prokaryotes (PERMANOVA, R2 = 0.65, p < 0.001) (Figure S1) and eukaryotes (PERMANOVA, R2 = 0.70, p < 0.001) (Figure S2) at the genus level [30]. Linear discriminant analysis with effect size (LEfSe) analysis [31] revealed that the phyla Actinobacteriota, Bacteriodota, Chloroflexota, and Myxococcota were significantly enriched in biofilms grown in 30–80% WW (Fig. 1). These results align with previous 16S rRNA gene-based analysis showing that exposure to WWTP effluent alters microbiome composition [27] mirroring changes in microbial communities downstream of WWTPs relative to communities upstream [14]. However, these previous field studies did not include WW ultrafiltration treatment. We therefore compared UF WW and non-UF WW-grown biofilms and identified two bacterial phyla that were enriched in biofilms grown in UF WW, namely the Planctomycetota and Verrucomicrobiota (Fig. 1). These results show abiotic factors (UF treatment) specifically altered microbial community composition, e.g., by enriching taxa with known relevance for WWTP processes. For example, Planctomycetota is the only known phylum including taxa capable of anaerobic ammonia oxidation (anammox) for nitrogen removal from WW [25, 55, 56].

Fig. 1
figure 1

Taxonomic analysis of WW-grown stream biofilms. A Cladogram of the top 200 prokaryotic taxa based on relative abundance and shaded by significant differences in abundances between treatment groups. Node size corresponds to overall taxon abundance while node color indicates significant enrichment in a given experimental condition according to LefSE analysis. Branches are colored by phylum. B Barplot of relative taxonomic abundances of the significantly different phyla across treatments

The taxonomic composition of biofilms grown in UF WW was more similar to that of biofilms grown in 0% WW than to biofilms grown in non-UF WW, suggesting that biotic factors from WWTPs play a greater role in changes in the stream biofilm community than abiotic factors (Fig. 2, Figure S3). The metagenomic taxonomic analysis demonstrated that WW alters stream biofilm composition, supporting previous studies that some microorganisms from WW are able to colonize and modify biofilm community composition [27]. To leverage the functional diversity encoded in shotgun metagenomic data, we investigated whether the treatment also altered functional enzyme composition (see Methods). PCoA revealed that the functional enzyme composition (Figure S4) was also significantly associated with the experimental treatment (PERMANOVA, R2 = 0.58, p < 0.001).

Fig. 2
figure 2

Procrustes analysis shows a structural alignment between microbiome and resistome composition. Points indicate positions on PCoA axes 1 and 2 of prokaryotic composition (microbiome), while arrow tips indicate positions on PCoA axes 1 and 2 of resistome composition following Procrustes rotation. The correlation coefficient in a symmetric Procrustes rotation was 0.901, the sum of squares (m12 squared) was 0.1889, and the p value < 0.001

Resistome composition correlates with microbiome composition

We next examined how the total set of ARGs detected in each sample—the ‘resistome’—varied with experimental treatment (e.g., WW% and non-UF vs. UF). Resistome composition was significantly affected by experimental treatment (PERMANOVA, R2 = 0.60, p < 0.001), similar to taxonomic and functional enzyme composition (Figure S5). Biofilms exposed to 30% or 80% UF WW clustered with 0% WW samples as opposed to samples exposed to equivalent amounts of non-UF WW (Fig. 2). Although flow-through channels were assigned treatments using a randomized design aimed to avoid potential block effects or systematic influences, we still observed significant variability in resistome composition even within biofilm replicates. Such biological variability was also observed in variable biofilm biotransformation rates [11] and antibiotic resistance markers measured in other flume system studies [57].

To test for a relationship between the microbiome and resistome composition, we used Procrustes superimposition of resistome PCoA vectors onto prokaryotic microbiome PCoA vectors (Fig. 2). Procrustes analysis revealed that the composition of the resistome was structurally related to the composition of the prokaryotic microbiome. The low Procrustes sum of squares of 0.1889 and calculated correlation value of 0.9 in a symmetric Procrustes rotation (p value < 0.001) indicated an alignment between the resistome and microbiome. Similarly, functional enzyme composition was associated with the resistome composition (Figure S6, Procrustes corr. coeff. = 0.89, p value < 0.001). Alignments between the microbiome and resistome composition have been previously observed in microbial communities downstream of WWTPs [13, 14], but local environmental features at different sampling sites are potential confounding factors that may have determined both microbiome and resistome composition in the field. Here, we demonstrated a microbiome-resistome relationship in a controlled experiment with biofilms grown in a defined mixing ratio of stream water to WW and ultrafiltered WW. Generally, the rotated PCoA vectors for the resistome are closer to the prokaryotic microbiome PCoA vectors for biofilm samples than they are for water samples (Fig. 2). These results further suggest that changes in resistome composition in WW-impacted stream biofilms are influenced by biotic factors such as colonization by WWTP bacteria [27, 57].

WW ultrafiltration treatment has ARG class-specific effects on biofilm resistomes

A total of 570 ARGs across 41 different higher-level ARG classes were identified in the biofilm metagenomes using a short reads pipeline for ARG detection [34]. In biofilms, the abundances of ARGs conferring resistance to beta-lactams, fluoroquinolones, peptides, rifamycins, and tetracyclines, as well as multidrug resistance ARGs, differed significantly between conditions (Table S5). The normalized read counts of several, but not all, ARG classes increased with increasing WW%, as shown in Fig. 3. In contrast, the metagenomic abundances of ARG classes did not increase in biofilms grown in comparable percentages of UF WW with the exception of sulfonamide ARGs (Fig. 3). These UF treatment results therefore suggest biotic factors primarily alter the abundances of ARGs in WW-grown biofilms.

Fig. 3
figure 3

Significant ARG classes profiled using a metagenomic short read pipeline. Abundance (16S-normalized read counts) of the DeepARG hits assigned to ARG classes varied significantly with adjusted p values ≤ 0.001 (Kruskal–Wallis test + Benjamini–Hochberg multiple hypothesis correction) between different treatment conditions. MLS = macrolide-lincosamide-streptogramin

As a notable exception, sulfonamide ARGs exhibited increasing abundance with increasing WW% yet they displayed the highest abundance under 80% UF WW conditions, suggesting a contribution from abiotic factors such as nutrients or MPs in WW. The high abundance of sulfonamide ARGs in UF-grown biofilms may also be related to gene family-specific mobility mechanisms such as extracellular genetic material which is small enough to pass through the UF membrane. In support of this hypothesis, the sulfonamide ARGs sul1-4 are known to be highly mobilized in WWTP effluent microbiomes [58]. Although many ARGs are mobile, sul1 ARGs are known for their pervasive colocalization within class 1 integrons [44] enabling bacteria to acquire and swap sulfonamide resistance genes [59, 60]. Analysis of the assembled metagenomic contigs in our samples revealed that sul1 genes were frequently colocalized on the same metagenomic contigs as mobile genetic elements, including int1 integron integrases, as well as recombinases and transposases (Figure S7). Critically, the presence of intI1 genes in environmental samples has been suggested as a gene proxy for estimating the levels of anthropogenic pollution, including MPs, as proposed by Gillings et al. [60]. The high abundance of sulfonamide ARGs observed in UF WW-grown biofilms supports the relationship between int1 genes, sulfonamide ARGs typically associated with int1, and abiotic factors, including pollutants [59].

Associations between assembled ARG abundances and WW percentages

Metagenomic assembly-based approaches are likely to have lower sensitivity but higher specificity than short-read-based approaches in the detection of ARGs [61, 62]. Therefore, as a complementary analysis to the DeepARG short read-based detection methods, we used stringent criteria (see “Methods” section) to query our assembled metagenomes against the Structural Antibiotic Resistance Genes (SARG) database [43]. The query yielded 2240 hits that clustered into 102 ARG cluster representatives at the 70% amino acid sequence identity level (a threshold previously used [62] for ARGs identified from short-read metagenomic data). Analysis based on ARG class resulted in the following ARG cluster counts: 15 beta-lactam, 14 multidrug, 9 aminoglycoside, 5 tetracycline, 3 macrolide-lincosamide-streptogramin (MLS), 2 sulfonamide, 1 spectinomycin, and 1 vancomycin. In addition, 51 cluster representatives were classified by SARG into the “other” ARG category, a category which includes multiple antibiotic classes with few characterized representatives (e.g., bacitracin, triclosan, pyrazinamide).

To quantify relative ARG abundances across treatment conditions, we mapped metagenomic reads back to the 102 ARG cluster representatives. Ten ARG cluster representatives differed significantly (Kruskal–Wallis test, p < 0.05) in abundance between at least one experimental treatment condition (Table S6). In addition to two beta-lactam resistance ARGs, one each of the multidrug, MLS, and aminoglycoside ARG representatives displayed the trend of increasing abundance with increasing WW%. The fold change increase in metagenomic abundance was the greatest for two beta-lactam resistance ARG clusters (Class-A beta-lactamase gene, p value = 0.0054; and a Class-D beta-lactamase bla OXA-4 gene, p value = 0.0054; Fig. 4). In contrast, a third beta-lactam ARG, belonging to the Belgium extended-spectrum beta-lactamase (BEL) family, (p value = 0.0488) displayed high coverage exclusively in the 80% UF WW condition (Fig. 4). These results suggest that the abundance of this cluster of BEL family ARGs may be influenced by abiotic factors such as MPs, the effects of which may have been out-competed by biotic effects in the non-UF WW-grown biofilms. Enabled by our metagenomic assembly-based ARG detection approach, we examined genes localized on the same metagenomic contig as the BEL family beta-lactamase hit. The BEL family beta-lactamase is directly flanked by a Pseudomonas TniC resolvase. This TniC resolvase shares 100% sequence identity with a TniC resolvase (ADC80450.1) in integron/gene cassette systems previously shown to be influenced by abiotic factors, e.g., heavy metal concentrations [63]. Although the causality must be experimentally demonstrated, genetic mobility influenced by abiotic factors may contribute to the observed abundance of BEL family beta-lactamases in 80% UF WW-grown biofilms.

Fig. 4
figure 4

Significant ARG classes using a metagenomic assembly-based approach. The corrected normalized average coverage of ARG cluster representatives was identified using an assembly-based ARG detection method. ARG clusters with significant differences in abundance between at least one different treatment condition are displayed (adjusted p values < 0.05 using Kruskal–Wallis test + Benjamini–Hochberg multiple hypothesis correction)

Overall, our experimental design separated biotic and abiotic factors and shed light on the gene family-specific mechanisms underlying how WWTP effluent may alter the resistome composition of stream biofilms. Our combined results from short read mapping and assembly-based ARG analysis generally aligned with the exception of the sulfonamide ARGs, which did not show a significant difference in the assembly-based approach. This may be due to limitations in the metagenomic assembly and highlights the value of using two different, complementary approaches. For example, the significance of beta-lactamase ARGs identified via the short read-based approach was also captured via the assembly-based approach. Notably, a consensus for beta-lactamase abundance emerged from the comparison of our study to previous studies in Switzerland [13] and globally [9, 14, 64, 65]. Beta-lactamase ARGs consistently remain the most abundant class of ARGs detected in WWTP effluent and aquatic environments [7, 66,67,68,69]. Their prevalence can be attributed to a variety of factors, including the early discovery and continued widespread overuse of beta-lactam antibiotics in humans and animals [70, 71]. However, the unexpectedly high abundance of BEL family beta-lactamases exclusively in 80% UF WW samples (Fig. 4) prompted us to focus on beta-lactamases for further experimental investigation and analysis of abiotic factors.

Micropollutants as abiotic factors and their potential influence on ARG transfer

We examined MP concentrations as abiotic factors capable of influencing gene family-specific ARG abundances, such as the sul1-4 and BEL family beta-lactamase ARGs showing unexpectedly higher abundances in the UF WW compared to non-UF conditions (Figs. 3 and 4). Based on previous chemical analysis by liquid chromatography-tandem mass spectrometry [11], the MP concentrations measured in the flume biofilm samples increased with increasing WW% (Table S4), but the concentrations of MPs extracted from biofilms did not significantly differ between the UF and non-UF WW. Many of the antibiotics measured, including sulfamethoxazole and sulfapyridine (Table S4), were present at concentrations below or near the limit of quantification in WW-grown biofilms. Among all antibiotics detected, the macrolide clarithromycin exhibited the highest concentrations in 80% WW and 80% UF WW-grown biofilms (45–95 ng per mg biofilm, Table S4). In line with the general MP trends, clarithromycin concentrations increased with WW% and did not vary significantly between UF and non-UF treatments. In contrast, the metagenomic abundances of macrolide-lincosamide-streptogramin (MLS) ARGs, which typically confer resistance to clarithromycin, were significantly lower in UF than in non-UF conditions (Fig. 4). This represents a case of an antibiotic-ARG pair where ARG abundances did not correlate with the expected abiotic factor (clarithromycin). However, this does not take into account the potential effects of other MPs which are not classified as antibiotics but nonetheless influence ARG abundance [3].

A growing body of evidence indicates that many anthropogenic pollutants that do not have known antimicrobial bioactivities nonetheless influence ARG transfer and proliferation [3, 72,73,74,75,76,77]. However, there is generally a lack of knowledge on the mechanisms underlying how MPs without antibiotic activities act as selective agents for ARG abundance. In our study, the 75 different MPs measured [11] included artificial sweeteners, non-antibiotic pharmaceuticals, corrosion inhibitors, and pesticides (Table S4). Compared with the other compound classes, artificial sweeteners were the most abundant MPs, with concentrations in the µg/L range detected in WW samples [11] (Table S4). Notably, the artificial sweeteners measured in the biofilms (acesulfame K, cyclamate, and saccharin) all contain sulfonamide moieties which are chemically analogous to the moieties found in sulfonamide antibiotics (Figure S8). Sulfonamides were among the few ARG classes with significantly greater abundances in UF WW-grown biofilms than in non-UF WW-grown biofilms (Fig. 3). These observations led us to speculate that abiotic factors, such as sulfonamide-containing MPs, may play a role in promoting the abundance of such gene family-specific ARGs. In support of this, several experimental studies have shown that artificial sweeteners promote high rates of horizontal gene transfer and natural uptake of ARGs in bacteria [76, 78]. Further work has demonstrated that acesulfame-K and other artificial sweeteners directly enhance the rate and quantity of horizontal plasmid-mediated transfers of sul genes [79]. Interestingly, a plasmid-encoded microbial hydrolase capable of cleaving the sulfonamide moiety of acesulfame-K was recently identified and shown to have a metallo-beta-lactamase structural fold [80]. Taken together, these findings prompted us to conduct a deeper functional investigation of beta-lactamase fold genes whose abundance increased with WW%.

Planctomycetota beta-lactamases display differential abundances in response to sulfonamide-containing MPs and WWTP effluent

To investigate the relationship between ARG abundances and sulfonamide-containing MPs, we consulted previous reports of enrichment studies of activated sludge communities grown on sulfonamide-containing MPs including artificial sweeteners [81, 82]. Multiple independent enrichment studies on acesulfame-K [82], saccharin [81], and cyclamate [81] reported significant increases in the relative abundance of Planctomycetota and particularly the family Phycisphaeraceae. Previously, Desiante et al. [11] also identified significant differences in the relative abundances of Planctomycetota and Phycisphaeraceae in biofilms grown in 80% WW and 80% UF WW in two independent flume biofilm experiments (log2-fold change > 4; p < 0.05). We therefore conducted differential abundance analysis and observed statistically significant differences in the Planctomycetota relative abundances between WW and UF WW treatments (Fig. 5A). Taxa assigned to the phylum Planctomycetota constitute up to 2.1% of the relative abundance of our biofilm, WW, and stream water metagenomes, confirming previous reports of the Planctomycetota as a numerically important constituent of both freshwater and wastewater microbial communities [25, 26, 56]. Many Phycisphaeraceae are algae-associated with the metabolic potential to impact the fitness of their algal hosts [83], but the current unculturability of most planctomycetes has limited experimental validation of Phycisphaeraceae ARGs. To address this knowledge gap, we chose to target Phycisphaeraceae ARGs for experimental characterization.

Fig. 5
figure 5

Functional characterization of an active Planctomycetota beta-lactamase. A Volcano plot displaying the log2-fold change of differentially abundant taxa between 30% WW UF vs. 30% WW. B Counts of hits showing significant homology to the putative beta-lactamase MBX3358097.1 in biofilm metagenomes. C Purified enzyme assays demonstrated beta-lactamase activity of MBX3358097.1 with the colorimetric cephalosporin analog, nitrocefin, compared to boiled enzyme and buffer-only controls. D Two views of a high-confidence AlphaFold2 model (mean pLDDT = 96.5) of MBX3358097.1 structurally aligned with its closest structural homolog in the Protein Data Bank (PDB100), as identified by FoldSeek [84] in November 2023, the Bacillus cereus beta-lactamase (PDB ID: 3KNS)

Due to known challenges in planctomycete cultivation [56], we used heterologous expression techniques to functionally validate two Phycisphaeraceae ARGs. We focused specifically on Phycisphaeraceae beta-lactamases for four reasons: (1) the beta-lactamases were identified in our ARG analysis as the most abundant and significantly different class of ARGs affected by the percentage of wastewater effluent, (2) beta-lactamase fold enzymes were also previously identified as plasmid-encoded enzymes active in the biotransformation of the sulfonamide MP acesulfame-K [80], revealing the promiscuity of a versatile metallo-beta-lactamase fold capable of cleaving both sulfonamide- and beta-lactam moieties, and (3) previous evidence for sulfonamide cleavage in the sweetener cyclamate was shown by metal-dependent hydrolases but their sequences were not determined [85]. Finally, (4) full-length homologs with greater than 50% amino acid identity to the sulfonamide resistance genes sul1-4 were not detected in the Phycisphaeraceae. Since classical sulfonamide resistance genes were lacking, we hypothesized that metallo-beta-lactamase-fold enzymes may instead play an overlooked role in sulfonamide cleavage in the Phycisphaeraceae.

Functional characterization of a new beta-lactamase ARG from the Planctomycetota

To identify metallo-beta-lactamase fold enzymes assigned to the Phycisphaeraceae, we queried the 75 metagenomes sequenced in this study as well as additional publicly-available metagenomes from previous sulfonamide-containing MP enrichment studies [81, 82]. Through bioinformatic prioritization (Methods), we selected two full-length metallo-beta-lactamase fold protein hits from sulfonamide-enriched Phycisphaeraceae MAGs for functional characterization [82]. The abundances of both proteins increased with increasing WW% (Fig. 5B and Figure S9). We heterologously expressed each of the Phycisphaeraceae genes in E. coli. One protein did not express in E. coli and despite troubleshooting was not amenable to protein production or purification with nickel-affinity chromatography. The second protein (MBX3358097.1) was soluble and purified to homogeneity with an average yield of 500 µg/mL protein (Figure S10) therefore we proceeded with the functional analysis of this protein. We tested the purified Phycisphaeraceae protein for beta-lactamase activity and observed enzyme-mediated cleavage of the cephalosporin analog, nitrocefin relative to boiled enzyme controls (Fig. 5C). We next tested its substrate specificity for sulfonamide-containing compounds. By high-performance liquid chromatography analysis and comparison to authentic standards, we determined that the Phycisphaeraceae enzyme did not biotransform the sulfonamide-containing MPs acesulfame-K, saccharin, or cyclamate, nor did it transform sulfamethoxazole, thus disproving our hypothesis for its role in sulfonamide cleavage. Nonetheless, this analysis resulted in the first biochemical characterization of a new, active beta-lactamase from the Phycisphaeraceae.

MBX3358097.1 displayed a maximum 30% amino acid identity and between 20–50% query coverage with proteins of known function in the UniProtKB-SwissProt [86] and the Beta-Lactamase database (BLDB) [87]. These low sequence similarities limit reliable functional annotation and beta-lactamase classification based on homology alone. Therefore, we generated an AlphaFold2 model (Fig. 5D, model available on GitHub) for the structural alignment and identification of beta-lactamase subfamily-specific motifs [88, 89]. MBX3358097.1 encodes a zinc-binding motif 2 (HxHxDH), which has a terminal histidine characteristic of the subclass B3 metallo-beta-lactamases. However, in another known motif 4, this enzyme contains the residues “GD”, which are more characteristic of related glyoxalase II enzymes rather than B1–B3 family beta-lactamases [88, 89]. These results are suggestive of promiscuous activity or additional enzymatic functions of MBX3358097.1 to investigate which are beyond the scope of this current study. Nonetheless, this analysis highlights the benefits of structural modeling of metagenomic sequences [90, 91], particularly in cases where low sequence identity prevents reliable functional annotation.

There are few previous reports of ARGs from the phylum Planctomycetota [92]. Planctomycetes were previously reported to have an unusual proteinaceous cell wall composition devoid of peptidoglycan [93]; however, later findings demonstrated that planctomycetes do indeed produce peptidoglycan and are thus susceptible to cell-wall targeting antibiotics including beta-lactams [94]. In fact, a study of six culturable planctomycetes revealed that all of the strains tested were resistant to beta-lactam antibiotics although beta-lactamases could not be detected [95]. Importantly, MBX3358097.1 was also not detected by standard computational methods for detecting ARGs, including queries to the SARG [43] and CARD [96] databases due to its low sequence identity with known ARGs. This work therefore highlights the size of the knowledge gap and invites further inquiry into the role of sequence-diverse ARGs from uncultivated microorganisms and their contributions to environmental ARG reservoirs.

Conclusions

Our results from a controlled flume system showed that WWTP effluent altered both the microbiome and resistome composition of stream biofilms. Through ultrafiltration treatment of WWTP effluent, we could attribute changes in ARG abundances primarily to biotic rather than abiotic factors. However, the abundance of several gene family-specific ARGs, most notably sulfonamide and beta-lactamase resistance genes, suggested an influence by abiotic factors. We biochemically characterized the beta-lactamase activity of a new Phycisphaeraceae protein with a low sequence identity to known ARGs. This provided proof-of-principle of a functional ARG from an uncultivated microorganism that was overlooked by standard computational ARG detection methods. This further highlights the need to characterize new ARGs from understudied taxa and to create a more comprehensive ARG catalog for querying complex microbiomes. Overall, our study exemplified a workflow to disentangle antibiotic resistance determinants in aquatic environments.

Availability of data and materials

Metagenomic sequencing data generated during the current study have been deposited in the European Nucleotide Archive with the study accession code: PRJNA1008123. Scripts and additional data used in analysis and figures are available at: https://github.com/MSM-group/AMR-ecoimpact-paper.

Abbreviations

ARGs:

Antibiotic resistance genes

DNA:

Deoxyribonucleic acid

MPs:

Micropollutants

ORFs:

Open reading frames

PCoA:

Principal coordinate analysis

PCR:

Polymerase chain reaction

rRNA:

Ribosomal ribonucleic acid

SARG:

Structured antibiotic resistance genes database

UF:

Ultrafiltration

WW:

Wastewater effluent

WWTPs:

Wastewater treatment plants

References

  1. Marti E, Variatza E, Balcazar JL. The role of aquatic ecosystems as reservoirs of antibiotic resistance. Trends Microbiol. 2014;22:36–41.

    Article  CAS  PubMed  Google Scholar 

  2. Singer AC, Shaw H, Rhodes V, Hart A. Review of Antimicrobial Resistance in the Environment and Its Relevance to Environmental Regulators. Front Microbiol. 2016;7:1728.

    Article  PubMed  PubMed Central  Google Scholar 

  3. Di Cesare A, Sabatino R, Sbaffi T, Fontaneto D, Brambilla D, Beghi A, et al. Anthropogenic pollution drives the bacterial resistome in a complex freshwater ecosystem. Chemosphere. 2023;331:138800.

    Article  PubMed  Google Scholar 

  4. Robins K, McCann CM, Zhou XY, Su JQ, Cooke M, Knapp CW, Graham DW. Bioavailability of potentially toxic elements influences antibiotic resistance gene and mobile genetic element abundances in urban and rural soils. Sci Total Environ. 2022;847:157512.

  5. Shallcross LJ, Davies DSC. Antibiotic overuse: a key driver of antimicrobial resistance. Br J Gen Pract. 2014;64:604–5.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Alexander J, Hembach N, Schwartz T. Evaluation of antibiotic resistance dissemination by wastewater treatment plant effluents with different catchment areas in Germany. Sci Rep. 2020;10:8952.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Rodriguez-Mozaz S, Chamorro S, Marti E, Huerta B, Gros M, Sànchez-Melsió A, et al. Occurrence of antibiotics and antibiotic resistance genes in hospital and urban wastewaters and their impact on the receiving river. Water Res. 2015;69:234–42.

    Article  CAS  PubMed  Google Scholar 

  8. Rogowska J, Cieszynska-Semenowicz M, Ratajczyk W, Wolska L. Micropollutants in treated wastewater. Ambio. 2020;49:487–503.

    Article  PubMed  Google Scholar 

  9. Mukherjee M, Laird E, Gentry TJ, Brooks JP, Karthikeyan R. Increased Antimicrobial and Multidrug Resistance Downstream of Wastewater Treatment Plants in an Urban Watershed. Front Microbiol. 2021;12:657353.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Carles L, Wullschleger S, Joss A, Eggen RIL, Schirmer K, Schuwirth N, et al. Impact of wastewater on the microbial diversity of periphyton and its tolerance to micropollutants in an engineered flow-through channel system. Water Res. 2021;203:117486.

    Article  CAS  PubMed  Google Scholar 

  11. Desiante WL, Carles L, Wullschleger S, Joss A, Stamm C, Fenner K. Wastewater microorganisms impact the micropollutant biotransformation potential of natural stream biofilms. Water Res. 2022;217:118413.

    Article  CAS  PubMed  Google Scholar 

  12. Perveen S, Pablos C, Reynolds K, Stanley S, Marugán J. Growth and prevalence of antibiotic-resistant bacteria in microplastic biofilm from wastewater treatment plant effluents. Sci Total Environ. 2023;856:159024.

    Article  CAS  PubMed  Google Scholar 

  13. Lee J, Ju F, Beck K, Bürgmann H. Differential effects of wastewater treatment plant effluents on the antibiotic resistomes of diverse river habitats. ISME J. 2023;17:1993–2002. Available from: https://doi.org/10.1038/s41396-023-01506-w.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Matviichuk O, Mondamert L, Geffroy C, Gaschet M, Dagot C, Labanowski J. River Biofilms Microbiome and Resistome Responses to Wastewater Treatment Plant Effluents Containing Antibiotics. Front Microbiol. 2022;13:795206.

    Article  PubMed  PubMed Central  Google Scholar 

  15. Depetris A, Tagliavini G, Peter H, Kühl M, Holzner M, Battin TJ. Biophysical properties at patch scale shape the metabolism of biofilm landscapes. NPJ Biofilms Microbiomes. 2022;8:5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Romero F, Acuña V, Font C, Freixa A, Sabater S. Effects of multiple stressors on river biofilms depend on the time scale. Sci Rep. 2019;9:15810.

    Article  PubMed  PubMed Central  Google Scholar 

  17. Haenelt S, Richnow H-H, Müller JA, Musat N. Antibiotic resistance indicator genes in biofilm and planktonic microbial communities after wastewater discharge. Front Microbiol. 2023;14:1252870.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Ju F, Lee J, Beck K, Zhang G, Gekenidis M-T, Hummerjohann J, et al. Phenotypic Metagenomics tracks Wastewater-Associated Clinically Important Beta-lactam Resistant Bacteria Invading River Habitats. Research Square. 2022. Available from: https://www.researchsquare.com/article/rs-1589365/latest.pdf. Accessed 23 Nov 2023.

  19. Reddington K, Eccles D, O’Grady J, Drown DM, Hansen LH, Nielsen TK, et al. Metagenomic analysis of planktonic riverine microbial consortia using nanopore sequencing reveals insight into river microbe taxonomy and function. Gigascience. 2020;9. Available from: https://doi.org/10.1093/gigascience/giaa053.

  20. Brumfield KD, Huq A, Colwell RR, Olds JL, Leddy MB. Microbial resolution of whole genome shotgun and 16S amplicon metagenomic sequencing using publicly available NEON data. PLoS ONE. 2020;15:e0228899.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Tremblay J, Singh K, Fern A, Kirton ES, He S, Woyke T, et al. Primer and platform effects on 16S rRNA tag sequencing. Front Microbiol. 2015;6:771.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Klindworth A, Pruesse E, Schweer T, Peplies J, Quast C, Horn M, et al. Evaluation of general 16S ribosomal RNA gene PCR primers for classical and next-generation sequencing-based diversity studies. Nucleic Acids Res. 2013;41:e1.

    Article  CAS  PubMed  Google Scholar 

  23. Oyarzúa P, Bovio-Winkler P, Etchebehere C, Suárez-Ojeda ME. Microbial communities in an anammox reactor treating municipal wastewater at mainstream conditions: Practical implications of different molecular approaches. J Environ Chem Eng. 2021;9:106622.

    Article  Google Scholar 

  24. Zhang Y, Zhao Z, Xu H, Wang L, Liu R, Jia X. Fate of antibiotic resistance genes and bacteria in a coupled water-processing system with wastewater treatment plants and constructed wetlands in coastal eco-industrial parks. Ecotoxicol Environ Saf. 2023;252:114606.

    Article  CAS  PubMed  Google Scholar 

  25. Andrei A-Ş, Salcher MM, Mehrshad M, Rychtecký P, Znachor P, Ghai R. Niche-directed evolution modulates genome architecture in freshwater Planctomycetes. ISME J. 2019;13:1056–71.

    Article  PubMed  PubMed Central  Google Scholar 

  26. Cai X, Yao L, Sheng Q, Jiang L, Wang T, Dahlgren RA, et al. Influence of a biofilm bioreactor on water quality and microbial communities in a hypereutrophic urban river. Environ Technol. 2021;42:1452–60.

    Article  CAS  PubMed  Google Scholar 

  27. Carles L, Wullschleger S, Joss A, Eggen RIL, Schirmer K, Schuwirth N, Stamm C, Tlili A. Wastewater microorganisms impact microbial diversity and important ecological functions of stream periphyton. Water Res. 2022;225:119119.

    Article  CAS  PubMed  Google Scholar 

  28. Li D, Liu C-M, Luo R, Sadakane K, Lam T-W. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics. 2015;31:1674–6.

    Article  CAS  PubMed  Google Scholar 

  29. Hyatt D, Chen G-L, Locascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 2010;11:119.

    Article  PubMed  PubMed Central  Google Scholar 

  30. Salazar G, Ruscheweyh H-J, Hildebrand F, Acinas SG, Sunagawa S. mTAGs: taxonomic profiling using degenerate consensus reference sequences of ribosomal RNA genes. Bioinformatics. 2021;38:270–2.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett W, Huttenhower C. Metagenomic biomarker discovery and explanation. Genome Biol. 2011;12:R60.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Mansfeldt C, Achermann S, Men Y, Walser J-C, Villez K, Joss A, et al. Microbial residence time is a controlling parameter of the taxonomic composition and functional profile of microbial communities. ISME J. 2019;13:1589–601.

    Article  PubMed  PubMed Central  Google Scholar 

  33. Buchfink B, Reuter K, Drost H-G. Sensitive protein alignments at tree-of-life scale using DIAMOND. Nat Methods. 2021;18:366–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Arango-Argoty G, Garner E, Pruden A, Heath LS, Vikesland P, Zhang L. DeepARG: a deep learning approach for predicting antibiotic resistance genes from metagenomic data. Microbiome. 2018;6:23.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Wickham H, Averick M, Bryan J, Chang W, McGowan L, François R, et al. Welcome to the tidyverse. J Open Source Softw. 2019;4:1686.

    Article  Google Scholar 

  36. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010;11:R25.

    Article  PubMed  PubMed Central  Google Scholar 

  37. Pereira MB, Wallroth M, Jonsson V, Kristiansson E. Comparison of normalization methods for the analysis of metagenomic gene abundance data. BMC Genomics. 2018;19:274.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

    Article  CAS  PubMed  Google Scholar 

  39. Oksanen J, Kindt R, Legendre P, O’Hara B, Stevens MHH, Oksanen MJ, et al. The vegan package. Community Ecol Package. 2007;10:719.

    Google Scholar 

  40. Paradis E, Claude J, Strimmer K. APE: Analyses of Phylogenetics and Evolution in R language. Bioinformatics. 2004;20:289–90.

    Article  CAS  PubMed  Google Scholar 

  41. Benjamini Y, Hochberg Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J R Stat Soc. 1995;57:289–300.

    Article  Google Scholar 

  42. McGinnis S, Madden TL. BLAST: at the core of a powerful and diverse set of sequence analysis tools. Nucleic Acids Res. 2004;32:W20–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Yin X, Jiang X-T, Chai B, Li L, Yang Y, Cole JR, et al. ARGs-OAP v2.0 with an expanded SARG database and Hidden Markov Models for enhancement characterization and quantification of antibiotic resistance genes in environmental metagenomes. Bioinformatics. 2018;34:2263–70.

    Article  CAS  PubMed  Google Scholar 

  44. Lee J, Ju F, Maile-Moskowitz A, Beck K, Maccagnan A, McArdell CS, et al. Unraveling the riverine antibiotic resistome: The downstream fate of anthropogenic inputs. Water Res. 2021;197:117050.

    Article  CAS  PubMed  Google Scholar 

  45. Fu L, Niu B, Zhu Z, Wu S, Li W. CD-HIT: accelerated for clustering the next-generation sequencing data. Bioinformatics. 2012;28:3150–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Li W, Godzik A. Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences. Bioinformatics. 2006;22:1658–9.

    Article  CAS  PubMed  Google Scholar 

  47. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012;9:357–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Albertsen M, Hugenholtz P, Skarshewski A, Nielsen KL, Tyson GW, Nielsen PH. Genome sequences of rare, uncultured bacteria obtained by differential coverage binning of multiple metagenomes. Nat Biotechnol. 2013;31:533–8.

    Article  CAS  PubMed  Google Scholar 

  49. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. 2009;10:421.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Mirdita M, Schütze K, Moriwaki Y, Heo L, Ovchinnikov S, Steinegger M. ColabFold: making protein folding accessible to all. Nat Methods. 2022;19:679–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Teufel F, AlmagroArmenteros JJ, Johansen AR, Gíslason MH, Pihl SI, Tsirigos KD, et al. SignalP 6.0 predicts all five types of signal peptides using protein language models. Nat Biotechnol. 2022;40:1023–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Bradford MM. A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding. Analytical Biochem. 1976;72:248–54.

    Article  CAS  Google Scholar 

  53. O’Callaghan CH, Morris A, Kirby SM, Shingler AH. Novel method for detection of beta-lactamases by using a chromogenic cephalosporin substrate. Antimicrob Agents Chemother. 1972;1:283–8.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Ryu BH, Ngo TD, Yoo W, Lee S, Kim B-Y, Lee E, et al. Biochemical and Structural Analysis of a Novel Esterase from Caulobacter crescentus related to Penicillin-Binding Protein (PBP). Sci Rep. 2016;6:37978.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Wu P, Chen J, Garlapati VK, Zhang X, Wani Victor Jenario F, Li X, et al. Novel insights into Anammox-based processes: A critical review. Chem Eng J. 2022;444:136534.

    Article  CAS  Google Scholar 

  56. Fuerst JA. Planctomycetes: Cell Structure, Origins and Biology. New York: Springer Science & Business Media; 2013. https://link.springer.com/content/pdf/10.1007/978-1-62703-502-6.pdf.

  57. Gionchetta G, Snead D, Semerad S, Beck K, Pruden A, Bürgmann H. Dynamics of antibiotic resistance markers and Escherichia coli invasion in riverine heterotrophic biofilms facing increasing heat and flow stagnation. Sci Total Environ. 2023;893:164658.

    Article  CAS  PubMed  Google Scholar 

  58. Ju F, Beck K, Yin X, Maccagnan A, McArdell CS, Singer HP, et al. Wastewater treatment plant resistomes are shaped by bacterial composition, genetic exchange, and upregulated expression in the effluent microbiomes. ISME J. 2019;13:346–60.

    Article  PubMed  Google Scholar 

  59. Gillings MR, Boucher Y, Labbate M, Holmes A, Krishnan S, Holley M, et al. The evolution of class 1 integrons and the rise of antibiotic resistance. J Bacteriol. 2008;190:5095–100.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Gillings MR, Gaze WH, Pruden A, Smalla K, Tiedje JM, Zhu Y-G. Using the class 1 integron-integrase gene as a proxy for anthropogenic pollution. ISME J. 2015;9:1269–79.

    Article  CAS  PubMed  Google Scholar 

  61. Tamames J, Cobo-Simón M, Puente-Sánchez F. Assessing the performance of different approaches for functional and taxonomic annotation of metagenomes. BMC Genomics. 2019;20:960.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Berglund F, Österlund T, Boulund F, Marathe NP, Larsson DGJ, Kristiansson E. Identification and reconstruction of novel antibiotic resistance genes from metagenomes. Microbiome. 2019;7(1):52. https://doi.org/10.1186/s40168-019-0670-1.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Rosewarne CP, Pettigrove V, Stokes HW, Parsons YM. Class 1 integrons in benthic bacterial communities: abundance, association with Tn402-like transposition modules and evidence for coselection with heavy-metal resistance. FEMS Microbiol Ecol. 2010;72:35–46.

    Article  CAS  PubMed  Google Scholar 

  64. Thornton CN, Tanner WD, VanDerslice JA, Brazelton WJ. Localized effect of treated wastewater effluent on the resistome of an urban watershed. Gigascience. 2020;9. Available from: https://doi.org/10.1093/gigascience/giaa125.

  65. Czekalski N, Sigdel R, Birtel J, Matthews B, Bürgmann H. Does human activity impact the natural antibiotic resistance background? Abundance of antibiotic resistance genes in 21 Swiss lakes. Environ Int. 2015;81:45–55.

    Article  CAS  PubMed  Google Scholar 

  66. Laht M, Karkman A, Voolaid V, Ritz C, Tenson T, Virta M, et al. Abundances of tetracycline, sulphonamide and beta-lactam antibiotic resistance genes in conventional wastewater treatment plants (WWTPs) with different waste load. PLoS ONE. 2014;9:e103705.

    Article  PubMed  PubMed Central  Google Scholar 

  67. Stoll C, Sidhu JPS, Tiehm A, Toze S. Prevalence of clinically relevant antibiotic resistance genes in surface water samples collected from Germany and Australia. Environ Sci Technol. 2012;46:9716–26.

    Article  CAS  PubMed  Google Scholar 

  68. Osińska A, Korzeniewska E, Harnisz M, Felis E, Bajkacz S, Jachimowicz P, et al. Small-scale wastewater treatment plants as a source of the dissemination of antibiotic resistance genes in the aquatic environment. J Hazard Mater. 2020;381:121221.

    Article  PubMed  Google Scholar 

  69. Fadare FT, Okoh AI. The abundance of genes encoding ESBL, pAmpC and non-β-lactam resistance in multidrug-resistant Enterobacteriaceae recovered from wastewater effluents. Front Environ Sci Eng China. 2021;9. Available from: https://www.frontiersin.org/articles/10.3389/fenvs.2021.711950/full.

  70. Sta Ana KM, Madriaga J, Espino MP. β-Lactam antibiotics and antibiotic resistance in Asian lakes and rivers: An overview of contamination, sources and detection methods. Environ Pollut. 2021;275:116624.

    Article  CAS  PubMed  Google Scholar 

  71. Agarwal V, Tiwari A, Varadwaj P. An Extensive Review on β-lactamase Enzymes and their Inhibitors. Curr Med Chem. 2023;30:783–808.

    Article  CAS  PubMed  Google Scholar 

  72. Wang Y, Lu J, Zhang S, Li J, Mao L, Yuan Z, et al. Non-antibiotic pharmaceuticals promote the transmission of multidrug resistance plasmids through intra- and intergenera conjugation. ISME J. 2021;15:2493–508.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Kimbell LK, LaMartina EL, Kohls S, Wang Y, Newton RJ, McNamara PJ. Impact of corrosion inhibitors on antibiotic resistance, metal resistance, and microbial communities in drinking water. mSphere. 2023;8:e0030723.

  74. Alderton I, Palmer BR, Heinemann JA, Pattis I, Weaver L, Gutiérrez-Ginés MJ, et al. The role of emerging organic contaminants in the development of antimicrobial resistance. Emerg Contaminants. 2021;7:160–71.

    Article  CAS  Google Scholar 

  75. Qiu D, Ke M, Zhang Q, Zhang F, Lu T, Sun L, et al. Response of microbial antibiotic resistance to pesticides: An emerging health threat. Sci Total Environ. 2022;850:158057.

    Article  CAS  PubMed  Google Scholar 

  76. Yu Z, Guo J. Non-caloric artificial sweeteners exhibit antimicrobial activity against bacteria and promote bacterial evolution of antibiotic tolerance. J Hazard Mater. 2022;433:128840.

    Article  CAS  PubMed  Google Scholar 

  77. Zhang H, Song J, Zheng Z, Li T, Shi N, Han Y, et al. Fungicide exposure accelerated horizontal transfer of antibiotic resistance genes via plasmid-mediated conjugation. Water Res. 2023;233:119789.

    Article  CAS  PubMed  Google Scholar 

  78. Yu Z, Wang Y, Henderson IR, Guo J. Artificial sweeteners stimulate horizontal transfer of extracellular antibiotic resistance genes through natural transformation. ISME J. 2022;16:543–54.

    Article  CAS  PubMed  Google Scholar 

  79. Li Z, Gao J, Guo Y, Cui Y, Wang Y, Duan W, et al. Enhancement of antibiotic resistance dissemination by artificial sweetener acesulfame potassium: Insights from cell membrane, enzyme, energy supply and transcriptomics. J Hazard Mater. 2022;422:126942.

    Article  CAS  PubMed  Google Scholar 

  80. Bonatelli ML, Rohwerder T, Popp D, Liu Y, Akay C, Schultz C, et al. Recently evolved combination of unique sulfatase and amidase genes enables bacterial degradation of the wastewater micropollutant acesulfame worldwide. Front Microbiol. 2023;14:1223838.

    Article  PubMed  PubMed Central  Google Scholar 

  81. Deng Y, Wang Y, Xia Y, Zhang AN, Zhao Y, Zhang T. Genomic resolution of bacterial populations in saccharin and cyclamate degradation. Sci Total Environ. 2019;658:357–66.

    Article  CAS  PubMed  Google Scholar 

  82. Huang Y, Deng Y, Law JC-F, Yang Y, Ding J, Leung KS-Y, et al. Acesulfame aerobic biodegradation by enriched consortia and Chelatococcus spp: Kinetics, transformation products, and genomic characterization. Water Res. 2021;202:117454.

    Article  CAS  PubMed  Google Scholar 

  83. Rambo IM, Dombrowski N, Constant L, Erdner D, Baker BJ. Metabolic relationships of uncultured bacteria associated with the microalgae Gambierdiscus. Environ Microbiol. 2020;22:1764–83.

    Article  CAS  PubMed  Google Scholar 

  84. van Kempen M, Kim SS, Tumescheit C, Mirdita M, Lee J, Gilchrist CLM, et al. Fast and accurate protein structure search with Foldseek. Nat Biotechnol. 2023. Available from: https://doi.org/10.1038/s41587-023-01773-0.

  85. Nimura T, Tokieda T, Yamaha T. Partial purification and some properties of cyclamate sulfamatase. J Biochem. 1974;75:407–17.

    Article  CAS  PubMed  Google Scholar 

  86. UniProt Consortium. UniProt: the universal protein knowledgebase in 2021. Nucleic Acids Res. 2021;49:D480–9.

    Article  Google Scholar 

  87. Naas T, Oueslati S, Bonnin RA, Dabos ML, Zavala A, Dortet L, et al. Beta-lactamase database (BLDB) - structure and function. J Enzyme Inhib Med Chem. 2017;32:917–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  88. González JM. Visualizing the superfamily of metallo-β-lactamases through sequence similarity network neighborhood connectivity analysis. Heliyon. 2021;7:e05867.

    Article  PubMed  PubMed Central  Google Scholar 

  89. Denakpo E, Arlet G, Philippon A, Iorga BI. Metallo-β-lactamases. In Metalloenzymes. London: Academic Press, Elsevier; 2024. p. 157–84.

  90. Robinson SL. Structure-guided metagenome mining to tap microbial functional diversity. Curr Opin Microbiol. 2023;76:102382.

    Article  CAS  PubMed  Google Scholar 

  91. Ruppé E, Ghozlane A, Tap J, Pons N, Alvarez A-S, Maziers N, et al. Prediction of the intestinal resistome by a three-dimensional structure-based method. Nat Microbiol. 2019;4:112–23.

    Article  PubMed  Google Scholar 

  92. Ivanova AA, Miroshnikov KK, Oshkin IY. Exploring Antibiotic Susceptibility, Resistome and Mobilome Structure of Planctomycetes from Gemmataceae Family. Sustain Sci Pract Policy. 2021;13:5031.

    CAS  Google Scholar 

  93. König E, Schlesner H, Hirsch P. Cell wall studies on budding bacteria of the Planctomyces/Pasteuria group and on a Prosthecomicrobium sp. Arch Microbiol. 1984;138:200–5.

    Article  Google Scholar 

  94. Jeske O, Schüler M, Schumann P, Schneider A, Boedeker C, Jogler M, et al. Planctomycetes do possess a peptidoglycan cell wall. Nat Commun. 2015;6:7116.

    Article  CAS  PubMed  Google Scholar 

  95. Cayrou C, Raoult D, Drancourt M. Broad-spectrum antibiotic resistance of Planctomycetes organisms determined by Etest. J Antimicrob Chemother. 2010;65:2119–22.

    Article  CAS  PubMed  Google Scholar 

  96. Alcock BP, Huynh W, Chalil R, Smith KW, Raphenya AR, Wlodarski MA, et al. CARD 2023: expanded curation, support for machine learning, and resistome prediction at the Comprehensive Antibiotic Resistance Database. Nucleic Acids Res. 2023;51:D690–9.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

We thank Dr. Bogdan Iorga for his export support in beta-lactamase classification and bioinformatic analysis. Dr. Robert Niederdorfer is acknowledged for his support and expertise in the metagenomic assembly through computational resources and in collaboration with the Genetic Diversity Centre (GDC), ETH Zurich. We acknowledge Niklas Ferenc Trottmann for his training and support for the cloning, expression, and protein purification. We thank the Eawag Partnership Program together with the Institute of Water Education (IHE Delft) under the auspices of UNESCO for supporting the exchange of Mustafa Attrah to Eawag. We thank Elia Ceppi for assistance in interpretation of chemical micropollutant data and Thierry Marti for support in biological interpretation. We further acknowledge Stuart Dennis for his computational infrastructure support at Eawag. We acknowledge the EcoImpact teams 1.0 and 2.0 for fruitful discussions as well as conceptualizing and running the flume experiments which provided the samples necessary for this study.

Funding

S.L.R. acknowledges support from the Swiss National Science Foundation (grant no. 501100001711-209124, project no. PZPGP2_209124), the Vontobel Foundation, and the Pierre Mercier Foundation. K.F. acknowledges the Swiss National Science Funding project no. 200021L_201006. H.B. and G.G. were supported by the ANTIVERSA project (BiodivERsA, European Union), Swiss National Science Foundation (Switzerland) grant no. 186531. Student funding for an advanced class at IHE Delft Institute of Water Education was provided to M.A. through IHE Delft’s DUPC2 programme project ‘Supporting integrated and sustainable water management in Iraq through capacity development and research’ grant no. 109070.

Author information

Authors and Affiliations

Authors

Contributions

M.A. contributed to the conceptualization, methodology, investigation, writing—editing and revisions of the manuscript, visualization, and funding acquisition. M.S. contributed to conceptualization, methodology, visualization, investigation, data curation, writing—editing and revisions of the manuscript. M.E. contributed to the methodology, investigation, and revision and editing of the manuscript. G.G. contributed to methodology and editing and revision of the manuscript. H.B. contributed to the conceptualization, and editing and revision of the manuscript. P.N.L.L. contributed to the revision and editing of the manuscript and project administration. K.F. contributed to the conceptualization, revision and editing of the manuscript and funding acquisition. J.v.d.V. contributed to the methodology, writing- editing and revision of the manuscript and project administration. S.L.R contributed to conceptualization, methodology, Investigation, writing- editing and revision of the manuscript, visualization, project administration and funding acquisition.

Corresponding author

Correspondence to Serina L. Robinson.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

None.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

40168_2024_1879_MOESM1_ESM.pdf

Supplementary Material 1. Figure S1. Principal coordinate analysis (PCoA) of Bray-Curtis dissimilarity between samples based on prokaryotic genus level abundance tables inferred from the 75 metagenomes using mTAGs [30], with samples colored according to treatment with wastewater effluent (WW) either with or without ultrafiltration (UF). Samples with % WW or % WW UF indicated are biofilm samples, while remaining samples are water samples (95% confidence ellipses for different treatments). Figure S2. A) Mean relative abundance of eukaryotic phyla in each experimental treatment. Taxa which were in the upper decile of abundance in at least one experimental treatment were retained, while remaining phyla were binned as other. B) PCoA of Bray-Curtis dissimilarity between samples based on eukaryotic genus level abundance table inferred from the 75 metagenomes using mTags [30], with samples colored according to treatment. Samples with % WW or % WW UF indicated are biofilm samples, while remaining samples are water samples (95% confidence ellipses for different treatments). Figure S3. Top discriminative prokaryotic taxa as determined by LEfSe analysis out of the top 200 taxa shown. Node size corresponds to overall taxon abundance while shading indicates significant enrichment under a given treatment condition as determined by LEfSe analysis using workflows described by Segata et al. [31] Figure S4. PCoA of functional enzyme counts across treatment conditions based on a Hellinger-transformation of counts of enzyme classes annotated with Enzyme Commission (EC) number fourth-level annotations (see Methods). The Bray-Curtis dissimilarity metric was used. Samples with % WW or % WW UF indicated are biofilm samples, while remaining samples are water samples (95% confidence ellipses for different treatments). Figure S5. PCoA of DeepARG [34] subtype counts across treatment conditions. ARGs were normalized by 16S read counts. The Bray-Curtis dissimilarity metric was used (95% confidence ellipses for different treatments). Figure S6. PCoA-based Procrustes analysis showing a structural alignment between resistome composition and functional enzyme composition (EC number, fourth-level) and antibiotic resistance profiles. Points indicate position on PCoA axes 1 and 2 of resistome composition, while arrow tips indicate position on PCoA axes 1 and 2 of functional genes corresponding to fourthlevel EC counts were subjected to a Hellinger transformation followed by a Procrustes rotation. The correlation coefficient in a symmetric Procrustes rotation was 0.8901, the sum of squares (m12 squared) was 0.2078, and p-value was less than 0.001. Figure S7. Representative examples of sul1 genes detected on biofilm metagenomic contigs colocalized with mobile genetic elements such as recombinases, integron-integrase, and transposases. Figure S8. Chemical similarity between example A) sulfonamide-containing antibiotics and other B) sulfonamide-containing micropollutants detected in WW. Figure S9. Metagenomic hits to the Phycisphaeraceae metallo-beta-lactamase fold protein MBX3359331.1 shows higher counts in WW80 rather than WW30 and lower counts in the UF treated samples. No hits were obtained in the WW00 (stream water) samples. Figure S10. SDS-PAGE gel analysis of purified MBX3358097.1. The protein markers of 250 kDa (Thermo Scientific™ PageRuler™, cropped) is in lane 1. Lanes 2-9 show the elution fractions after expressing the 34 kDa protein. Table S1. Metagenomes sequenced and published in this study (n=75) from two independent experiments described by Desiante et al. [11] including experiment number, treatment type, sampling day (D = day), and number of raw reads per sample. Table S2. Parameters used for read mapping with bowtie2. Table S3. Information about the synthesized Phycisphaeraceae beta-lactamase genes in this study. Table S4. Mean concentrations of micropollutants (MPs) in biofilms normalized by biomass (ng MP/mg AFDW biofilm ± standard error) measured across all WW% summarized from Desiante et al. [11]. The selected group of 27 MPs were retained from the original dataset that consisted of 50 MPs, the selection criteria was discarding the MPs with more than 3 concentrations lower than the quantification limit (LOQ). Table S5. ARGs with highly significant differences (Kruskal-Wallis test + Benjamini-Hochberg p-value correction) in biofilms across all conditions.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Attrah, M., Schärer, M.R., Esposito, M. et al. Disentangling abiotic and biotic effects of treated wastewater on stream biofilm resistomes enables the discovery of a new planctomycete beta-lactamase. Microbiome 12, 164 (2024). https://doi.org/10.1186/s40168-024-01879-w

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40168-024-01879-w

Keywords