The current study investigated the resistome and taxonomy of the microbial community in municipal wastewater and the effects of biological treatment on these factors. The Village of Dunnottar Wastewater Treatment Facility (50° 26′ 39.3″N, 97° 01′ 04.0″W), located near the west shore of Lake Winnipeg’s south basin (approximately 1 h north of Winnipeg, Manitoba, Canada), is a municipal sewage lagoon servicing 989 permanent residents and a summertime (cottage) population of approximately 3,000. The lagoon consists of a primary, secondary and tertiary treatment cell (Supplementary Figure S2). The lagoon works through continuous flow and raw sewage is delivered to the primary cell by septic truck. The Dunnottar Wastewater Treatment Facility was used to study microbial resistome through both in situ sampling and as the source influent for a mesocosm experiment described below.

Supplementary Table S1 provides an outline of all samples taken throughout the study along with their location, material composition, date of collection and analysis performed.

Water and plant tissue were analyzed for metal concentration using an Agilent 7850 inductively coupled plasma mass spectrometer (ICP-MS following sample digestion with HCl and HNO₃ according to U.S. EPA SW-846 Method 3010A and Thermo Scientific AN 43446²⁰). Along with measuring metal concentrations in the mesocosms, metal concentrations were measured at each metagenomic sampling point in the primary and secondary lagoon on the corresponding day.

The results derived from the metagenomic analyses were imported into R (RStudio Team, 2020) for statistical analysis and data visualization. The average values and the standard deviations (when n > 1) of the different parameters evaluated were reported. A Bayesian-multiplicative replacement was applied to reduce the sparsity of the microbiome data using the R package zCompositions (Palarea-Albaladejo and Antoni Martín-Fernández, 2015). Data visualization was performed using ggplot2 (Wickham, 2016). The R package Metacoder (Foster et al., 2017) was used to generate differential heat-trees of the changes in the resistome and mobilome profiles after the treatments. For this, the log₂ of the fold change (log2 FC) of the total abundance of the different genes (ARGs, MRGs, and MGEs) was estimated and then encoded in the color of the plot. The number of genes in each level of the trees was used to define the size of the nodes and edges.

Spearman’s rank correlation analyses were used to evaluate the relationship between heavy metal concentrations, resistomes and mobilomes. These results were then visualized with correlograms created using appropriate R packages. Treatment effects were assessed using non-parametric tests (e.g., Kruskal-Wallis, Wilcoxin Ranked Sum). Only differences where statistical tests had p-values <0.05 were stated as significant.

Network analysis of ARGs, MRGs, biocide resistance genes, and bacterial genera was conducted using ggraph in RStudio (RStudio Team, 2020). The network was based on Bray-Curtis similarity, with minimum requirements for edges set to P = 0.70. The minimum number of connections per node was set to 2, and nodes present in less than 50% of sample locations were excluded from the network. These parameters ensured significant edge connections.

To scale bacterial abundance to reads (RPKM), the Hellinger scaling method was applied using the decostand function in the “vegan” R package (Oksanen et al., 2022). Genera from phyla with no edges (P > 0.70) linked to resistance genes were removed from the network analysis.

Total concentrations of each metal measured in wastewater varied over the growing season and the effect of each treatment varied across metals as well (Figure 2). Overall, there was a trend of metal removal throughout the growth season followed by a re-establishment of initial metal concentrations, likely due to biomass die-off and decomposition. Emerged macrophytes (FTW) treatment led to significantly (p < 0.05) greater Cu, Hg, Pb, As, and Cd removal than the other two treatment types, achieving an average metal removal on the final day (across all metal elements measured) of 59.49% (±7.45%) of initial concentration, as opposed to algae treatment [−117.02% (±69.53%)] and Lemna treatment [−32.51% (±8.31%)], and were the least affected by end of season metal release. Total metal concentration increased at the end of the growing season to levels higher than the initial concentration for mesocosms treated by algae. It is likely that this is due to the death and decomposition of the algae biomass releasing metals into the water column that were not initially suspended (Deng et al., 2016).

Lemna treatment was effective for certain elements (Cr, Cu, Fe, Ni, Pb, and Zn), even at the end of the growth season. It out-performed FTW treatment in iron (Fe) and chromium (Cr) removal, achieving 93.22% (±2.63%) and 95.97% (±1.59%), respectively. Metals were also measured in macrophyte root, shoot, and Lemna tissue at the end of the growth season (see Supplementary Figure S1). These results corroborate the findings of the water-metal analysis as greater concentrations of both iron and chromium were found in Lemna tissue (Fe = 2006.28 ± 229.77, Cr = 19.67 ± 1.16 mg/kg) than in macrophyte roots (Fe = 1574.11 ± 329.29 mg/kg, Cr = 8.86 ± 1.04 mg/kg) or shoots (Fe = 383.68 ± 53.97 mg/kg, Cr = 5.12 ± 3.18 mg/kg). With the exception of Hg, macrophyte roots accumulated greater concentrations of metals than the shoots.

The resistome and mobilomes within the microbial communities of the different samples are shown in Figure 3. A total of 138 different ARGs were found in the primary lagoon. ARGs encoding resistance to beta-lactams (47 genes, 25.1 rpkm), multidrug (24 genes, 30.5 rpkm) and aminoglycosides (22 genes, 6.78 rpkm) were the most abundant and diverse in this sample (Figure 3A). In terms of MRGs, chromium (Cr, 23.96 rpkm), multi-metals (11.1 rpkm), and mercury (Hg, 9.1 rpkm) resistance were the dominant in the primary lagoon sample (Figure 3B). The highest diversity and total levels of MGEs were also observed in the primary lagoon sample (981 genes, 75.5 rpkm), confirming the potential for the transmission of ARGs and MRGs, and the emergence of ARB in wastewater systems.

Interestingly, both the diversity (100 genes) and the total levels of ARGs (14.7 rpkm) and MRGs were substantially lower in the secondary lagoon, suggesting that the biological activity and environmental conditions in the lagoon systems resulted in the reduction of these genes. The most abundant ARGs in this sample were multidrug (34 genes, 5.4 ± 1.2 rpkm), beta-lactams (25 genes, 1.1 ± 0.35 rpkm), and aminoglycosides (17 genes, 2.2 ± 1.2 rpkm). The total levels of MRGs in the secondary lagoon (8.1 ± 2.2 rpkm) were 86% lower than those in the primary lagoon (59.4 rpkm), although a total of 46 different genes were detected. In terms of MGEs, the mobilome was more diverse (1,299 genes) but showed less abundance (899 ± 153 rpkm) in the secondary lagoon (Figure 3C). Cumulatively, these results suggest the capability of the biological activity in this lagoon system to potentially reduce the risk of spreading antimicrobial resistance (AMR).

The three biological treatments evaluated in this study reduced the diversity and total levels of ARGs (p = 0.05) and MRGs (p = 0.02) present in the secondary lagoon, which served as the wastewater source for the experiments. Mesocosms with floating macrophyte mats exhibited the lowest abundance of MRGs (p = 0.04) and MGEs (p = 0.04), while mesocosms treated with Lemna had the lowest abundance of ARGs (p = 0.10).

In the FTW mesocosms, only 12 distinct ARGs were detected, with total ARG levels (5.0 ± 0.33 rpkm) reduced to 34% of those in the secondary lagoon (Figure 2A). ARGs conferring resistance to aminoglycosides, MLS, and fluoroquinolones were the most reduced groups. This treatment also decreased the diversity of MRGs (9 genes) and their total abundance (1.0 ± 0.6 rpkm, p = 0.01) by 20% and 88%, respectively (Figure 2B). Specifically, MRGs encoding resistance to Ag, As, Zn, Hg, and Fe showed lower levels, while resistance to Cu and multi-metal elements increased (Figure 3B). A total of 578 MGEs were detected in the FTW Mesocosm, representing a 55% reduction compared to the source. Although total MGE levels were also reduced (515 ± 96.2 rpkm), the decrease was not statistically significant (p = 0.15, Figure 2C). MGEs related to transfer mechanisms and phages exhibited the highest reductions (Figure 3C). In contrast, the abundance of MGEs involved in integration and excision functions increased, suggesting the presence and proliferation of microbial groups carrying these MGEs in this treatment.

In the Lemna mesocosms, only 12 ARGs were detected, reflecting an 88.5% reduction from the secondary lagoon. Furthermore, total ARG levels decreased by 87.5%, with beta-lactam, aminoglycoside, and MLS ARGs being the most reduced groups (Figure 4A). The diversity of MRGs (9 genes) decreased by 80%, and their total abundance (4.2 ± 1.0 rpkm) was reduced by 48%, although this reduction was not statistically significant (p = 0.11). MRGs encoding resistance to Hg, Ag, As, Cu, Fe, and Zn were lower in abundance (Figure 4B). Fewer MGEs (677 genes) were detected compared to the secondary lagoon, though total MGE levels (1,318 ± 492 rpkm) tended to be higher. This increase, while not statistically significant, appeared to result from the enrichment of MGEs related to integration and excision (Figure 4C).

The Algae mesocosms had 20 detectable ARGs, reducing total levels by 70.5%. ARGs encoding resistance to MLS, beta-lactams, and aminoglycosides were found in lower concentrations (Figure 4A). For MRGs, 18 genes were detected, with total levels trending lower (6.18 ± 4.3 rpkm), although this was not statistically significant (p = 0.26). MGE diversity (756 genes) was reduced by 42%, but total MGE levels (1,043 ± 257 rpkm) showed a 14% increase, though this change was also not statistically significant (p = 0.545). Phage-associated MGEs and those involved in replication, recombination, and repair were generally reduced, while integration/excision MGEs showed enrichment (Figure 4C).

Overall, the results showed that the total levels of ARGs and MRGs in all three treatments were lower than those observed in the secondary lagoon control, suggesting that the treatments effectively reduced the risk of AMR transmission by decreasing the abundance of these resistance genes. While past research has shown the statistical comparability of metagenomics and qPCR (Flores-Orozco et al., 2022), further research should be conducted to validate and expand on these findings using qPCR and antibiotic susceptibility testing.

The microbial communities at the phylum levels are shown in Figure 5, whereas the abundance of the top 40 most abundant genera is displayed in the Supplementary Figure S4. The most abundant phylum in the primary lagoon was Pseudomonadota (91.5%), which was mainly composed of the genera Aeromonas (68.5%), Pseudomonas (6.8%), Shewanella (4.6%), and Escherichia (1.9%). The microbiomes of the secondary lagoon samples were also dominated by the Pseudomonadota (61.2%), although Bacteroidota (13.7% ± 10.4%), Cyanobacteriota (8.9%), and Actinomycetota (8.5%) also represented a significant part of the community. At the genus level, Flavobacterium (11.3%) had the highest abundance, followed by Acidovorax (6.5%), Sphingopyxis (6.5% ± 6.0%), and Enterobacter (5.9%).

The microbial communities showed distinct shifts during the different treatments, especially an increase in Cyanobacteriota at the expense of Pseudomonadota. For example, the microbiome in the Algae mesocosms was characterized by high levels of the phylum Cyanobacteriota (58.2%), mainly composed of Cyanobium (23.9%), Nostoc (17.8%), and Trichormus (8.5%). Pseudomonadota was the second most abundant phylum, accounting for about 30.3% of the microbial communities, and Aeromonas (7.1%) and Delftia (3.7%) were the dominant genera in this phylum. In the Lemna mesocosms, the phyla Cyanobacteriota (67.9% ± 11.3%) and Pseudomonadota (23.7%) also dominated the microbial communities. However, in this treatment the genus Nostoc (35.3% ± 12.3%) was the most abundant, followed by Trichormus (16.9%) and Cyanobium (9.2%). On the other hand, Pseudomonadota (55.4%), Bacillota (16.4%), and Actinomycetota (15.3%) were the most abundant phyla in the FTW treatment samples. At the genus level, Delftia (13.5%), Cutibacterium (10.4%), and Staphyloccocus (7.9%) had the highest abundance in this treatment.

Across all the metagenomic samples in this study a Spearman’s correlation between the abundance of ARGs and MRGs was r = 0.70 (including rhizome biofilm samples). Co-occurrence and co-selection of ARGs and MRGs have been found in many past studies (e.g., Li et al., 2016; Amarasekara et al., 2023; Li et al., 2022). There are three accepted models that explain ARG/MRG co-selection, which are co-resistance, cross-resistance, and co-regulation (Baker-Austin et al., 2006; Gillieatt and Coleman, 2024). Cross resistance refers to ARGs and MRGs that confer resistance to both antimicrobials and metals, co-resistance refers to ARGs and MRGs that share close genetic linkages, and co-regulation refers to resistance genes that share the same regulatory systems (Gillieatt and Coleman, 2024). Co-resistance may lead to resistance genes that are maintained in clusters as adjacent genetic elements (Li et al., 2016). The resulting co-residence or clustering of ARGs and MRGs on MGEs (such as integrons, transposons, and plasmids) allow for HGT and further explain the widespread co-occurrence observed across the microbial population (Amarasekara et al., 2023). In the current study, a significant positive correlation (r = 0.62, p < 0.05) was observed between MRGs and MGEs, suggesting that MRGs were to some extent present and carried on MGEs. There was strong correlation between specific ARGs and MRGs-most notably: Fe resistance genes and antibiotic target replacement (r = 0.90, p < 0.05) and antibiotic target alteration (r = 0.86, p < 0.05) genes; Hg resistance genes and antibiotic efflux pumps genes (r = 0.89, p < 0.05); and Organic-mercury (Hg) resistance genes and antibiotic inactivation (r = 0.89, p < 0.05), antibiotic target replacement, (r = 0.90, p < 0.05) and antibiotic target alteration genes (r = 0.91, p < 0.05).

Past research has shown specific correlation of the presence of certain metal class MRGs and antimicrobial class ARGs. These include the co-occurrence of zinc (Zn) resistance with beta-lactam resistance, copper (Cu) resistance with beta-lactam resistance, and arsenic (As) resistance with beta-lactam resistance (Li et al., 2016). Figure 6 presents correlations between specific metal compound resistance and ARG abundance categorized by molecular mechanism. In this study, Cu resistance was significantly correlated with antibiotic inactivation (Figure 6), the main mechanisms by which beta-lactams act, providing more evidence for the co-occurrence of Cu resistance and beta-lactam resistance.

In the present study, there was a strong correlation between Cu resistance and aminoglycoside resistance (r = 0.76), organo-mercury resistance and aminoglycoside resistance (r = 0.95), mercury resistance and aminoglycoside resistance (r = 0.81), and iron resistance and aminoglycoside resistance (r = 0.79). Beta-lactam resistance was correlated to mercury (r = 0.91) and organo-mercury resistance (r = 0.80). Multi drug resistance was found to be correlated to iron (r = 0.89), mercury (r = 0.85), and organo-mercury resistance (r = 0.79). Additionally, across all samples, iron resistance had a strong correlation to rifamycin resistance (r = 0.96). A strong correlation existed between iron resistance and antibiotic target protection (r = 0.87), and antibiotic target replacement (r = 0.86) (Figure 6). Specific gene correlations exist between iron resistant genes acn and brfA, and antibiotic inactivation genes conferring resistance to aminoglycosides (aaDA12, aaDA17, aaDA2) and macrolides (mphE, mphG). Mercury resistance genes merD, merT, and merR also had strong correlation to antibiotic inactivation genes aaDA12, aaDA17, aaDA2, aaDA3, aaDs, mphE, mphF, mphG, EreD.

Corelograms were used to visualize Spearman’s rank correlation of heavy metal concnetration in wastewater and resistance gene abundance. Figure 7 shows the correlation of the 6 different heavy metals measured throughout the study to the abundance of total ARGs, MRGs, and MGEs.

Strikingly, a strong and statistically significant correlation (r = 0.95, p < 0.05) between total ARG abundance and Fe concentration was observed (Figure 6A). This high correlation suggests Fe concentration could influence the presence of ARGs in the microbial communities in these types of environments. Unlike other measured metals, such as Cu, As, Cr, and Hg, Fe does not have a high level of toxicity and is generally innocuous at environmental concentrations. Iron has been found to cause both increased antibiotic resistance and susceptibility in lab experiments, depending on the bacteria and antibiotic (Ezraty and Frédéric, 2016).

Iron plays a role in certain processes related to antimicrobial mechanisms such as efflux pumps, respiration and Fenton process-dependent biocidal activity. Choi et al. (2022) found that Fe supplementation promoted antibiotic resistance in Streptomyces coelicolor and Mycobacterium smegmatis by upregulating genes involved in metabolism and respiration, even under anaerobic conditions. This resistance involves a potassium cyanide-insensitive respiration type, suggesting a strategy that helps bacteria avoid oxidative damage from bactericidal antibiotics. To the best of our knowledge, the present study is the first to look at the effect of Fe concentration on ARGs in an environmental setup. The strength of the correlation suggests that another effect, beyond the previously discussed mechanism, may be responsible for the observed increase in ARG abundance in the presence of higher Fe concentrations. Rather than Fe directly influencing the expression of ARGs, it may promote the growth of bacteria genera and species more likely to harbor and procure resistance, such as certain facultative Pseudomonadota bacteria (further discussed below). Further studies are needed to understand the molecular and evolutionary mechanisms by which Fe influences ARGs.

Of all the phyla identified in the study, Pseudomonadota was the most strongly correlated with ARGs (r = 0.77, p < 0.05), followed by campylobacteria (r = 0.66, p < 0.05). To further investigate the relationship between bacteria and ARGs, a network analysis was performed (Figure 8).

There was a close relationship between clusters of ARGs, MRGs and biocides and labelled bacteria genus (Figure 8). Tightly clustered groupings of resistance genes represent genes that were similar (P > 0.70) in presence and abundance across the different samples, indicating that they were related, possibly through chromosomal proximity on MGEs or through presence in certain bacterial species. Bacteria genera in the network analysis that are connected to resistance clusters are also similar in presence and abundance across the different samples, suggesting that they could harbour the resistance genes with which they are clustered. It is noteworthy that the majority of labelled bacteria (those with the strongest similarity to resistome abundance) were genera belonging to Pseudomonadota, followed exclusively by Campylobacteria. These results concur with the relationship identified through Spearman’s rank correlation, despite using different analysis methods (Bray-Curtis’s similarity).

Of the eight genera identified through the network analysis to be clustered with resistance genes, four (Aeromonas, Enterobacter, Laribacter, and Rhodoferax) were facultative Pseudomonadota. Aeromonas, Enterobacter, and Laribacter contain pathogenic bacterial species often associated with resistance, such as A. hydrophila and A. salmonicida ( Pintor-Cora et al., 2023), and L. hongkongensis (Lau et al., 2011). Aliarcobacter, belonging to the phylum Campylobacteria, is known to harbour metal and antibiotic resistance genes (Müller et al., 2020). Sulfurosprillum spp., also belonging to the Campylobacteria phylum, are known to reduce sulfur and metal species such as arsenate, and for living in heavy metal-rich environments (Goris and Diekert, 2016).

Neighbour (direct node connections) and next of neighbour (secondary node connections) for each of the eight bacterial connections were identified. These genes represent the resistance genes most similar in abundance and occurrence to each genus. Aeromonas had five neighbours (CEPH-A3, OXA-12, cph-A6, OXA-724, cph-A4) and one next of neighbour gene (cph-A7) belonging to the beta-lactam resistance class. Aeromonas is known to commonly harbour resistance to beta-lactam antimicrobials, including cefepime and other extended-spectrum cephalosporins (Aravena-Román et al., 2012). Six out of 10 neighbour and next of neighbour resistance genes clustered with Aeromonas from the network analysis were beta-lactam ARGs and four of these six confer resistance to cephalosporins.

The congruence between the findings from the present study and known Aeromonas resistance in literature provides an indication of validation for the network analysis. Both Laribacter and Enterobacter bacteria were clustered with a high percentage (6 out of 6 and 5 out of 8, respectively) efflux pump ARGs and MRGs. L. hongkongensis has been found to harbour efflux pump resistance genes, including Acr eflux system genes (Lau et al., 2011), which were present in both Laribacter and Enterobacter resistance clusters. Aliarcobacter was directly connected to only one resistance gene, fetB/ybbM, an MRG conferring resistance to Fe, while its next of neighbours included other MRGs conferring resistance to As (pstB) and to multi-metals (czcA). Aliarcobactor’s next of neighbours also included beta-lactam resistance genes.

Given the strong correlation between Fe concentration and ARGs across the wastewater samples, the relationship between the eight genera connected to resistance clusters identified from the network analysis and Fe was further explored (Table 1).

The three genera most strongly correlated with ARGs (Enterobacter, Laribacter, and Aliarcobacter) were also the three genera most strongly correlated with Fe resistance and Fe concentration (Table 1). The majority of the genera identified by the network analysis were facultative or aerotolerant bacteria that can utilize Fe in their respiration systems. The strong correlation of Fe and ARGs across the wastewater samples (r = 0.95, p < 0.05) and the association of facultative, Fe reducing bacteria with ARGs in this study suggests that Fe promotes the proliferation of Pseudomonadota and Campylobacteria harboring ARGs.

Metal ions are used in anaerobic respiration employed by facultative bacteria. Fe³⁺ and Mn⁴⁺ can be used as alternate electron acceptors (Fenchel et al., 2012). Respiratory Fe reduction has been found to play significant roles in anaerobic carbon transformations in freshwater sediments (Fenchel et al., 2012). Iron reduction [Fe(III) to Fe(II)] pathways have been identified in bacterial species belonging to Aeromonas, Enterobacter, Laribacter and Rhodoferax genera (Table 1).

Sulfur metabolization in sulfur-reducing bacteria is also associated with anaerobic respiration. Sulfur-metabolizing genes (SMGs) were identified utilizing the Sulfur Metabolizing Data Base (SMDB) and were found to be correlated with ARGs (r = 0.63, p < 0.05), and Pseudomonadota (r = 0.76, p < 0.05) and Campylobacteria (r = 0.87, p < 0.05). Additionally, correlations between MRGs (r = 0.51, p < 0.05), and Fe resistance (r = 0.75, p < 0.05), Fe concentration (r = 0.78, p < 0.05), and total metal concentration (of measured metals) (r = 0.83, p < 0.05) exist with sulfur-metabolizing genes. Reduction of ferric iron (Fe³⁺) and sulfate (SO₄ ²⁻) are the most common electron-accepting processes for anaerobic respiration in anoxic freshwater environments (Berg et al., 2019). Iron also plays a role in Fe(III)-driven sulfide oxidation (Berg et al., 2019). It is therefore probable that Fe provides advantages to each genus identified in Table 1, especially in anoxic or low oxygen environments.

A relationship between facultative and aerotolerant bacteria, and ARGs is strongly suggested by the results of this study. Facultative bacteria have been shown to play a role in ARG dissemination in municipal and livestock wastewater (Hembach et al., 2019; Flores-Orozco et al., 2024). and are overrepresented on the WHO antibiotic-resistant Priority Pathogens list, accounting for 8 of 12 listed pathogens (André et al., 2021).

The use of phytoremediation systems such as emerged floating macrophytes and Lemna corrals may therefore reduce ARGs in wastewater through the removal of metals, including Fe, limiting facultative bacterial growth that are dependent on electron acceptors for anaerobic respiration.

The molecular and evolutionary mechanisms that may explain why facultative bacteria seem more likely to harbor higher-level ARGs than other bacteria remain unknown. Further genomic investigation is recommended to study this relationship.

The results of the present study suggest that phytoextraction of metals in wastewater could reduce the abundance of MRGs and ARGs in wastewater effluents. This is of importance to wastewater treatment lagoon and plant operators as phytoremediation systems such as FTWs or Lemna treatment corrals, provide an affordable and practical augmentative treatment to improve lagoon performance (McCorquodale-Bauer et al., 2023). Additionally, the results are relevant beyond the scope of phytoremediation systems, as metal removal through more traditional constructed wetlands or other treatment systems may result in a decrease in MRG and ARG abundance. More research is needed to understand how FTWs would affect metals, MRGs and ARGs at a full-scale level.