Twenty-four solid end products (biosolids) from two municipal wastewater treatment plants (WWTPs) located at city A (WWTP-A) and city B (WWTP-B), within 100 km distance, were collected at 10-year intervals in 2009 and 2019, respectively. The WWPT-A at a more prominent city treats waste-activated sludge by anaerobic digestion and dewatering process with the end product of wet pellets. WWTP-B treats waste-activated sludge in a dewatering/pelletization process involving the belt-filter press and final pelletization process by thermal drier (250–450°C at the entry, and 80–130°C at the exit). Samples were taken in triplicate twice a year with an interval of 5-6 moths per plant (n = 24 in total). 500 g per sample were taken, and 10 g were used for DNA extraction. The biosolids were either wet or dry pellets in irregular shape with approximately diameters between 0.1—0.5 center meter and the moisture contents approximately ranging from 59—73% for WWTP A and 3—11% for WWTP B. The samples were transported on ice and were stored at −20°C until further use.

In addition, for the samples used in the current study, live coliforms, generic E. coli and Salmonella were detected for the biosolids collected in 2009 as described in a separate publication for which the WWTPs were named as WWTP A and C (Yergeau et al., 2016; Duceppe et al., 2019). Similar results for coliforms and generic E. coli have been obtained for the samples used in 2019 with the WWTPs referred as WWTP A and B respectively) (data not shown).

Biosolids were first ground manually with mortar and pestle to a fine powder. DNA was extracted using the PowerMax^® Soil DNA Isolation Kit following the manufacturer’s instructions (QIAGEN, Germany) and quantified using Qubit™ four Fluorimeter (Life Technologies, United States). The DNA extracts were stored at −20 °C and retrieved for ARGs detection using qPCR array and for bacterial host analysis using 16S ribosomal RNA (rRNA) gene sequencing.

Clinically relevant 84 ARGs (Supplementary Table S1) were detected with a commercial ARG qPCR array kit from QIAGEN (Cat# 330261 BAID-1901ZR) using ABI 7500 Fast real-time PCR system (Applied Biosystems, United States), according to the manufacturer’s instruction. Briefly, the 25 μl reaction volume consisted of 2x microbial qPCR master mix (QIAGEN, Germany) and 52 ng DNA for each well. The manufacturer’s instruction recommended minimum of 500 ng/plate or 5.2 ng/well. In order to increase the detection sensitivity, various concentrations of DNA extracts were tested in our preliminary study which showed the increased numbers of AMR genes detected with the increase of the DNA concentrations. 5,000 ng/plates or 5.2 ng/well (or reaction) was used based on the numbers of AMR gene detected and concentrations of the DNA extract obtained from the same DNA extraction. Plates were incubated at 95 °C for 10 min followed by 40 cycles of 95°C for 15 s and 60°C for 2 min. Threshold cycle (Ct) values were validated and normalized using an Excel data analysis template provided by QIAGEN. A Ct value below 34 was considered positive based on the manufacturer’s instruction, indicating an ARG was present in the sample. The heatmap was plotted using the R package pheatmap v1.0.12 (https://github.com/raivokolde/pheatmap). In addition, the binary data of ARG presence and absence were imported into R v4.1.2, and principal coordinates analysis (PCOA) was performed using Vegan v2.5-7 (https://github.com/vegandevs/vegan).

16S rRNA gene amplicon sequencing was performed as described previously (Laskey et al., 2020). Briefly, the V3–V4 region of the 16S rRNA gene was amplified by PCR. Libraries were prepared and sequenced using a MiSeq system (Illumina, United States). Raw read data was demultiplexed and then analyzed by Qiime2 (Bolyen et al., 2019) using a modified version of the Qiime2 pipeline (AmpliconPipeline, https://github.com/forestdussault/AmpliconPipeline). The biom, tree, and metadata files generated from this pipeline were combined into a phyloseq object using the R package microbiomeMarker v1.0.1 (https://github.com/yiluheihei/microbiome Marker). R packages phyloseq v1.38.0 (McMurdie and Holmes, 2013), ggplots2 v3.3.5 (https://github.com/tidyverse/ggplot2), and vegan v2.5-7 (https://github.com/vegandevs/vegan) packages were used for further analysis and visualization.

Biomarker analysis was performed using the linear discriminant analysis (LDA) effect size (LEfSe) method (Segata et al., 2011) using the R package microbiomeMarker v1.0.1 (https://github.com/yiluheihei/microbiomeMarker). A threshold score of LDA score = 4 and Kruskal-Wallis test α value = 0.05 were set as thresholds to detect differences.

A correlation matrix and significance levels were generated using the R Hmisc package (https://github.com/harrelfe/Hmisc) to calculate Pearson correlation with significance levels (p-value) between 259 bacterial genera and 71 ARGs presented in 24 samples. A correlation was considered statistically robust if the Pearson’s correlation coefficient (r) was ≥0.8 and the p-value was ≤0.01 (Junker and Schreiber, 2008). The significant and robust correlation visualization was performed using GraphPad Prism 7.0 software (United States).

Differences in the alpha diversity (Chao1 and inverse Simpson index) between the groups were analyzed using the Wilcoxon tests in the R package ggpubr v 0.4.0 (https://github.com/kassambara/ggpubr). The difference in Beta-diversity (PCOA plot) between different years, locations, or years and location combined were analyzed using PERMANOVA followed by Betadisper tests (Vegan v 2.5-7). The data should be homoscedastic (p-value was not significant in the Betadisper test) to apply a PERMANOVA. A p-value < 0.05 was considered statistically significant.

A total of 71 out of 84 targeted ARGs, representing nine different antibiotic classifications, were detected with various abundances (Figure 1A and Supplementary Table S1), which was consistent with a previous study that reported the detection of 70 ARGs (Agga et al., 2015). Moreover, only 12 ARGs (aadA1, bla _GES , bla _VEB , bla _OXA-2 , ereB, qnrB-5, qnrB-8, qnrS, ermB, mefA, tetA, and tetB) were detected in all individual samples (Figure 1A and Supplementary Table S1). Five macrolide resistance genes (ermA, ermB, ermC, mefA, and msrA) were detected in most samples. Particularly, ermB has the highest target copy numbers (the lowest average Ct value), a typical characteristic gene detected in urban sewages in North America and Europe (Hendriksen et al., 2019). Numerous studies on the characterization of ARG or ARB in wastewater have been conducted over the past decades using various methods to facilitate the detection of ARGs in biosolids or sludge (Viau and Peccia, 2009a; Novo et al., 2013; Agga et al., 2015; Yergeau et al., 2016; Li et al., 2020; Law et al., 2021; Zielinski et al., 2021; Jankowski et al., 2022). In contrast to these studies, we investigated the changes of the antibiotic resistomes at two different facilities over the 10-year interval and detected a total of 58 and 53 ARGs in samples collected from WWTP-A in 2009 and 2019, respectively. Whereas 66 and 53 ARGs were detected in samples from WWTP-B in 2009 and 2019, respectively. A slight reduction of the total number of ARGs from 69 (in 2009) to 63 (in 2019) was observed in the two facilities combined. Only two genes, vanC and mecA believed to be exclusively found in Gram-positive cocci were detected in the samples collected in 2009 at both WWTPs (Figure 1A and Supplementary Table S1), indicating that Gram-positive bacteria containing these ARGs were probably present in both WWTPs at the early sampling date but not 10 years later or the wastewater treatment efficiencies improved over time and subsequently eliminated the bacteria. Moreover, several ARGs, such as bla _BES-1, bla _IMP-5, bla _CMY-10, and bla _FOX, were consistently detected in the samples from WWTP-B, while the bla _ACC-1 group was only present in the samples from WWTP-A (Figure 1A and Supplementary Table S1).

β-diversity analysis was used in this study to estimate the overall similarity in the antibiotic resistomes between the samples collected from two different sites on two occasions 10 years apart. There was no clear distinction between the antibiotic resistomes clustered for the grouped samples based on the year and site combined (Figure 1B), years (Figure 1C), or sites alone (Figure 1D), indicating that the antibiotic resistomes are similar to each compared group. Additionally, antibiotic resistome in samples collected from WWTP-A in both 2019 and 2009 have a longer distance to the centroid, indicating antibiotic resistome in samples collected from WWTP-A are more variable than WWTP-B (Supplementary Figure 1S). Compared to antibiotic resistomes of WWTPs in the world (Hendriksen et al., 2019), both WWTP-A and -B show very similar ARG profiling to that of WWTPs in Canada and United States, which further confirm diversity of antibiotic resistomes is consistent in the close geographical regions.

In recent years, effective actions, including a drastic reduction in antibiotics consumption, reduction in the uncontrolled use in humans and animals, and the ban of the sub-clinical dosage growth promoter) usage of antibiotics for animal and crop productions to promote growth have been in place to reduce the spread of ARGs globally (WHO, 2015; Merlin, 2020) and in Canada (McCubbin et al., 2021). It seems that the measures to reduce the use of AMRs in animals were effective in reducing ARGs in these animals (Dorado-Garcia et al., 2016; McDonnell et al., 2017; Bennani et al., 2020). However, it is not clear whether this translated into ARG reduction in humans. In this study, a slight reduction of the total number of ARGs was observed, but overall, there were no apparent shifts of ARGs in biosolids over 10 years in two WWTPs in Canada. ARGs or related genetic elements have existed in the environment for a very long time, even before the antibiotic era. These genes can be transferred from the environment to human-related activities or the opposite direction, but it is not easy to get rid of them from our world. Eventually, their concentrations on earth will increase if more clinically important ARGs or ARBs show up due to the misuse in non-clinical settings. Our results showed no significant increase in the total number of ARGs, and there were no apparent shifts of ARGs in biosolids over 10 years in two WWTPs in Canada, suggesting a reduction or slowing down or stabilization of the speed of ARGs dissemination from human sources into the environment (Viau and Peccia, 2009b). More investigations are needed to evaluate the potential sources that contribute to the circulation and accumulation of ARGs in biosolids, and the potentially transfer to food crops.

V3-V4 hypervariable regions of the 16S rRNA gene amplicons sequencing were used to characterize bacteria communities in WWTPs at different sites and years. Overall, compared to the WWTP-A, higher microbial richness (Chao1 index) was observed in the samples collected from WWTP-B in 2009 and 2019, indicating more diverse bacterial species present in the biosolids samples collected from WWTP-B (Figure 2A), compared to that of WWTP-A.

Neither samples collected from the two WWTPs showed a significant difference for each site between 2009 and 2019 based on the Chao1 index (Figure 2A). Similar to the Chao1 index findings, microbial diversity (Inverse Simpson index) was similar in the samples collected from WWTP-A within 10 years. However, for the WWTP-B, samples collected in 2019 had higher microbial diversity than that of 2009, indicating that the microbiota community is more heterogeneous in the 2019 samples than those in 2009. The samples collected in 2019 from WWTP-B also showed different microbial diversity compared to the samples collected from WWTP-A in 2019 (Figure 2A). The two WWTPs used different treatment procedures resulting in different moisture contents as explained above, and different bacteria counts, both higher for samples from WWTP-A (Yergeau et al., 2016; Materials and methods—Study site and sample collect). These results indicate that the different WWTP procedures could be one of the factors contributing to the differences of microbiota and the subsequent differences of ARBs from the samples of these two facilities.

Using Bray-Curtis distances, the PCoA of bacterial microbiota divided the samples into four separate clusters based on years and sites (Figure 2B). The bacterial microbiota compositions in the samples collected from WWTP-A and -B differed in 2009 (Figure 2B). After 10 years, there were significant parallel shifts in the microbiota composition in both facilities, especially samples collected from WWTP-A (Figure 2B). Additionally, a permutational multivariate ANOVA (PERMANOVA) exhibits a significant difference in the microbiota composition in biosolid samples collected from different sites and years, respectively (Table 1).

A total of 34 phyla were observed across all biosolid samples, with most sequences classified to one of three phyla: Proteobacteria, Firmicutes, and Bacteroidetes (Figure 2C), which were also the most frequently identified in a variety of samples collected from different WWTPs (Nascimento et al., 2018; de Celis et al., 2020; Niestepski et al., 2020). Within these dominant phyla, taxa typically associated with human gut microbiota (e.g., Aeromonas, Acidovorax, Dechloromonas, Acinetobacter, Arcobacter, Flavobacterium, Bacteroides, Bacillus, Christensensellaceae, Clostridium, Eubacterium, Erysipelotrichaceae, and Ruminococcaceae) existed in our samples, which is in agreement with previous publications, showing the detection of human gut microbiota in the urban sewage (Newton et al., 2015; Su et al., 2017). The prevalence of human gut microbiota in biosolid samples is not surprising as most human wastes enter WWTPs directly. Additionally, bacteria associated with the environment, such as water and soil (e.g., Zoogloea, Devosia, Ferruginibacter), were also found in our biosolids samples. These results suggest that microbial communities in municipal biosolids are predominantly composed of human and environmental-associated bacteria.

Characterizing microbial communities between different years and sites is key to identifying unique microbial contributions. Predominant taxa differentially associated with each group were determined by LEfSe. The results showed that the bacterial community in samples from WWTP-A mainly belonged to the phyla Cloacimonetes, Verrucomicrobia, and Acidobacteria for 2009 and the phyla WS6, Planctomycetes, and Atribacteria for 2019 (Figure 2D). In contrast, biosolids samples collected from WWTP-B were better characterized by the presence of Bacteroidetes, Actinobacteria, and Saccharibacteria for 2009 and Proteobacteria and Fibrobacteres for 2019 (Figure 2D). Further LDA effect size-based cladogram showed Proteobacteria and Planctomycetes phyla, especially Planctomycetacia class, enriched in biosolid sample from 2019, while bacterial species enriched in biosolid samples collected from 2009 were taxa in the Spirochaetae phylum and Actinobacteria class (Figure 2E). The bacterial communities inhabiting biosolid collected from WWTPs can vary between different treatment facilities over time, due to the differences in treatment procedures between various facilities and the sources of WWTP influents (Yergeau et al., 2016; de Celis et al., 2020; Wolters et al., 2022). However, interestingly, given the fact that there is no significant differences in antibiotic resistome between these two facilities (Figure 1D), observation of significantly different microbiota compositions between these two facilities indicated potential ARG transfer between different bacterial species or existence of large amount of cell-free ARGs in biosolids.

It is well known that a particular microbial community in biosolids-derived fertilizers is thought to play an essential role in enabling and facilitating particular functions to enrich nutrition in the soil (Hu et al., 2019; Pan et al., 2021). In our data, a significant increase in Proteobacteria and Planctomycetes phyla populations was found in 2019 at both facilities. The members of Planctomycetes and Proteobacteria phyla are capable of several metabolisms and biosynthesis, such as nitrogen-fixation and ammonia oxidation (Delmont et al., 2018; Rivas-Marin et al., 2019; Stultiens et al., 2020), and provide an advantage to biosolids as fertilizers to enrich soil microbial community. Future studies predicting the functions of these populations or communities can potentially improve the efficiency of WWTP processes.

Many studies have already demonstrated the presence of a broad range of human bacterial pathogens in wastewater and biosolids samples such as Escherichia coli, Mycobacterium tuberculosis, Pseudomonas aeruginosa, Salmonella enterica, Shigella flexneri, Staphylococcus aureus, Yersinia enterocolitica, etc. (Law et al., 2021; Nguyen et al., 2021; Young et al., 2021; Bai et al., 2022; Jankowski et al., 2022). In our study, due to the limitations of the partial regions of 16S rRNA amplicon sequencing, it was impossible to achieve microbiota’s taxonomic resolution at the species level. However, the microbial community in these biosolids samples consists of several genera, such as Bacteroides, Clostridium, Escherichia, Mycobacterium, which contain some human pathogens. Further studies, such as the full-length 16S rRNA gene amplicon sequences and analysis using Nanopore sequencing platform or bacterial culture, needed to be done to assess the potential environmental risks.

To identify potential ARG hosts, correlation analysis of 71 ARGs abundance and 259 bacterial genera abundance was evaluated based on strong (r > 0.8) and significant (p < 0.01) correlations (Figure 3). Co-occurrence patterns of ARGs and microbial genera in the biosolids could indicate possible host information of ARGs if the ARGs and microbial taxa display strong and significant positive correlations (Nguyen et al., 2021). Exploration of connections between ARGs and bacterial genera showed that 17 out of 20 genera were from the top three phyla: Firmicutes, Bacteroidetes, and Proteobacteria, and most of ARGs belonged to the beta-lactam resistance group (Figure 3). Interestingly, seven genera, including Ruminococcaceae NK4A214 group, Ruminococcaceae UCG 013, Fluviicola, Bacteroides, Fibrobacter, Aeromonas, and Dechloromonas, associated with human gut microbiota, contain ARGs related to aminoglycoside and beta-lactam resistance. In addition, members of the Ruminococcaceae NK4A214 genus are possible hosts for multiple genes responsible for the resistance of β-lactamase-resistant enzymes (Figure 3). Previous studies have revealed the presence of ARB, including specific pathogens in WWTPs (Nguyen et al., 2021). Similar to our results, Proteobacteria have been shown to contain the highest number of ARB, followed by Firmicutes and Bacteroidetes (Liu et al., 2019). In the Firmicutes phylum, Ruminococcaceae NK4A214 and Ruminococcaceae UCG 013 are members of the Ruminococcaceae family, a significant butyrate producer and essential bacteria involved in intestinal health (Vital et al., 2017). However, the Ruminococcaceae family in WWTP might be the host of qnrS, bla _TEM, bla _SHV (Narciso-da-Rocha et al., 2018). Moreover, members of Aeromonas, Bacteroides, and Mycobacterium genera containing human pathogens have been reported to carry several ARGs (Barancheshme and Munir, 2017; Hultman et al., 2018; Stalder et al., 2019). Our results also revealed typical environmental bacteria species belonging to Ferruginibacter, Rhodobacter, and Arenimonas phyla might carry genes that resist multiple types of beta-lactam antibiotics (Figure 3). Further Hi-C (a high throughput sequencing technique which captures in DNA-DNA proximity) and shotgun metagenomics approaches (Stalder et al., 2019) will be used to confirm the correlated hosts in the near future.

Several studies have shown the presence of live bacteria and ARBs in biosolids, indicating that treatments in WWTPs cannot completely inactivate all the bacteria (Alam et al., 2007; Viau and Peccia, 2009b; Yergeau et al., 2016). For the samples used in the current study, live coliforms, generic E. coli and Salmonella were detected for the biosolid samples collected in 2009 as described in separate publications for which the WWTPs were named as WWTP A and C (Yergeau et al., 2016; Duceppe et al., 2019). Similar results for coliform and generic E. coli have been obtained for the samples used in 2019 with the WWTPs referred as WWTP-A and -B respectively) (data not shown). Live ARBs have an opportunity to interact with soil microbiota and spread ARGs into the soil bacterial community via the horizontal gene transfer (Woegerbauer et al., 2020). With the annual application of biosolid fertilizers, the ARGs from biosolids may persist (Mays et al., 2021) or accumulate in the soil (Xie et al., 2016; Zhuang et al., 2021), which may continue to disseminate ARGs to humans potentially through various potential ways including food chains (EFSA BIOHAZ Panel et al., 2021). Because AMR is an increasingly severe threat to public health globally, concurrent actions from all the countries in the world using a One Health approach are necessary (Prestinaci et al., 2015; Murray et al., 2022). Currently, most developed countries have established surveillance programs for ARGs or ARBs. Still, these programs only investigate limited and common ARBs, such as pathogenic E.coli and Salmonella (PHAC, 2020), and do not cover the majority of AMR hosts. Our study has detected many environmental and gut commensal bacteria potential carrying ARGs, pressing the need to conduct comprehensive investigation or detection of ARBs, including pathogenic and non-pathogenic ARBs in the monitoring programs to assess the risks to the food chain. Compared to stable antibiotic resistome in these biosolid samples (Figure 1), significant changes of biosolid microbiota composition between the two WWTPs were observed in our study (Figure 2), which may explain why only a few possible host information of ARGs were identified using the correlation approach (Figure 3). Further studies are required to explore the linkage between ARGs and possible hosts using a combination of proximity ligation data with shotgun sequencing data, as explained above (Stalder et al., 2019) and further more by isolation.