Biosolids for safe land application: does wastewater treatment plant size matters when considering antibiotics, pollutants, microbiome, mobile genetic elements and associated resistance genes?

Soil fertilization with wastewater treatment plant (WWTP) biosolids is associated with the introduction of resistance genes (RGs), mobile genetic elements (MGEs) and potentially selective pollutants (antibiotics, heavy metals, disinfectants) into soil. Not much data are available on the parallel analysis of biosolid pollutant contents, RG/MGE abundances and microbial community composition. In the present study, DNA extracted from biosolids taken at 12 WWTPs (two large-scale, six middle-scale and four small-scale plants) was used to determine the abundance of RGs and MGEs via quantitative real-time PCR and the bacterial and archaeal community composition was assessed by 16S rRNA gene amplicon sequencing. Concentrations of heavy metals, antibiotics, the biocides triclosan, triclocarban and quaternary ammonium compounds (QACs) were measured. Strong and significant correlations were revealed between several target genes and concentrations of Cu, Zn, triclosan, several antibiotics and QACs. Interestingly, the size of the sewage treatment plant (inhabitant equivalents) was negatively correlated with antibiotic concentrations, RGs and MGEs abundances and had little influence on the load of metals and QACs or the microbial community composition. Biosolids from WWTPs with anaerobic treatment and hospitals in their catchment area were associated with a higher abundance of potential opportunistic pathogens and higher concentrations of QACs.

Introduction

Organic wastewater solids or biosolids are not only valuable organic fertilizers improving soil characteristics but also one of the main anthropogenic sources of chemical pollutants, resistance genes (RGs) and mobile genetic elements (MGEs) in agricultural soil (Singh and Agrawal, 2008; Rizzo et al., 2013; Jechalke et al., 2014; Verlicchi and Zambello, 2015; Bondarczuk et al., 2016; Wolters et al., 2019). Here we define biosolids as treated sewage sludge (Spinosa and Vesilind, 2007). Several studies revealed the presence of substances such as heavy metals, residues of antibiotics or disinfectants in biosolids that are known to exert selective pressure on bacteria (Baker‐Austin et al., 2006; Pal et al., 2015; Wales and Davies, 2015). Although RGs naturally occur in soil (Nesme et al., 2014; Nesme and Simonet, 2015) application of biosolids to soil not only contributes to changes in the soil microbial community composition, due to the addition of biosolids microbiome, nutrients, and pollutants, but also to an increased antibiotic RG load (Wolters et al., 2019). Due to the presence of microbial communities, originating from different sources, together with selective compounds and a high nutrient availability, facilitating metabolic activities and high cell densities, wastewater treatment plants (WWTPs) are regarded as hot spots of horizontal gene transfer. Therefore, WWTP sewage sludge is an important potential transfer gateway between clinical and environmental bacteria (Zhang et al., 2011; Su et al., 2015; Jacquiod et al., 2017). The contribution of biosolid application to the spread of antibiotic RGs in soil has been previously recognized (Rizzo et al., 2013; Bondarczuk et al., 2016). Once introduced into soil, the fate and bioavailability of antibiotics, disinfectants and metal pollutants depend on their respective physicochemical properties, soil characteristics and other factors, including biotic parameters such as the associated soil microbial communities or the crops grown (Lozano et al., 2010; Wuana and Okieimen, 2011; Jechalke et al., 2014; Yang et al., 2014). It has been shown that disinfectants, metals and antibiotics can enhance horizontal gene transfer at sub‐inhibitory concentrations (Kim et al., 2014). Notably, several non‐antibiotic compounds or water disinfection products such as chlorine have been shown to highly enhance transformation efficiency, thereby promoting competence in exposed cells, via increased membrane permeability due to overproduction of reactive oxygen species (Jin et al., 2020; Wang et al., 2020).

While antibiotics and disinfectants might be degraded or sequestered, and thus their effects may vanish over time, heavy metals are persistent and usually accumulate in the top soil layer (Yang et al., 2014). This is especially worrisome in light of observations from soil microcosms that metals can remain bioavailable for longer than antibiotics, such as tetracycline, and may therefore pose stronger selective pressure on antibiotic‐resistant bacteria than the antibiotics themselves (Song et al., 2017). Disinfectants can also accumulate in wastewater irrigated soils over time (Heyde et al., 2021). Thus, it is highly relevant to determine the presence of selective pollutants (heavy metals, antibiotics and disinfectants) and the abundance of RGs and MGEs in biosolids applied to soil. The use of biosolids as organic fertilizer is contentious – on the one hand, they provide valuable fertilization while on the other they introduce pollutants that probably affect the soil resistome and its transferability. Thus Murray et al. (2019) reported that root vegetables grown in soil amended with biosolids displayed higher antibiotic RG abundance compared to those grown in untreated soils.

Recently, regulations were implemented that made phosphorus recovery (>20 g kg−1 of dry matter, DM) and incineration mandatory for biosolids from larger sewage treatment plants by the Ordinance on the Reform of Sewage Sludge Utilization established by the German Ministry for the Environment (Bundesministerium für Umwelt, 2017). These regulations will come into force by 2029 for WWTPs with a catchment scale of more than 100 000 inhabitant equivalents (IEs) and by 2032 for WWTPs with more than 50 000 IEs. However, biosolids from smaller WWTPs can still be used as fertilizer, although there are no studies showing that biosolids from smaller WWTPs present lower environmental risks. Thus, these risks should be properly assessed in order to inform European Union policies like Directive 86/278/EEC and its amendments currently regulating the usage of sewage sludge in agriculture in the EU (European Commission, 1986).

Here we report on the comparison of RG and MGE copy number in DNA extracted from biosolids (biosolid DNA) obtained from two large‐scale (AL, CL) six medium‐scale (BM, EM, FM, H, JM, KM) and four small‐scale WWTPs (DS, GS, IS, LS; see Table 1) determined by means of quantitative PCR (qPCR). We also measured concentrations of pollutants associated with the measured RGs, namely: antibiotics, heavy metals and the disinfectants triclosan, triclocarban and quaternary ammonium compounds (QACs). In addition, the composition of microbial communities present in the biosolids was analysed via Illumina sequencing of 16S rRNA gene fragments amplified from biosolid DNA. Here, we aimed to elucidate whether biosolids pollutant load is influenced by WWTP size, which is typically linked to the treatment process. Anaerobic sludge stabilization is used in larger WWTPs while aerobic treatment processes are employed by smaller WWTPs. We observed a correlation between pollutant concentrations in biosolids and the biosolid resistome and mobilome, resulting in higher abundance of MGEs and RGs in biosolids containing higher concentrations of selective substances. This information is crucial for proper evaluation of current regulations on sewage sludge disposal on soil. Risks, such as dissemination of RGs and potential pathogens, together with benefits regarding their usage as fertilizer need to be evaluated.

Table 1

Size ranges of WWTPs in terms of inhabitant equivalents (IEs), state of biosolids, presence of hospital in catchment area and disposal methods for biosolids.

WWTP	IEs a	Size	WWTP treatment	State of biosolids	Hospital in catchment area	Biosolid disposal	Other features of catchment area
AL	>100 000	Large	Anaerobic digestion	Dewatered	Yes	Land use	Hospital
BM	<50 000	Medium	Anaerobic digestion	Liquid	Yes	Land use	Hospital, food industry
CL	<100 000	Large	Aerobic stabilization	Liquid	No	Incineration	Food industry
DS	<10 000	Small	Aerobic stabilization	Liquid	No	Land use
EM	<50 000	Medium	Aerobic stabilization	Liquid	No	Land use
FM	<50 000	Medium	Anaerobic digestion	Dewatered	Yes	Land use	Hospital
GS	<10 000	Small	Aerobic stabilization	Liquid, non‐stabilized b	No	Land use b
HM	<50 000	Medium	Aerobic stabilization	Liquid	No	Incineration
IS	<10 000	Small	Aerobic stabilization	Liquid	No	Land use
JM	<50 000	Medium	Aerobic stabilization	Liquid	No	Land use
KM	<50 000	Medium	Aerobic stabilization	Liquid	No	Land use
LS	<10 000	Small	Aerobic stabilization	Liquid	No	Land use	Amino‐acid production factory

Due to confidentiality agreement, only selected size ranges (<10 000; <50 000; <100 000; >100 000) are given. WWTP scale is categorized based on catchment size as follows: >50 000 IEs are categorized as ‘large’; between 10 000 and 50 000 as ‘medium’; <10 000 as ‘small’. The superscript letters (L, M, S) denote the classification.

Biosolids from this WWTP are further treated in a larger WWTP included in our study and before application to land.

Results

Resistance genes and mobile genetic elements in biosolids' DNA

All targeted genes [sulfonamide RGs sul1 and sul2, tetracycline RGs tet(A), tetA(P), tet(M), tet(W), tet(Q), streptomycin RGs strA and aadA, fluoroquinolone RG qnrS, and QAC RG qacE/qacEΔ1] were detected at quantifiable levels via qPCR in biosolid DNA obtained from each WWTP, with the exception of tet(A) in biosolid DNA from WWTP HM and tetA(P) in WWTP GS (Fig. 1). Normalized log10 transformed abundances of copy number per kg of dry matter (kg DM) ranged from 9.32 ± 0.19 to 13.88 ± 0.08 for qnrS and sul2 respectively (Fig. 1, Fig. S1). Also, most of the targeted MGE‐specific genes were detected in all biosolids (Fig. 1). qPCR with primers specific for korB and trfA (IncP‐1 and IncP‐1ε plasmids respectively) and intI1 genes (class 1 integrons) indicated their high abundance in all biosolids (Fig. 1; Fig. S1). Class 2 integrons (intI2) were detected only in dewatered sludge samples at an abundance (log10 copies per kg of DM) of 11.70 ± 0.20 in WWTP AL and 12.63 ± 0.12 in WWTP FL (Fig. 1; Fig. S1). IncI1 plasmids were detected in biosolid DNA originating from WWTPs DS, EM, GS, IS and JM with abundances (log10) ranging from 9.64 ± 0.26 to 10.03 ± 0.10 copies per kg of DM. In DNA from dewatered biosolid samples copies per kg of DM for IncP plasmid specific genes korB and trfA were one to two orders of magnitude lower compared to liquid biosolids. In contrast, the class 2 integrons integrase gene intI2 was only detectable in dewatered samples, compared to liquid biosolids where intI2 was below detection level. A higher abundance of tetA(P) was observed in biosolids obtained from WWTPs with anaerobic sludge processing and hospitals in their catchment area (AL, BM, FM) than in other biosolids (Wilcoxon rank‐sum test, p = 4.6 × 10−10) and this effect is still significant when tested only within liquid sludge samples (p = 2.4 × 10−5). The target genes specific for plasmids of the incompatibility groups IncF (traI) and IncI2 (traI) as well as LowGC plasmids (traN) were not detected.

Fig. 1

Heat map displaying abundance of target gene biosolids. Colour scale displays the mean values (n = 4) of target gene copy numbers per biosolids kg of dry of mass in each WWTP. Values below detection limit are indicated in grey. The letter ‘h’ indicates the presence of a hospital in the WWTP's catchment area.

The results of the semi‐quantitative detection of the plasmid groups PromA, IncW and IncN by PCR‐Southern blot hybridization in the biosolid‐DNA are given in Table 2. PromA plasmids were detected with strong signals even after only short exposure time (10 min) in all biosolid‐DNAs analysed, with a tendency for stronger signals occurring in DNA obtained from liquid biosolids. IncW plasmids were detected only in DNA from biosolids taken from WWTPs ES, FM and LS while plasmids belonging to the IncN group were detected after a long exposure time with very weak to weak hybridization signals indicative of low abundance in all biosolids, except for those obtained from WWTPs AL and KM.

Table 2

Detection of sequences specific for transferable plasmids of the incompatibility groups PromA, IncN and IncW in environmental DNA of biosolids via PCR and Southern blot hybridization (exposure time: PromA: 10 min, IncW: 1 h, IncN: 2 h 30 min; −: not detected, (+): weak signal, +: positive signal, ++: strong signal, +++: very strong signal).

Inc Group	AL	BM	CL	DS	EM	FM	GS	HM	IS	JM	KM	LS
PromA	+	+++	++	+++	+++	+	+++	++	++	+++	++	+
IncN	−	(+)	(+)	+	+	+	+	(+)	+	(+)	−	(+)
IncW	−	−	−	−	+	+++	−	−	−	−	−	+

Concentrations of heavy metals, antibiotics and disinfectants in biosolids

The concentrations of heavy metals, antibiotics and disinfectants determined per kg DM are summarized in supplemental data Table S2. The measured concentrations of heavy metals per kg DM ranged from 335 μg kg−1 DM for Hg to 1302.5 mg kg−1 DM for Zn (supplemental data Table S2a). Thallium (Tl) was detected at very low level in WWTPs CS and EM. In all biosolids analysed in this study, no cephalosporins, penicillins or erythromycin antibiotics were detected. The concentrations of the detected antibiotics in the biosolids (per unit DM) ranged from 1.68 μg kg−1 to 6.75 mg kg−1 for roxithromycin and ciprofloxacin respectively (supplemental data Table S2b). The highest concentrations of triclosan and triclocarban were 2.59 mg kg−1 DM and 0.11 mg kg−1 DM respectively. Triclocarban was exclusively detected in dewatered biosolids (supplemental data Table S2c). QACs were detected in most biosolids analysed (supplemental data Table S2d). The measured concentrations of QACs ranged from 0.4 μg kg−1 for chlormequat in biosolids from WWTPs HM, IS and JM up to 6.94 mg kg−1 for BAC‐C12 in biosolid from WWTP AL. The concentrations of Zn, Pb, Cd, As, triclosan, triclocarban, vancomycin, alkyltrimethylammonium compounds (ATMACs), benzyldimethylammonium compounds (BACs) and benzethonium per kg/DM tended to be higher in dewatered than liquid biosolids. Furthermore, for Cr, Pb, Ni, Cu, vancomycin, ATMACs, BACs, dialkyldimethylammonium compounds (DADMACs) and benzethonium, higher concentrations were found in biosolids obtained from WWTPs AL, BM and FM using anaerobic sludge stabilization that had hospitals within their catchment area than in biosolids from other WWTPs (Fig. S2, supplemental data Table S2d).

Correlation between contaminant concentrations, normalized gene abundances per kg of dry matter in biosolids and WWTP catchment size

As the target compounds CrVI, Tl, triclocarban, clindamycin sulfoxides, trimethoprim, sulfomethoxazol and vancomycin were undetected in most samples (supplemental data Table S2a‐d), these were excluded from the correlation analysis. Furthermore, we used total concentrations of different QAC homologues for correlation analysis. Positive and significant (p < 0.05) correlations between WWTPs size (number of inhabitants equivalent) and pollutants could only be verified for the heavy metals Cd and Cr (Fig. 2). Concentrations of heavy metals Pb, Cu, Ni, Zn and As, together with triclosan, were not significantly correlated to the number of IEs. For antibiotics, the number of IEs was negatively correlated with concentrations of all antibiotics and all but one is significant (p < 0.05). Carbamazepine concentrations were not significantly correlated with the number of IEs.

Fig. 2

Correlogram displaying all correlations between measured environmental variables in liquid biosolid samples. Spearman's rank correlation coefficients between measured antibiotics, metals and chemical pollutants, resistance genes, MGEs and wastewater catchment size are shown and only significant correlations (p < 0.05) are coloured. All values used in correlations are expressed in mass or gene copies per kg of dry matter of biosolids, except for WWTP catchment size expressed in inhabitant equivalents.

While no significant correlations (Fig. 2) were revealed with the number of IEs, the determined QAC concentrations were significantly higher for biosolids from anaerobic stabilized biosolids that originated from WWTPs with a hospital in the catchment (Fig. S2). The abundances of all remaining RGs, except for tet(A), tet(M), tet(W) and tetA(P) in the biosolids were negatively correlated with the number of IEs (Fig. 2). Furthermore, MGE target genes specific for IncP‐1, the subgroup IncP‐1ε plasmids and integrase gene intI1 were all negatively correlated to the number of IEs. As the target genes for intI2 integrase and IncI1 plasmids were not or only sporadically detected in the biosolids, these genes were excluded from the analysis.

Correlation between contaminant concentrations and RG and MGE abundance in biosolids

Significant (p < 0.05) correlations were found between concentrations of heavy metals and abundance of the targeted RG and MGE (Fig. 2). Specifically, the concentration of Cu was significantly and positively correlated with the abundance of rrn, sul1, tetA(P) and qacE/qacEΔ1 as well as intI1. Zn was significantly and positively correlated with the abundance of sul1, strA, aadA, qacE/qacEΔ1, all tested tetracycline RGs, and MGE markers intI1 and trfA. Among these, the correlation coefficients between Zn and sul1, tetA(P) and qacE/qacEΔ1 were strong (|ρ| ≥ 0.7). For some heavy metals, i.e. Pb, Cr and As, negative and significant correlations with RG abundance were observed [e.g. with tet(Q) or aadA].

Interestingly, all RG and MGE marker genes were significantly and positively correlated with the concentration of the antibiotics levofloxacin and ciprofloxacin (two fluoroquinolones) and with the tetracycline antibiotic doxycycline but not for RG tet(A) and tet(M). The macrolide antibiotics roxithromycin and clarithromycin were significantly correlated with abundances of aadA and qnrS RGs and the intI1 integrase gene (Fig. 2). In contrast, the antiepileptic carbamazepine was negatively correlated with trfA, aadA, tet(M) and tet(Q) (Fig. 2).

Significant correlations were also revealed between QAC concentrations and some target genes (Fig. 2). Especially tetA(P) was significantly and positively correlated to the concentrations of all analysed QACs. Also gene copy numbers of rrn, sul1, sul2, tet(M) and tet(W) but also MGE marker genes korB and intI1 were significantly and positively correlated to the concentrations of QACs. The QAC RG qacE/qacEΔ1 was, however, not correlated to any of the measured concentrations of QACs (Fig. 2).

Microbial community composition in biosolids

We identified 2104 unique amplicon sequence variants (ASVs) and 859 670 reads after curation, across samples corresponding to 87.9% of the input dataset reads (per sample: min. 6471; max. 31 782; mean 18 026). Relative abundance aggregated at the family level for each WWTP is reported in supplemental data Table S3. Rarefaction analysis showed that the number of reads obtained from each sample sufficiently covered their diversity (Fig. S3). Mean species richness (total observed ASVs) in biosolids of each WWTP, based on four replicate samples per WWTP, ranked from 120 ASVs (SD = 20.9) in WWTP FM to 416 ASVs (SD = 23.3) in WWTP KM (Fig. S4). The Shannon's diversity index, accounting for species richness and evenness, varied between 3.38 in WWTP AL (SD = 0.14) and 5.26 in WWTP KM (SD = 0.05). The class of the WWTP catchment size (large, medium, small) had no influence on estimated diversity indices when considering the large variation within a size class and that only two large WWTPs were sampled. Medium‐sized WWTPs had significantly higher Shannon's diversity than large or small WWTPs (Fig. S5).

When using weighted UniFrac distance to compare different biosolid microbial communities, a clear clustering of biosolid samples based on their WWTP of origin was observed in the ordination plot (Fig. 3A). A second pattern emerging was the clustering of samples based on anaerobic sludge stabilization or the presence of a hospital or food industry in the catchment area, indicating a strong influence of these parameters on the recovered community (Fig. 3A, Fig. S6). The significance of this clustering (p‐value <0.001; R 2 = 23.6%) was verified using permutational ANOVA (Dixon et al., 2003) with 10 000 random permutations.

Fig. 3

Dissimilarity and composition of the WWTP biosolid samples' microbial communities.

A. PCoA ordination of weighted UniFrac distance between all samples. Principal coordinate analyses based on weighted UniFrac distance using between samples at the ASV level. Ellipses show the 95% confidence intervals delineating the two groups of samples from WWTPs with a hospital or food industry in their catchment area (triangles) or without (circles). Percentages display the amounts of variance explained by each axis. The clustering in the two groups is supported by permutational ANOVA (p‐value <0.001; R 2 = 23.6%) based on 10 000 re‐samplings.

B. Stacked barplot of family relative abundance of WWTP biosolids (mean of four replicates). Families below 1% relative abundance are grouped in ‘Other’. For both panels, the letter ‘h’ after the WWTP name indicates the presence of a hospital in the WWTP's catchment area.

Most taxa identified in the biosolid samples were typical of WWTP communities. Most treatment plants were dominated by filamentous Actinobacteria identified as Candidatus Microthrix sp. of the Microthricaceae family (Fig. 3B). Candidatus Microthrix sp. was the most prevalent and abundant taxa identified across the dataset, although its relative abundance varied a lot between WWTP biosolids (max. 34.3% WWTP GS; min. 2.7% WWTP BM) and was not detected in WWTPs with anaerobic sludge stabilization AL, CL and FM. The WWTP AL displayed a very large relative abundance of genus Coprothermobacter from the Coprothermobacteraceae family (mean = 41.36%; SD = 4.49%) that was not detected in any other WWTP (Figs 3B and 4).

Fig. 4

Heatmap displaying the log10 transformed relative abundances of all genera with a mean >2% in any WWTP biosolid samples (n = 4 per WWTP). The amount of reads represented in this heatmap accounts for 69.5% of all reads in the dataset. The colour scale displays log10 transformed relative abundance for each sample. A pseudo‐count equal to the lowest non‐zero relative abundance in the table divided by 2 was added to zero values before transformation. The samples are clustered by Ward algorithm using the weighted UniFrac distance matrix used in principal coordinate analysis. The letter h behind the WWTP name indicates the presence of a hospital in the catchment area. Row names indicate first the phylum and then the lowest taxonomic identification, up to the genus level.

When considering all genera with a mean relative abundance above 2% in any WWTP (hereby defined as the core microbiome), the same pattern was observed and samples were grouped first by origin (WWTP) and then by the treatment process (anaerobic vs. aerobic treatment – Fig. 4). The size of the WWTP catchment did not completely explain the grouping but all small‐scale plants had a more similar microbial community than Medium and Large ones. The heatmap shows that dewatering treatment had a strong influence and these samples (WWTPs AL and FL) were specifically enriched with members of the Dysgonomonadaceae family, as recently described (Chaumeil et al., 2019), based on metagenome‐assembled genome from an anaerobic digester sludge sample. Paeniglutamicibacter, Caldatribacterium, Sporosarcina, Pedobacter and Brevundimonas tended to show higher abundance in dewatered and anaerobically stabilized biosolid samples (Fig. 4).

We identified 84 ASVs belonging to five ‘high‐risk genera’ (see methods): Acinetobacter, Enterococcus, Pseudomonas, Stenotrophomonas and Streptococcus that could potentially pose substantial health risks, either as opportunistic pathogens or efficient spreader of mobile RGs A detailed species identification of these 84 ASVs is available in supplemental data Table S4. These taxa were abundant in biosolid samples (Fig. 5A) and their overall abundance was significantly higher in WWTPs with a hospital or food industry in their catchment area (p‐value <1e−6; Wilcoxon test; Fig. 5B). All displayed differential abundances determined by the metagenomeSeq algorithm were statistically significant (p < 0.01) after false discovery rate (FDR) adjustment for multiple testing of p‐values (23 out of 84 ASVs – Fig. 5C). Interestingly, we found that the response was not genus wide as ASVs from the same genus could either had significantly increased or decreased relative abundances. However, different ASVs belonging to the same species showed a common response as the relative abundance of two ASVs affiliated with Acinetobacter indicus increased while it decreased for three ASVs belonging to Acinetobacter johnsonii decreased. Overall, most high‐risk species were enriched by more than 10‐fold in biosolids of WWTPs from a catchment area with a hospital, using anaerobic sludge stabilization, or a food industry associated with high emission of QACs.

Fig. 5

Normalized abundance of high‐risk genera in each sample and differential abundance testing at species level.

A. Stacked barplots of 16S rRNA gene qPCR corrected abundance of ASVs per kg of dry matter in each individual biosolid sample. Only members from the displayed five genera were identified in this dataset.

B. Boxplot of the sum of these normalized abundances for each individual and a Wilcoxon test show that this normalized abundance in biosolid samples is significantly different if the sample is taken from a WWTP with either a hospital or a food industry in its catchment area.

C. Differential abundance testing of manual BLAST curated species abundance using metagenomeSeq algorithm. Positive logFC values indicate the species is found in higher abundance in biosolid samples from WWTPs with a hospital or a food industry in their catchment area, and negative values indicate the opposite. For each risk ASV (rASV) the identified species, phylum and FDR adjusted p‐value are displayed (all adjusted p‐values <0.01).

Discussion

The present study provides unique insight into the pollutant load (antibiotics, heavy metals and disinfectants), the microbial community composition and the abundance of RGs and MGEs in biosolids that can be used as organic fertilizers (except for CL and HM). A selection of antimicrobial RGs and MGEs found in WWTPs or involved in RG dissemination between the environment and humans was analysed (Heuer et al., 2012; Blau et al., 2018; Liu et al., 2019; Wolters et al., 2019). The results of this study underline the role of biosolids as important sources of RGs and MGEs together with potentially selective pollutants. Our results indicate that dewatering of sewage sludge (AL, FL), anaerobic digestion (AL, BM, FL) and the presence of a hospital (AL, BM, FL) or a food industry (BM, CL) in the catchment area had a strong influence on the microbial community composition and abundance of some of the tested RGs and MGEs in the biosolids. However, the experimental design does not allow us to disentangle the WWTP specific processes, such as sludge stabilization, and presence of a hospital in the catchment.

In particular, the hosts of IncP‐1 plasmids seemed to be affected by dewatering as they were substantially less abundant in dewatered rather than in liquid biosolids. This might indicate that the respective bacterial hosts were outcompeted by other taxa or damaged during the dewatering and storage of biosolids.

PromA plasmids were shown to have a broad host range and effectively move between different bacteria via retro‐transfer (van Elsas et al., 1998). Thus, although not carrying RGs, PromA plasmids are considered important facilitators of gene spread via mobilization, e.g. of IncQ plasmids in bacterial communities (Tauch et al., 2002; Heuer and Smalla, 2012). Accordingly, these plasmids may further foster the gene spread between biosolids and soil or plant‐associated bacteria. PromA plasmids were detected in all biosolid‐DNA and IncN plasmids in 10 out of the 12 WWTPs tested by PCR–Southern blot hybridization. PromA, IncN and IncP‐1 all have a broad host range and might serve as shuttles for RG dissemination between species. This is the first study reporting the potentially concerning prevalence of PromA plasmids in WWTP biosolids, because PromA plasmids possess mobilization and self‐transfer capabilities potentiating transfer of resistance carrying plasmids (Werner et al., 2020). Presence of IncW and IncN plasmids was reported in previous WWTPs studies (Akiyama et al., 2010; Hutinel et al., 2021). IncW plasmids have been isolated from different species of Enterobacteriaceae together with IncN plasmids and are of medical significance because many can confer resistance to clinical antibiotics such as carbapenems and third‐generation cephalosporins (Fernández‐López et al., 2006; Eikmeyer et al., 2012). However, the evaluation of the load of accessory genes on these plasmids or on IncP‐1 plasmids, known to carry antibiotic and mercury RGs (Smalla et al., 2006; Schlüter et al., 2007), requires either capturing of these plasmids in recipients or cultivating their hosts. The capture approach was employed in a parallel study using the protocol described by Wolters et al. (2015) and interestingly most of the exogenously captured plasmids from biosolids belonged to IncP‐1 or IncN groups. Plasmids exogenously captured directly from biosolids typically conferred multiple resistance, e.g. towards combinations of sulfadiazine, tetracycline, streptomycin and trimethoprim (Hauschild et al., unpublished). Previous studies proposed that the diversity of IncP‐1 plasmids was driven by gene cassettes inserted into the class 1 integrons (Heuer et al., 2012; Wolters et al., 2015), and indeed all IncP‐1ε plasmids captured carried intI1, sul1 and qac genes, typical of clinical class 1 integrons.

For several monitored selective pollutants, dewatered biosolids tend to contain higher concentrations of most heavy metals and triclosan in comparison to other biosolid samples. Furthermore, triclocarban and the antibiotic vancomycin were exclusively detected in dewatered biosolids. These effects are probably due to pollutant concentrations during the dewatering process. Except for Cd in sample AL, all measured heavy metal concentrations were below the legal limit values for the reuse of biosolids on land (Collivignarelli et al., 2019). For the other analysed substances, no such legal limits exist in Germany (Bundesministerium für Umwelt, 2017).

Biosolids originating from WWTPs with anaerobic sludge stabilization that had also received wastewater from hospitals contained higher levels of QACs and higher relative abundances of the tetA(P) RG typically occurring in Clostridia. The recently developed qPCR system targeting tetA(P) was tested in the present study as this gene was proposed as an indicator for soils treated with organic fertilizers (Blau et al., 2019). Interestingly, doxycycline concentration in biosolids was positively correlated to tet(Q), tet(W), and tetA(P), while it was not significantly correlated with tet(A) or tet(M). Also, Song et al. (2017) reported in a soil microcosm study that the antibiotic tetracycline, even spiked in unrealistically high concentrations (up to 100 mg kg−1 soil), did not enhance the bacterial community tolerance against tetracycline compared to controls. Strikingly, in their study they showed instead that the heavy metals Zn and Cu were bioavailable longer than tetracycline and were able to increase the bacterial community tolerance against tetracycline. Unfortunately, their study did neither determine nor quantify the type of tet‐genes that caused the resistance observed. We show in our results that abundance of Zn and Cu are correlated with higher relative abundance of tetracycline RGs (Fig. 2). Direct selection by tetracyclines in soil is indeed unlikely since tetracyclines are quickly adsorbed to soil particles (Liu et al., 2020) reducing their bioavailability and accessibility. It is however likely that the resistance mechanisms to toxic levels of Zn, Cu or tetracycline are the same and mediated by efflux pumps, explaining the correlation.

Although the integrase gene of clinical class 1 integrons was proposed as a proxy of diverse anthropogenic pollution (Gillings et al., 2015), the abundance of intI1 gene did not show strong correlation with any pollutants concentrations in biosolids, except for Zn and Cu. The int1 gene was, however, strongly correlated with the typical genes also present on clinical class 1 integron: qacE∆1/sul1, aadA, qnrS and strA [or aph(6)‐la gene] genetic platform. This suggests that many different substances can promote the abundance of class 1 integrons and associated RGs by co‐selection.

To our surprise, hardly any correlation was revealed between any of the tested QACs and the abundance of specific QAC RG qacE/qacEΔ1. While these results might be impacted by the challenging QACs trace analysis, it is also conceivable that the level of resistance conferred by the incomplete qacEΔ1 is low, as reported by Paulsen et al. (1993). The truncated gene cannot be distinguished from the full gene with the qPCR assay used. The observed significant correlations between QACs and some targeted RGs might be due to potential co‐selection via other QAC genes such as, i.e. qacA, qacB or qacC, as mentioned in Mulder et al. (2018).

The high abundance of sequences assigned to Coprothermobacter exclusively in biosolids from WWTP AL is because this was the only WWTP that performed thermophilic anaerobic digestion of biosolids. Coprothermobacter is strictly anaerobic and often abundant in microbial communities present in thermophilic anaerobic digesters throughout the globe, thus providing optimal conditions for Coprothermobacter bacteria with optimal growth temperatures between 55°C and 70°C (Gagliano et al., 2015; Kunath et al., 2019). In the context of this study, we were particularly interested in investigating the introduction of potential opportunistic pathogens via soil application of biosolids. We could show that potentially high‐risk taxa were prevalent and abundant in samples from biosolids obtained from WWTPs with anaerobic sludge stabilization. Despite the lack of proper experimental to test this it is tempting to speculate that the presence of a hospital or a food industry in the catchment area might have also contributed to the increased relative abundance of potential pathogenic bacteria or spreaders of resistance plasmids. The anaerobically treated biosolids also all have a hospital in their catchment area and a significantly higher level of QACs compared to all other WWTPs (Wilcoxon's test p‐value <3.6e‐5 – Fig. S2). This suggests a selective effect of these substances in sewage waters when combined with anaerobic treatment.

The present study focus on a chosen subset of RGs, MGEs and samples were collected only once per WWTP. Thus potential seasonal fluctuations, as described by Martín et al. (2012) were unaccounted for. It would be interesting to complement these targeted approaches with more recent holistic methods such as shotgun metagenomics (Michán et al., 2021) which could, at the cost of sensitivity, yield a more complete picture of the correlations between the total resistome and pollutants concentration. Additional information, such as the presence of heavy industry or agricultural farms in the catchment area was not available. Although nearly all measured heavy metal concentrations were below the legal limit of values for the reuse of biosolids on land, heavy metals, antibiotics and disinfectants were detected in biosolids of all tested WWTPs at considerable levels. We also showed that biosolids from small‐scale WWTPs were not less contaminated with selective pollutants, RGs or MGEs than those from large‐scale facilities. Rather, the sizes of the WWTPs in terms of IEs were negatively correlated to most target genes abundance tested in this study. This was probably influenced by the changes in the microbial community composition during anaerobic digestion. In order to minimize the introduction of RGs and MGEs and potential pathogens via biosolids into the environment, the results of the present study do not support the regulation planned for biosolids based solely on the size of the catchment. Our results suggest that biosolids from small‐size WWTPs have a considerable pollutant load that should be considered before using these biosolids for soil fertilization.

Experimental procedures

Characteristics of biosolid samples used in this study

Biosolids from WWTPs are the final and usable product obtained after the aerobic or anaerobic treatment of sewage sludge for stabilization (Spinosa and Vesilind, 2007). We refer to dewatered biosolids when most of the sewage sludge water content is removed by a mechanical process in an industrial decanter. The details of all biosolid samples used are reported in Table 1.

Sample collection and processing

Biosolids from a total of 12 different WWTPs located in the South East of Lower Saxony, Germany, were kindly provided by the WWTP staffs in spring 2018. Due to a confidentiality agreement, only selected detailed information on the WWTPs is given in Table 1.

Four samples of 4L to 6L of biosolids from each plant were received and immediately stored at −20°C before quantification of heavy metals, disinfectants and antibiotics. For extraction of biosolid DNA at the day of receipt, 45 ml of liquid biosolids from each sample were centrifuged at 4000g for 10 min, supernatants were discarded and pellets were stored at −20°C until use. Samples of pelleted (n = 40) or dewatered biosolids (n = 8) (WWTPs AS and FM) were homogenized manually with sterile spatula and 0.3 g was used for DNA extraction from biosolids using the FastDNA™ SPIN Kit for Soil (MP Biomedicals, Heidelberg, Germany). Biosolid‐DNA preparations were purified using Geneclean™ Spin Kit (MP Biomedicals). For qPCR using TaqMan‐based assays, aliquots of the purified DNA were diluted 1:10 with Tris‐EDTA buffer (10 mM Tris, 1.25 mM EDTA, pH 8.0) to dilute potential PCR inhibitors that would affect correct quantification. For conventional PCR assays, undiluted DNA was used.

Quantification of 16S rRNA genes, RGs and MGEs in DNA from biosolids

The abundances of bacterial 16S rRNA gene (rrn) copies, marker genes specific for various plasmid groups (IncF: traI; IncI1: traI; IncI2: traI; IncP‐1: korB; IncP‐1ε: trfA; LowGC: traN), class 1 and 2 integron integrase genes (intI1, intI2) and selected RGs (sul1, sul2, tet(A), tetA(P), tet(M), tet(W), tet(Q), strA, qacE/qacEΔ1) were quantified in biosolid‐DNA by 5′‐nuclease qPCR assays (TaqMan) in a CFX96 real‐time PCR detection system (Bio‐Rad, Hercules, CA, USA). Detailed information on the assays used is given in supplemental data Table S5. The tetA(P) primers and probe were developed based on an alignment of six tetA(P) sequences related to Clostridium using the program CLC Main Workbench version 8.1.3. (CP019580.1; NC_019259.1; AB054980.1; NC_010937.1; NG_048314.1; NG_048316.1). Standard dilutions for tetA(P) were generated from a cloned 131 bp PCR product obtained from a pig manure sample. Gene copy abundances per sample were calculated per kg of DM.

Detection of sequences specific for plasmid groups PromA, IncW and IncN in DNA from biosolids by PCR and Southern blot hybridization

PCR assays used for the detection of PromA, IncN and IncW plasmids were previously described and primers and references are listed in supplemental data Table S5. To increase specificity of detection, PCR products were Southern blotted and hybridized with corresponding digoxygenin‐labelled probes. Probe generation and Southern blot hybridization were performed as detailed in Götz et al. (1996) and Binh et al. (2008). PCR products obtained from reference plasmids (IncW: R388; IncN: RN3; PromA: pIPO2) were digoxygenin‐labelled using the DIG‐labelling Kit (Roche Diagnostics, Mannheim, Germany).

Determination of heavy metal, antibiotic and disinfectant concentrations in biosolids

Quantification of the heavy metals CrVI, Pb, Cd, Cr, Cu, Ni, Hg, Zn, As and TI concentrations in biosolids was determined by AGROLAB Agrar und Umwelt GmbH (Kiel, Germany). The concentration of triclosan and triclocarban in the biosolids were determined at the LUFA Nord‐West (Hameln, Germany). If applicable, corresponding DIN – (‘Deutsches Institut für Normung’ – German Institute for Standardization) and ASU – (‘Amtliche Sammlung von Untersuchungsverfahren’ – Official Collection of Test Methods) norms for each measured pollutant are given in supplemental data Table S6.

Antibiotic concentrations were determined separately in the solid and liquid phase of the biosolids at the Institute of Clinical Pharmacology, Technical University of Dresden, as previously described (Rossmann et al., 2014; Marx et al., 2015). The analysis included the following antibiotic molecules, representative for different drug classes: sulfonamides (sulfamethoxazole), tetracyclines (doxycycline), fluoroquinolones (ciprofloxacin and levofloxacin), penicillins (penicillin V, piperacillin and amoxicillin), diaminopyrimidines (trimethoprim), lincosamides (clindamycin and clindamycin sulfoxide), macrolides (erythromycin, azithromycin, clarithromycin and roxithromycin), glycopeptides (vancomycin) and cephalosporins (cefuroxime, cefotaxime). Furthermore, the analysis included the anti‐epileptic carbamazepine, an emerging contaminant frequently detected in wastewater and biosolids due to its low elimination rate (Ramaswamy et al., 2011). The respective concentrations of antibiotics in the original substances were normalized using the known DM contents of the biosolids.

The concentrations of QACs were determined as described in Heyde et al. (2020). In brief, 2 g of freeze‐dried sample was extracted in three cycles with 10 ml acetonitrile containing 0.1% HCl in each cycle in an ultrasonic bath followed by centrifugation. Matrix clean‐up of the extracts was performed with solid‐phase extraction using CN cartridges (Chromabond, Macherey Nagel, Düren, Germany). A Waters 2690 liquid chromatograph coupled to a Waters Micromass Quattro Micro triple quadrupole mass spectrometer was used for the HPLC–MS/MS separation and analysis. Mobile phases were a formic acid/sodium formate buffer (A) and acetonitrile (B) while as a solid phase a C18 CSH column was used (Waters XSelect CSH Phenyl‐Hexyl column; 3.5 μm particle size, 2.1 × 150 mm). The column temperature was set at 37°C. A flow of 0.25 ml min−1 and a 20 μl injection volume were chosen. The gradient was 0 min 100% A, 3 min 25% A, 14.5 min 0% A, 20 min 0%, 21.01 min 100%, 26 min 100% A. The method included the analysis of the following QACs: ATMACs, BACs, DADMACs of the chain lengths C8–C16 for each group as well as chlormequat and benzethonium. All measurements are reported in supplemental data Table S2.

Correlation analysis between measured environmental parameters

Spearman's rank correlation coefficients (ρ) between the size of each WWTP catchment (in IEs), abundance of genes targeted by qPCR and concentrations of heavy metals, antibiotics and disinfectants were determined by means of the base R package using the cor() function and associated p‐values with the cor.test() function. Sizes of the WWTPs’ catchment area (in IEs) were used to calculate correlation coefficients. Due to confidentiality agreement these values are not explicitly reported, WWTPs were classified as follows: ‘large’ (>50 000 IEs); ‘medium’ (10 000–50 000 IEs) and ‘small’ (<10 000 IEs). Only significant correlations for which p < 0.05 are coloured in Fig. 2. As liquid and dewatered biosolids differed obviously, the analysis was restricted to liquid biosolids to avoid spurious correlations. Analytes for which a majority of samples had values below limits of detection were excluded from the correlation analysis: genes intI2 and traI IncI‐1; hexavalent chromium (CrVI) and thallium (Tl) metals; biocide triclocarban and antibiotics clindamycin sulfate, trimethoprim, sulfamethoxazole and vancomycin. Remaining values below the respective detection levels of each pollutant were treated as zeroes (U.S. Environmental Protection Agency, 2000). The correlogram was plotted using corrplot R package (Wei, 2017).

Illumina MiSeq sequencing of 16S rRNA gene amplicons

Microbial community compositions in biosolid samples were assessed by high‐throughput amplicon sequencing. A fragment spanning the hypervariable regions V3–V4 of the 16S rRNA gene was amplified from each biosolid DNA sample by an initial PCR step using primers Uni341F (5′‐CCTAYGGGRBGCASCAG‐3′) and Uni806R (5′‐GGACTACNNGGGTATCTAAT‐3′) originally published by Yu et al. (2005) and modified as described in Sundberg et al. (2013). In a second PCR reaction step the primers additionally included Illumina specific sequencing adapters and a unique dual index combination for each sample. After each PCR reaction, amplicon products were purified using HighPrep™ PCR Clean Up System (AC‐60500, MagBio Genomics, Gaithersburg, MD, USA) paramagnetic beads using a 0.65:1 (beads:PCR reaction) volume ratio to remove DNA fragments below 100 bp in size and primers. Samples were normalized using SequalPrep Normalization Plate (96) Kit (Invitrogen, Maryland, MD, USA) and pooled using a 5 μl volume for each sample. The pooled sample library was concentrated using DNA Clean and Concentrator™‐5 kit (Zymo Research, Irvine, CA, USA). The pooled library concentration was determined using the Quant‐iT™ High‐Sensitivity DNA Assay Kit (Life Technologies). Before library denaturation and sequencing, the final pool concentration was adjusted to 4 nM. Amplicon sequencing was performed on an Illumina MiSeq platform using Reagent Kit v2 (2 × 250 cycles) (Illumina, San Diego, CA, USA).

Generation of ASVs

Raw sequence reads were trimmed of primer sequences used in the first PCR, discarding read pairs for which any of the two primer sequences could not be detected using cutadapt version 2.3 (Martin, 2011). Primer‐trimmed sequence reads were error‐corrected, merged and ASVs identified using DADA2 version 1.10.0 (Callahan et al., 2016) plugin for QIIME2 (Guerrini et al., 2019) with the following parameters: truncL = 0, truncR = 0; trimL = 8, trimR = 8, a minimum overlap of 12 bp and otherwise using default parameters. Each ASV sequence was given a taxonomic annotation using q2‐feature‐classifier classify‐sklearn module trained with SILVA SSU rel. 132 database (Quast et al., 2012). Prior to training the classifier, the V3–V4 region of SILVA database sequences was extracted at the same primer position. A multiple sequence alignment of the ASVs was performed with mafft v7.407 with –retree 1 (Katoh and Standley, 2013) and used to build an approximate ML tree with FastTree v2.1.10 and rooted at midpoint (Price et al., 2010). We used decontam R package (Davis et al., 2018) to filter ASVs identified as PCR contaminants using ‘prevalence’ method and threshold =0.5. We filtered any ASVs shorter than 380 bp, with less than five reads across all samples or not present in at least two out of four replicate samples. Relative ASV abundances weighted by total gene copies per kg of DM were obtained by multiplying each ASV relative abundance by the total rrn copy number per kg of DM determined by 16S rRNA gene qPCR in that sample. We did not apply correction using ribosomal operon number as this can introduce systematic biases that are difficult to control as recently shown in a systematic study (Louca et al., 2018). As such, these transformed abundances remain estimates but allow us to express all measurements per kg of biosolid DM.

Microbial community composition

All data analyses were performed in R v3.6.1 (R Core Team, 2019) in RStudio v1.1.463 (RStudio Team, 2018) and made use of the following packages: ggplot2 (Wickham, 2011), ggpubr (Kassambara, 2018), phyloseq (McMurdie and Holmes, 2013), DAtest (Russel et al., 2018), MicEco, COEF (github.com/Russel88/), vegan (Dixon et al., 2003), RandomForest (Liaw and Wiener, 2002), cluster (Maechler et al., 2018) and dplyr (Wickham et al., 2019). Differential abundance analysis was performed using metagenomeSeq (Paulson et al., 2013) after selection using DAtest (Russel et al., 2018).

Identification of high‐risk taxa

ASVs confidently affiliated at the genus level after manual BLAST curation (100% query cover; >95% identity) to any defined high‐risk genera were included. These correspond to Enterococcus, Staphylococcus, Klebsiella, Acinetobacter, Pseudomonas, Enterobacter, Burkholderia, Herbaspirillum, Ochrobactrum, Ralstonia, Helicobacter, Campylobacter, Salmonella, Neisseria, Streptococcus, Haemophilus, Shigella or Stenotrophomonas. These ‘high risk taxa’ were selected both based on risk assessment by the WHO and its list of Critical Pathogens (WHO, 2017) and based on a list genera containing species striving both in plant‐associated soil and humans (Berg et al., 2005). The selected genera contain a number of opportunistic pathogens and known spreaders of resistance plasmids, hence their relative high risk for public health. We quantified the abundance of high‐risk taxa using their relative abundance obtained by amplicon sequencing and the normalized copy numbers of rrn copies per kg of DM in each sample.

Sequence data availability

Raw sequencing reads used in this study are available at ENA‐SRA under BioProject accession number PRJEB36651.

Fig. S1. Copy numbers of targeted gene copies normalized per kg of DM estimated by TaqMan qPCR assays. Each panel displays values obtained for given target genes. Samples (n = 4 per WWTP) are grouped based on their WWTP of origin class size. Difference in the means of WWTP size class groups are evaluated using Wilcoxon's rank‐sum test.

Fig. S2. Concentrations of the different quaternary ammonium compounds in mg per kg of DM quantified in this study. Differences in the means of the two groups, based on the presence of a hospital in the catchment area were evaluated using Wilcoxon's rank‐sum test. All samples (n = 48) are used in the comparison. For Alkyltrimethylammonium compounds (ATMACs), benzyldimethylammonium compounds (BACs), dialkyldimethylammonium compounds (DADMACs), the sum of all compounds of different carbon chain length is used.

Fig. S3. Rarefaction curves of observed amplicon sequence variants (ASVs) for wastewater biosolid samples. Rarefaction curves for each biosolid sample based on 1000 independent rarefactions of four replicate samples per individual WWTP showing saturation of observed ASVs for all samples indicating sequencing depth is covering most of diversity across the whole dataset. Each sample rarefaction curve (n = 4 per WWTP) point is plotted using the mean of 100 independent rarefaction events for each read depth. Error bars are indicating standard deviation of this measure.

Fig. S4. Observed amplicon sequence variants (ASVs) and Shannon's diversity index of wastewater biosolid samples. Boxplots of observed ASVs (total richness) and Shannon's diversity index estimators of alpha diversity are calculated using four individual samples per wastewater biosolids. Colours indicate the catchment size treated by each WWTP. The box displays the interquartile range, the thick line within the box display the median and the whiskers display the upper and lower quartiles of the distribution.

Fig. S5. Post‐hoc comparison tests of alpha‐diversity estimators between large, medium and small scale WWTP biosolid communities. Boxplots represent comparisons of observed richness (A) and Shannon's estimator of Diversity Index (B) of samples grouped by their respective WWTP class size. Significance levels between groups are calculated using Wilcoxon rank test and display absolute p‐values.

Fig. S6. Samples clustering based on the prokaryotic microbial community (A) and the concentration of resistance genes, MGEs and pollutants (B). Dendrograms showing the hierarchical classification of samples using Ward clustering based on weighted UniFrac distance matrix calculated on normalized ASV abundances (A) or Euclidean distance matrix on centre‐scaled (z‐score) of measured concentration of RGs, MGE markers by qPCR and chemical analysis of pollutants (B). All values were normalized per kg of dry matter of sample. Leaf labels display the WWTP; the class size of the WWTP and the type of biosolid samples and are written in blue if a food industry or a hospital is present in the catchment area. Coloured rectangles display the best two clusters of samples and respectively separate samples based on the presence of clinic or food industry (A) or only presence of a clinic (B). Clusterings were tested using PERMANOVA and 10,000 re‐samplings on their respective distance matrix (p‐value <1e‐5 for both).

Fig. S7. Comparisons between biosolid samples and negative controls. (A) Hierarchical classification of samples using Ward clustering based on Bray–Curtis dissimilarity calculated on normalized ASV abundances of all samples, including negative controls. Negative controls are clearly clustering outside the biosolid samples. (B) Number of reads identified as risky ASVs belonging to each sample, including negative controls. The ASVs identified in first PCR negative control do not belong to the same genera of those identified in biosolids samples except for Pseudomonas.

Click here for additional data file.

Supplemental Data File S1. Relative abundances of target genes (per rrn) in biosolids.

Supplemental Data File S2. Concentrations of heavy metals, biocides and ABs in biosolids per Kg of dry matter.

Supplemental Data File S3. Mean taxa relative abundance and standard deviation (sd), based on four replicates in each WWTP. Only Families where abundance >1% in any sample are shown. Each line represents the agglomerated relative abundance of all ASVs belonging to each Family. When ASVs could not be assigned down to the Family level, relative abundance is given for the lowest taxonomic rank that could be confidently annotated.

Supplemental Data File S4. BLAST alignment results for each amplicon sequence variant (ASV) and corresponding species‐level identification. When several equally good alignments are produced on different database reference sequences, all possible species names and corresponding accession number on NCBI GenBank are given for that ASV.

Click here for additional data file.

Supplemental Data File S5. Quantitative real time PCR assays.

Supplemental Data File S6. DIN‐ and ASU‐ regulations for contaminant trace analytics.

Supplemental Data File S7. Absolute copy numbers of each target gene measured using TaqMan quantitative PCR assays in each DNA sample (n = 4) from the 12 WWTPs biosolids. When the targeted gene was not detected in the sample (Ct below non template negative control), the copy number value is ‘NA’.

Click here for additional data file.