Annual dynamics of antimicrobials and resistance determinants in flocculent and aerobic granular sludge treatment systems

The occurrence and removal patterns of 24 antimicrobial agents and antimicrobial resistant determinants namely 6 antibiotic resistance genes (ARGs) and 2 mobile genetic elements (MGEs), and the fecal indicator E. coli were investigated in three full-scale wastewater treatment plants. Their waterlines and biosolids lines (including secondary treatment based on both granular and activated sludge) were sampled monthly throughout one year. Samples were analyzed by means of LC-MS/MS, qPCR and cell enumeration, respec- tively. The inﬂuence of rainfall, temperature, and turbidity on the occurrence and removal of the aforementioned agents was assessed through statistical linear mixed models. Ten of the antimicrobial agents (macrolides, ﬂuoroquinolones, tetracyclines, and sulfonamides) were commonly found in inﬂuent in concentrations of 0.1-2 μg L − 1 , and the predominant ARGs were ermB and sul1 (6.4 and 5.9 log 10 mL − 1 re- spectively). Warmer temperatures slightly reduced gene concentrations in inﬂuent whilst increasing that of E. coli and produced an uneven effect on the antimicrobial concentrations across plants. Rainfall di- luted both E. coli (-0.25 logs, p < 0.001) and antimicrobials but not genes. The wastewater treatment reduced the absolute abundance of both genes (1.86 logs on average) and E. coli (2.31 logs on average). The antimicrobials agents were also partly removed, but 8 of them were still detectable after treatment, and 6 accumulated in the biosolids. ARGs were also found in biosolids with patterns resembling those of inﬂuent. No signiﬁcant differences in the removal of antimicrobials, genes and E. coli were observed when comparing conventional activated sludge with aerobic granular sludge. Irrespective of the type of sludge treatment, the removal of genes was signiﬁcantly reduced with increasing hydraulic loads caused by rainfall (-0.35 logs per (cid:2) average daily ﬂow p < 0.01), and slightly decreased with increasing turbidity (-0.02 logs per (cid:2) 1 nephelometric turbidy unit p < 0.05) .


Introduction
The occurrence of anthropogenic antibiotic resistance genes (ARGs) in the environment seems to be strongly related to fecal pollution ( Karkman et al., 2019 ). Human fecal bacteria are transported through sewage networks and wastewater treatment plants (WWTPs) to the environment. WWTPs, primarily designed for nutrient removal, have therefore been hypothesized as key vectors in the environmental dissemination of antibiotic resistance ( Rizzo et al., 2013 ). The fate of ARGs in WWTPs has been assessed in recent years, providing a wide range of results. In some cases, WWTPs reduced the absolute concentration of ARGs ( Di Cesare et al., 2016 ;Rodriguez-Mozaz et al., 2015 ), while in others, the relative and even absolute abundances of ARGs increased after treatment ( Rafraf et al., 2016 ).
There is a growing interest in determining which technologies or operational conditions achieve greater ARGs removal. How-ever, due to the wide diversity of treatment processes, it is difficult to obtain general results that are applicable across the variety of WWTP systems and ARG types ( Korzeniewska and Harnisz, 2018 ;Novo and Manaia, 2010 ). Previously, broad sampling efforts to analyze the influent and effluent of more than 60 WWTPs in The Netherlands helped to determine that the removal of ARGs was significantly smaller when the baseline hydraulic load of the WWTPs increased because of rainfall events ( Pallares-Vega et al., 2019 ). Yet, even under dry weather conditions, process design and operational parameters did not explain the remaining variability. This suggests that other factors contributed to the detected differences.
One factor that may impact the occurrence and removal of ARGs is the presence of antimicrobial and disinfectant residues. These compounds are collected in the wastewater along with the feces and might accumulate in the sewage sludge ( Gao et al., 2012 ;Göbel et al., 2005 ). Antimicrobial and disinfectant residues may enhance the co-selection of resistance genes and promote horizontal gene transfer even at subtherapeutic concentrations ( Gullberg et al., 2014 ). In addition, fluctuations in other abiotic factors such as temperature might also influence the occurrence and removal of ARGs. Higher antimicrobial consumption in colder seasons ( Coutu et al., 2013 ;Marx et al., 2015 ) might increase the discharge of resistant bacteria and their ARGs into the sewer. In addition, cold temperatures reduce the efficiency of water treatment ( Johnston et al., 2019 ), and might also reduce the ability to remove ARGs. Literature addressing the effect of temperature on the occurrence and fate of ARGs during wastewater treatment is, however, scarce. A few studies have focused on quantifying the seasonal variation of ARGs in influent and their removal ( Caucci et al., 2016 ;Jiao et al., 2018 ). Moreover, temperature or seasonality have seldom been addressed in combination with rainfall or flow ( Schages et al., 2020 ). Thus, there is a need for a comprehensive approach studying the combined effect of these variables on ARG patterns.
Furthermore, biotic factors such as the accumulation of ARGs in the natural microbiome of the biological treatment deserve attention. In the biological treatment, the microbial community, which converts wastewater compounds, is aggregated into flocs (activated sludge), biofilms, or granular sludge. Such bioaggregates might provide a suitable environment for cell-to-cell interactions and genetic exchange of mobile genetic elements (MGE) containing ARGs, ( Manaia et al., 2018 ) potentially increasing the so-called resistome of the sludge. Diverse sludge types (flocculent, granules) might accumulate ARGs differently, derived from their characteristic physical structure and microbial community. Ultimately, the intrinsic sludge resistome may counterbalance a further elimination of ARGs during the wastewater treatment. The accumulation of ARGs in the biosolids is also of concern, given that sludge is often processed and applied as a fertilizer due to its high content of organic nutrients and phosphorus. As such, this practice comprises another possible route for the spread of anthropogenic ARGs to the environment ( Rahube et al., 2014 ).
Generally, studies covering the effect of a broad number of variables on ARGs are restricted to few time-point measurements.
Although cross-sectional studies provide relevant information, it is necessary to investigate whether measurements at single time points are valid throughout extended periods. Such information is needed to determine suitable sampling strategies to answer specific research questions. Hence, the aim of the present study was to investigate the occurrence and fate of ARGs, MGEs, and viable fecal bacteria over one year of operation at three different fullscale municipal WWTPs located in The Netherlands. These treatment plants performed biological nutrient removal, with three systems based on activated sludge and one system based on aerobic granular sludge. Besides gene determinants, we investigated the presence of selective agents as antimicrobial compounds and disinfectants in both the water and biosolids line. Analyses were aimed at determining the role of abiotic parameters such as the hydraulic load factor, seasonal temperature , and the effluent turbidity (as a surrogate for effluent TSS) on the gene removal capacity of WWTPs for an extended sampling period of one year, and to study the degree of temporal variability.

Characteristics of the selected WWTPs
Three Dutch WWTPs ( supplementary information Figure S1 ), of different sizes and process design, were considered. The sampling points are displayed in Fig. 1 . WWTP1 (Leeuwarden, 226'667 p.e.) processes 25 0 0 0 m 3 d −1 under dry weather conditions. The biological nutrient removal activated sludge process is operated on raw wastewater. To support the biological phosphate removal, iron (FeII) is added to the activated sludge tanks. The chemical oxygen demand load (COD) to the plant consists of 56% household wastewater and 26% industrial wastewater (from which half corresponds to a dairy industry). The catchment area includes a medium-size hospital (650 beds) with a load contribution of 1%.
WWTP2 (Harnaschpolder, Den Haag, 1'260'0 0 0 p.e) is the largest plant in The Netherlands. It treats an average of 150 0 0 0 m 3 d −1 . 84% of the COD load comes from households and 16% from industry. The catchment area includes several hospitals, totalling 2610 beds. The design of this WWTP consists of 8 identical parallel lines. Each line is composed of primary settling and a biological nutrient removal activated sludge process.
WWTP3 (Garmerwolde, Groningen, 340'146 p.e.) treats 71 800 m 3 d −1 (64% households,14% industry, and 1% hospital (totalling 1920 beds) in two separate treatment lines. Approximately 50% of the influent is treated in a two-stage activated sludge adsorption-bioxidation (AB) process ( Böhnke, 1997 ). The other half of the wastewater is treated by an aerobic granular sludge process (Nereda®). The AB system has been described in detail by De Graaff et al. (2016) . The raw influent undergoes two consecutive treatment steps. First, the organic content is removed in the highly loaded A-stage activated sludge process operated at a short solid retention time. Phosphate is removed in the A-stage by the addition of iron (Fe III). After the intermittent clarifier, the second B-stage activated sludge process is operated at a long solids retention time to allow nitrification and removal of the remaining biological oxygen demand (BOD). Nitrogen removal is limited by the low BOD content; therefore, methanol is added to promote denitrification. The parallel aerobic granular sludge line ( Pronk et al., 2015 ) includes a buffer tank to store the influent wastewater for up to 3 h before treatment, and two Nereda® reactors containing aerobic granular sludge. The activated sludge is mainly in granular form but also contains a fraction of flocs ( Ali et al., 2019 ). The HRT of this reactor is 6-h in dry weather conditions. All removal processes (BOD, denitrification, and phosphate removal) occur under the alternation of anaerobic feeding and aeration. Denitrification partly occurs as simultaneous nitrification/denitrification and partly by on/off aeration. Phosphate removal is essentially biological, with supplemental addition of iron (Fe III) only taking place during extreme rain weather flow.
In all three WWTPs, the surplus sludge is digested in a mesophilic reactor, and the digested sludge is subsequently dewatered and incinerated. In WWTP1, the digester receives sludge from other industrial physical-chemical treatment plants in the area.

Collection of samples from the wastewater line and the biosolids line
Waterline and biosolids line samples were taken every month over one year from April 2017 until March 2018 from the three aforementioned WWTPs.
For the waterline, volumes of 1 L of 24-h flow-proportional composite samples were collected in sterile plastic bottles (VWR, NL), except for the effluent of the AB line in WWTP3 that lacked an autosampler. In that case, a grab sample was collected instead. All autosamplers for the waterline collection had a refrigeration system to ensure samples were kept cold during the 24-h collection period. To account for possible daily variations, the waterline samples were collected three days in a row except for WWTP3, for which this was not possible due to technical reasons. All waterline samples were filtered for DNA and culture-based methods within 6 h after collection. Filters were stored at -20 °C upon DNA extraction. Waterline samples were backed up frozen at -20 °C for downstream chemical analysis.
For the biosolids line, samples were collected once per month as grab samples in 0.5 L sterile plastic bottles. The flocculent biomass samples (further referred to as "AS") were collected from the mixed liquors from the activated sludge tanks. The aerobic granular sludge samples (further referred to as "AGS") were taken from the purge of excess sludge. Digested sludge samples (further referred to as "DS") were taken as a grab sample from the digested sludge leaving the digestor. All samples were stored in cooling boxes and kept cold (4 °C) during transportation. Biosolids line samples were serially diluted and filtered for the culture-based method within -6 h after collection. Additionally, aliquots were backed up frozen at -20 °C until processing for the downstream chemical analysis and DNA extraction.

Filtration for E. coli enumeration
Monitoring the presence and removal of E. coli was used to evaluate whether fecal bacteria and ARGs would follow similar trends. E. coli was chosen as a surrogate for fecal indicators. For the waterline, samples were processed as previously described in Verburg et al., 2019 . For the biosolids line, 2 g of homogenous sample was re-suspended in 20 mL of the saline solution (NaCl) 0.85% (w/v). Serial dilutions were performed and filtered as indicated for the waterline.
After filtration, the resulting filters from both the water and the biosolids lines were plated on Tryptone Bile X-Glucuronide (TBX) selective media (Oxoid, Thermofisher, UK). The plates were incu-bated for 24 h at 37 °C, and The CFUs were enumerated following ISO guidelines ( ISO 89199:2005-12 ).

DNA extraction and qPCR analysis
In order to analyze the water and sludge samples for the presence of ARGs and MGEs, samples were pre-treated and extracted as follows.
For the waterline samples, volumes of 200 mL of effluent and 25 mL of influent (and samples with similar solids content) were filtered through 0.22 μm Durapore PVDF membranes (Merck-Millipore) in a Millipore-Sigma filtration system. The filters were frozen at -20 °C upon extraction. The DNA extraction for waterline samples was performed using the DNeasy kit Power Water from (Qiagen, NL) following the manufacturer's instructions.
For the biosolids line samples, 0.50 g of all types of AS, 0.05 g of DS and 0.10 g of AGS samples were extracted according to the MiDAS Field Guide to the Microbes of Activated Sludge and Anaerobic Digesters, versions 7.0 (AS), and 1.0 (DS) (ref: http: //www.midasfieldguide.org ) with small modifications: the beadbeating step was performed at 6800 rpm in a Precellys homogenizer (Bertin Technologies SAS, FR) which is equivalent to a speed mode of 5.5 in the FastPrep homogenizer. The final elution was reduced to 100 μL to achieve more concentrated DNA extracts.
After the extraction, the DNA extracts were diluted 10 or 100fold to avoid inhibition and ensure the target was within the range of quantification. The diluted DNA was analyzed by quantitative polymerase chain reaction (qPCR). The panel of studied genes ( Table 1 ) included six ARGs ( sul1, sul2, tetM, qnrS, ermB, and bla CTX-M ), two MGEs ( intI1 gene and korB gene). The 16S rRNA gene acted as a surrogate for total bacteria. The reactions were prepared and performed as previously described in Pallares-Vega et al. (2019) with small modifications for the korB assay. This information is specified in Appendix A .

Total solids content and antimicrobial and disinfectant compounds in water and biosolids samples
The total solids (TS) content of the sludge samples was determined by standard methods ( Clesceri et al., 1998 ). This information was required to express the results from the microbiological, molecular, and physicochemical analyses in biosolids as normalized to the total solids (TS).
The determination of antimicrobial and disinfectant residues in waterline and biosolids lines samples was performed using liquid chromatography-tandem mass spectrometry (LC-MS/MS). The analyzed compounds are compiled in Table 1 . Pretreatment and analysis of the waterline samples were performed, as described in Verburg et al. (2019) . Biosolids sample preparation, specifications of the device, and run are summarized in Appendix B. For each sample, a surrogate sample spiked with the known concentrations was used to calculate the recovery.

Sampling factors and statistical analysis
Information about abiotic factors during sampling collection was obtained as follows: data about daily flow and average annual flows for 2017 and 2018 were obtained from the WWTP operators. With this data, the hydraulic load factor (HLF) was calculated. This parameter stands for the ratio of the flow ( i.e. volume of water treated) on the day of sampling divided by the average daily flow (derived from the annual flow) of each WWTP ( Pallares-Vega et al., 2019 ).
Air temperature on the day of sampling was retrieved from https://weerstatistieken.nl . Turbidity in effluent samples was analyzed by means of a turbidimeter (2100 N IS, Hach). The influence of these factors in the incoming and removal of genes and E. coli, as well as the role of intermediate steps in the removal of genes, were analyzed by linear models and linear mixed models, further described in detail in Appendix C.

Antimicrobials, ARGs and E. coli in the influent: role of rain dilution and seasonal temperature
From the 24 antimicrobials included in this study, 10 were detected in the influent of the three WWTPs, in general, in concentrations within the ng L −1 scale ( supplementary information Table S1 ). Overall, the values observed in this study were consistent with several other works across different countries . The fluoroquinolone ciprofloxacin was the most abundant compound in the influent, consistently exceeding the predicted no inhibitory concentrations (PNEC) described by ( Bengtsson-Palme and Larsson, 2016 ) for selection of antimicrobial resistance for this compound (0.064 μg L −1 ), with maximum levels reaching 1.20 μg L −1 . The macrolides azithromycin and clarithromycin, and the therapeutic group of sulfonamides (sulfapyridine sulfamethoxazole)-trimethoprim were also commonly detected in influent and sometimes exceeding their corresponding PNEC levels. In contrast, tetracyclines (doxycycline and tetracycline), which are the second most consumed antibiotics, were prevalent in the influent but below their PNEC levels (2 and 1 ug L −1 respectively).
The quaternary ammonium compounds (QACs), benzalkonium chloride (BAC) 12, and 14 were also prevalent in the influent, especially in WWTP1, with values up to 8.22 μg L −1 . These QACs are used as surfactants in cleaning products and disinfectants. Qac resistance genes are frequently associated with class 1 integrons, and other ARGs included in those MGEs ( Gillings et al., 2009 ). QACs selective pressure might entail the co-selection of MGEs and their associated ARGs. Possibly, the use of disinfectants by a neighbouring dairy industry in the catchment area of WWTP1 (contributing to ~12% of the influent) could explain these higher levels.
The occurrence of a selected panel of ARGs and MGEs in the influent of the waterline followed similar trends across the three WWTPs (supplementary information Table S2). From the ARGs selected, ermB (6.39 log gene copies mL −1 ) and sul1 (5.85) had the highest average annual concentration values, while the βlactamase bla CTX-M had the lowest (4.05). The overall ARG patterns are in accordance with our previous study of more than 60 Dutch WWTPs ( Pallares-Vega et al., 2019 ) and in other recent works ( Di Cesare et al., 2016 ;Rodriguez-Mozaz et al., 2015 ). The high concentration of ermB gene in the influent cannot be associated with a direct antibiotic selective pressure within sewage, as erythromycin residues were not detected in the influent. The high occurrence of ermB might result from its location in Lactobacillales, which are common in the gut microbiome and, therefore, predominant taxa in the influent ( Ali et al., 2019 ;Cai et al., 2014 ). The high occurrence of sul1 may be explained by its extended use in the past and its association with MGE, such as integron class 1. Their presence could also be maintained by the persistence of sulfonamide antibiotic residues in wastewater ( Baran et al., 2011 ). As for the prevalence of MGEs, both korB (standing for incP-1 plasmids) and intI1 , encoding for the integrase of class 1 integron, were ubiquitous in the influent samples. Moreover, intI1 had a significantly higher concentration (p < 0.01) in the influent of WWTP1 (7.04 log gene copies mL −1 ) when compared to the other two WWTPs (6.16 logs on average), and above the range measured in our previous study ( Pallares-Vega et al., 2019 ). These results might be explained through co-selection events by the extended use of QACs within the dairy industry facilities. Further analysis addressing the presence of qac resistance genes and qac-intII1 relation would be necessary to confirm such a hypothesis.
To investigate the role of rainfall on the occurrence of antimicrobials, ARGs, and E. coli in influent, both their concentrations and their daily loads per population equivalent (pe) ( i.e. the absolute number of gene copies passing the WWTP per day divided by the population equivalent) were studied. Unlike concentrations, the daily load per pe should be constant despite differences in rainfall dilution if freshly discharged human feces were the only source of these compounds. Therefore, using daily loads or daily loads per pe (of both genes and E. coli ) for graphical representation and model response should be better suited to detect possible effects of seasonal temperature.
Increased rainfall led, as expected, to decreased concentrations of E. coli (-0.25 logs per average daily flow, p < 0.001) but not to decreased daily loads per pe ( + 0.06 logs per average daily flow, p = 0.15), Fig. 2 , supplementary information Figure S3 and Table  S3 models 5-7, confirming the dilution effect of E. coli upon heavier rainfall. In contrast, for the studied genes, the reduction in concentrations with increased rainfall was less clear (the best model did not include HLF as a determinant), and there was a significant positive effect of increased rainfall on the genes when the daily loads per pe were used as the response variable ( + 0.42 logs per average daily flow, p < 0.001), Fig, 2 , supplementary information Figure S3 and Table S3 models 1-4.
An increase of the daily loads per pe of ARGs and MGEs with increased rainfall might point to an additional source of genes besides freshly discharged feces. We hypothesized that such a source could consist of resident antibiotic resistant microbiota in the sewers, located for instance, within sewer biofilms or sewer sediments ( Auguet et al., 2017 ). These might also have been originally introduced into sewers with fecal microbiota. With increasing flow due to rainfall, a washout of the sewer microbiota could increase the incoming loads of resistance genes per pe, similar to washout of in-sewer stocks of, e.g. organic matter ( Gromaire et al., 2001 ). The contrasting behavior of the loads per pe of E. coli during rainfall events might indicate a minor accumulation of this organism in the sewer pipes. A limited accumulation of E. coli O157:H7 and gammaproteobacteria in sewer biofilm has been previously detected ( Auguet et al., 2017 ). However, the observed discrepancies between genes and E. coli might also be a consequence of a methodology bias. Unlike qPCR, culturable-based methods account for neither the dead nor a viable-but not-culturable fraction of bacteria.
With respect to seasonal temperature, the studied agents in influent (antimicrobials, genes, and E. coli ) showed an inconsistent response ( Fig. 2 , Table S3 models 1-7). The antibiotic loads varied per plant, with only WWTP2 showing an increase of antibiotic loads at colder temperatures. With respect to resistance genes, a slight but significant decrease of incoming genes (for both concentration and daily loads per pe) with increasing temperatures was observed for the set of genes as a whole (-0.02 logs per °C, p < 0.05). In contrast, increasing temperatures significantly enhanced both concentrations and daily loads per pe of E. coli in the influent of all three WWTP. The effect was mild as per degree increased (0.03 logs increase per °C, p < 0.001).

supplementary information Figures S2-S4 and
Overall, there is limited information regarding seasonal fluctuations in influent of both antimicrobials and ARGs. Some studies have observed an increase in loads of some antibiotics during winter months ( Coutu et al., 2013 ;Marx et al., 2015 ). This has been related to higher consumption of antibiotics used to treat winterrelated conditions ( Marx et al., 2015 ), which could increase the selective pressure and enhance the occurrence of resistant bacteria and ARGs in the sewage. This hypothesis is supported by three studies conducted in German and Chinese WWTPs. In these studies, either higher relative abundance of several ARGs ( Caucci et al., 2016 ;Jiao et al., 2018 ) or higher absolute concentration of numerous bla genes ( Schages et al., 2020 ) were found in cold seasons. In contrast, Karkman et al. (2016) did not observe any difference across seasons after quantifying the relative concentration of a broad set of ARGs in a Finnish WWTP. Differences in antibiotic prescriptions across countries and the degree of seasonality might explain the contradicting observations. Moreover, integrating flow variations of the sampling days into the studies might contribute to explain the uneven results.

Removal of E. coli and gene determinants through conventional water treatment and novel aerobic granular sludge
All three WWTPs significantly (p < 0.001) removed the fecal indicator bacteria and the tested genes. The removal efficiency varied across WWTPs and measured agent ( Fig. 3 , supplementary information Table S3 ). WWTP1 achieved the best removal for both E. coli and genes, with an average removal of 2.31 logs for E. coli ( supplementary information Figure S5 ) and ≥2 logs of tested genes (except for korB ). The other two WWTPs performed significantly worse ( + 0.4 logs p < 0.001) removing both ARGs and E. coli ( + 0.2-0.4 logs p < 0.05, supplementary information Table S3 , models 11 and 12), although within the range of removal previously observed for Dutch WWTPs ( Pallares-Vega et al., 2019 ). Moreover, the patterns for gene removals were, in general, similar in all three WWTPs and are following our previous study. The most successfully removed genes were ermB (2-3 logs on average), tetM , and bla CTX-M (2 logs on average). Bla CTX-M was undetectable or unquantifiable in 10-40% of the effluent samples of WWTPs 3 and 1, respectively. The most resilient ARGs of the panel were those relating to sulfonamides ( sul1 and sul2 ) with average removals that ranged from 1 to 2 logs. The two MGEs genes, intI1 and korB , were also more resilient to the treatment, with removals in the range of 0.6-1.5 logs on average, except for WWTP1 in which intI1 was significantly better removed (2.75 logs, p > 0.01). Therefore, although WWTP1 received a more considerable amount of intI1 gene, it succeeded in removing it to the same or lower levels than the other two WWTPs.
For the greatest part of the measured genes, the wastewater treatment did not exacerbate but rather decreased the relative abundance of the studied ARGs. For some of the genes ( intI1, sul1, or sul2) , a non-significant change or a slight relative increase was found in some of the WWTPs. In contrast, korB relative abundance increased significantly (p < 0.001) after the treatment in all three WWTPs ( supplementary information Figure S6 ). These data confirm our previous observations for one-time measurements in 60 WWTPs for more extended sampling periods.
A sampling of intermediate steps within the treatment was performed to evaluate the contribution of each treatment step to the removal of both E. coli and genes. The primary treatment step (WWTP2) did not affect the removal of genes and exhibited a moderate but significant effect in removing the fecal indicator (-0.11 logs, p < 0.05). The A-stage (AB line, WWTP3) moderately removed E. coli (-0.17 logs, p < 0.05) and genes (-0.36 logs, p < 0.001). Therefore, the greatest removal of both E. coli and the genes occurred in the biological nutrient removal stages. Activated sludge with short solid retention times and short clarification, as in the Astage, is likely not sufficient for the removal of pathogens or ARGs.
Lastly, we compared the removal efficiencies of two parallel lines-AB-line (based on flocculent sludge) and aerobic granular sludge-treating the same influent. Aerobic granular sludge is a modern water treatment technology requiring smaller space and footprint than conventional activated sludge systems ( Pronk et al., 2015 ). Aerobic granular sludge is based on bacterial aggregation in granules instead of flocs. This configuration comprises a different spatial distribution and bacterial community that could affect the removal of ARGs compared to flocculent sludge. The presence of ARG in granules has so far only been studied concerning accumulation during the granulation processes in bench-scale aerobic granular sludge reactors ( Li et al., 2020 ). However, information on the occurrence of ARGs in the sludge fraction of full-scale AGS installations as compared to conventional sludge is missing.
After a one-year of sampling, no significant differences were observed in removing ARGs and MGEs among the two parallel treatments. Exceptionally, ermB gene was better removed in the AB system than in the aerobic granular sludge system ( + 0.22 logs, p < 0.01). The removal of E. coli was also similar to the AB system ( Figure S5), in line with a recent study addressing the removal of fecal indicators ( Barrios-Hernández et al., 2020 ).

The removal of genes and E. coli is compromised by high hydraulic loads and effluent suspended solids but not by seasonal temperature
The effect of abiotic parameters (HLF, turbidity, and average temperature) on the removal of ARGs, MGEs, and E. coli , was investigated through linear mixed models ( supplementary information Table S3 models 8-10). Irrespective of the type of wastewater treatment, the removal of both E. coli and genes was hampered at high HLF ( Fig. 4 ). The removal capacity was modeled to decrease by 0.53 log CFUs (p < 0.001) and 0.35 log gene copies (p < 0.01) at double the average daily flow. This gene removal rate is in good agreement with that obtained in our previous study (-0.38 logs) based on single measurements across many plants. Higher turbidity in the effluent was also correlated with a minor but significant decrease in the removal of E. coli , (-0.01 logs per 1 Nephelometric turbidity unit −1 p < 0.05) and genes (-0.02 logs p < 0.05), Fig. 4 . In contrast, seasonal changes in the average air temperature on the day of sampling did not alter the removal of E. coli nor genes ( supplementary information Figure S7 ). Hence, opposite to what was observed for influent, variation in flow was the leading cause of variability, and the seasonal temperature had no contribution. The mechanisms by which the removal capacity of WWTPs might be disturbed with the increasing flow have been discussed previously ( Pallares-Vega et al., 2019 ). In short, increasing flow causes wastewater to spend a shorter time in the biological treatment and sedimentation steps.
The lack of effect of seasonal temperature in the removal capacity might appear unexpected, since fluctuations in seasonal tem- perature are known to significantly shape the bacterial community composition within the activated sludge ( Griffin and Wells, 2017 ) and alter the treatment performance, i.e., by nitrification failure during winter ( Johnston et al., 2019 ). Surprisingly, the impact of seasonal temperature on ARGs and E. coli removal during full-scale wastewater treatment is seldom reported and therefore remains poorly understood. From the available studies, no statistical differences can be found regarding the seasonal occurrence of E. coli in full-scale WWTPs effluents ( Lépesová et al., 2019 ;Osi ńska et al., 2020 ). Moreover, Barrios-Hernandez et al. (2020) described no effect of seasonal temperature on the removal of E. coli .
Seasonal peaks of absolute ARGs in effluent have been previously reported in winter and spring  or summer ( Jiao et al., 2018 ). In this study, seasonal fluctuations in absolute effluent concentrations (moderately higher with lower temperatures) were observed only in WWTP2 ( supplementary information Figure S8). Despite this trend, our results indicate that changes in the seasonal temperature did not influence the removal rates in any of the WWTPs ( supplementary information Figure  S7 ). In contrast, Jiao et al. (2018) reported a better removal of ARGs during summer. However, in Jiao ś study, the effect of temperature cannot be detached from that of flow dynamics (highly significant according to our results), because information about the flow was not included. The degree in which temperature or flow influence the removal efficiency of the treatment might vary across countries with different tem perature and precipitation regimes. Thus, additional studies in other locations accounting for both flow and temperature might be needed to understand further the role of tem-perature and flow in the removal of resistance determinants and E. coli .

Fate of antimicrobials and disinfectants during wastewater treatment
The fate of the different antimicrobials and disinfectant residues during wastewater treatment depended on the compounds studied ( Fig. 5 and S8 ). Some compounds were found both in effluent and biosolids (azithromycin, ciprofloxacin, sulfapyridine), while others were present either in the effluent samples (sulfamethoxazole, trimethoprim, clarithromycin) or in the biosolids line as the tetracyclines (tetracycline and doxycycline) and the disinfectants (BAC 12 and BAC 14 ). Last, although erythromycin was not detected in any of the influent samples, it was sometimes present in AGS and DS from WWTP3. All types of treatments, including that based on granular sludge ( Figure S9 ), reduced to a similar extent the antimicrobial concentrations (2-10-fold, depending on the compound). Specifics of the concentrations in each WWTP ś effluents are gathered in supplementary information Table S1 . Despite the partial decrease, eight of the tested antimicrobials were still detectable in some of the effluent samples, although only ciprofloxacin and azithromycin were above the PNEC levels. Most of those compounds have not been commonly detected either in the upstream or downstream surface waters of the WWTPs discharge points ( Sabri et al., 2018 ;Verburg et al., 2019 ). Hence, despite WWTPs discharging antimicrobials into the receiving waterbodies, the residues are diluted and/or sorbed to sediment, reducing their concentrations in the water stream below the limits of detection. The compounds that sorbed to the biosolids line were found in both DS and AS, although with higher concentrations in the DS samples than in the AS samples, likely derived from the difference in the solids content of each type of sample. Ciprofloxacin was again the most common antibiotic residue (2-4 mg kg −1 TS). Ofloxacin, while barely detected in influent, was often detected in the biosolids line but in lower quantities than ciprofloxacin, in agreement with previous studies in Europe ( Lindberg et al., 2005 ;Radjenovi ć et al., 2009 ). In general, the concentration of tetracyclines and sulfapyridine followed the trends observed elsewhere ( Göbel et al., 2005 ;Lindberg et al., 2005 ;Shafrir and Avisar, 2012 ). Tetracyclines were often detected in values ranging between 0.1-1.2 mg kg −1 TS for both AS and DS samples, which is around 5 to 40 μg kg −1 of fresh digested sludge for tetracycline and doxycycline respectively. Concentrations of 15 μg L −1 of tetracycline (150 times below the minimum inhibitory concentrations) have shown to enhance the growth of tet resistant bacteria ( Gullberg et al., 2014 ) and to stimulate horizontal gene transfer events in vitro ( Jutkina et al., 2016 ), although the bioavailability of these residues in the biosolids may reduce such an effect.
Besides the aforementioned antimicrobial residues, the two disinfectants tested were also highly sorbed onto sludge. Concentrations ranged between 1-14 mg kg −1 TS in CAS-like AS samples and 3-23 mg kg −1 TS in DS samples. The highest levels were reported for the AS in the A stage of WWTP3 with up to 49 and 101 mg kg −1 TS for BAC 12 and BAC 14, respectively ( supplementary information Figure S9 and Table S4 ). These concentrations meet with literature reports ( Martínez-Carballo et al., 2007 ) and reflect the important accumulation of these compounds in the sludge. A high occurrence of BACs has been shown to hamper methanogenesis in anaerobic digesters ( Zhang et al., 2015 ). Moreover, field amendments of BACs rich biosolids could result in the accumula- tion of these compounds, especially in clay soils, which could potentially lead to the selection of qac genes and co-selection of ARGs ( Mulder et al., 2018 ).

The occurrence of resistance determinants in biosolids mirrored those in the influent
The occurrence of the ARGs in the biosolids line of the three WWTPs reflected the patterns observed in the influent ( supplementary information, Figure S10 ). The highest concentrations were found for ermB, and sul1 with 9.4-9.9 log copies g −1 TS, respectively ( Fig. 6 ). On the lower rank were once again qnrS and bla CTXM , which often laid below the limit of detection or quantification (up to 86% of some of the DS and AGS samples). Contrarily, a recent study analyzing a broad range of ARGs suggested no contribution of influent ARGs to the recycled activated sludge resistome, which was richer in abundance but poorer in diversity when compared to the influent ( Quintela-Baluja et al., 2019 ). Since mixed liquors and not sedimented sludge were used in our study, a higher resemblance to influent ARGs patterns can be expected. Our data also demonstrate a high occurrence of MGE elements in the biosolid line, particularly in the AS systems, where the korB gene was found in similar ranges than the integrase ( intI1) gene. In contrast, intI1 was 1-2 log more abundant than korB in influent samples . The high prevalence of korB in the activated sludge might explain the poor removal of this gene and the subsequent equalization of intI1 and korB levels in the effluent. IncP-1 plasmids have been detected in biosolids of activated and digested sludge ( Dröge et al.,20 0 0 ). However, to the best of our knowledge, this is the first time incP-1 plasmids have been quantified in activated sludge samples showing a high occurrence, which confirm their relevance in studies addressing horizontal gene transfer events in biosolids-like systems.
Neither the flow nor the seasonal temperature seemed to homogenously alter the concentration of genes in the CAS-like activated sludge systems ( supplementary information Table S3, Fig-Fig. 6. Absolute abundance of 16S rRNA, ARGs, and MGEs in biosolids line of three Dutch WWTPs through a year. Abbreviations: AS: Activated sludge; AS A stage: Activated sludge from the A stage of the AB system, AS B stage: Activated sludge from the B stage of the AB system; AGS: Aerobic Granular Sludge treatment; DS: Digested Sludge ure S11 ); Even though some trends could be observed across plants with increasing temperatures, the effect varied per studied gene (i.e. a decrease of ermB and tetM and an increase of qnrS absolute concentrations). Overall, there was a modest variation of concentrations of ARG in biosolids between WWTPs (roughly 0.5 logs) and across the year. However, the absolute abundance of most of the genes was greater in AS than DS. The same effect was observed for the E. coli concentrations ( supplementary information Figure S5, Table S5 ). Lower absolute concentrations of ARGs, MGEs, and E. coli were also observed among AGS in contrast to the AS from the B stage (which is comparable to a conventional activated sludge system). These differences are likely due to normalization per gram of TS, which is roughly 10-fold higher in the DS and AGS samples compared to AS ones. When normalized to the 16S rRNA ( supplementary information, Figure S12 ), the relative concentration of several of the ARGs and MGEs in biosolids was similar among the aforementioned pairs (AS vs. DS and AGS vs. AS from the B stage). Only a slightly lower relative concentration of ermB was observed in AGS. A lower adhesion of bacteria harboring ermB gene to the granular sludge fraction could explain its poorer removal after AGS treatment in comparison with the AB line.
Consistently higher relative concentrations of ARGs in AS A stage were observed in comparison to AS from B stage, and another flocculent AS from WWTP1 and WWTP2, and AGS ( Figure  S12 ). This difference could be due to the operational conditions of A stage, where the solids retention time is significantly shorter (0.3 days) than that used for B stage (23 days), and another con- A persistent higher concentration of tetM gene after the anaerobic digestion was observed when compared with AS samples for both absolute and relative abundances ( Figure S12 ). This suggests that the anaerobic treatment might select for bacteria harboring this gene. As aforementioned, the concentrations of tetracycline residues in digested sludge might also contribute to the selection of tet genes, although similar effects were not found for quinolone resistance. A slight enrichment of the relative abundance after anaerobic digestion was also observed for ermB ( Figure S9 ). These findings are in accordance with the results of Ma et al. (2011) in bench-scale mesophilic digesters, where they even observed an increase in the absolute abundance of erm genes and some of the tested tet genes. In contrast, intI1 and sul genes decreased in both relative and absolute abundance. An increase of several ARGs, including erm, tet, and sul genes was also observed in two full-scale anaerobic digestors in China ( Tong et al., 2019 ), while in another full-scale study in the US, the relative abundance of three tet genes varied depending on the sampling dates ( Ghosh et al., 2009 ). Digested sludge is used in some countries as fertilizer for crops. The impact of pathogens and ARGs from sludge amendments in soil is still under debate ( Rahube et al., 2014 ;Rutgersson et al., 2020 ). In The Netherlands, digested sludge undergoes incineration. However, there is a growing interest in nutrient recovery from this byproduct; therefore, increasing the knowledge of possible hazards in the handling and downstream processing of digested sludge is important.

Depicting sampling strategies
If the results obtained in this studied are compared with those of our previous work (sampling multiple WWTPs but on a single occasion), the former managed to capture similar variability in ARGs occurrence and removal as in repeated sampling across one year ( supplementary information Figure S13 ). Thus, shorter sampling efforts might be enough to evaluate the removal abilities of a WWTP. If the objective is to evaluate variability in performance (and address possible solutions), rainy and dry periods might be more interesting to assess than seasons.
In the absence of molecular methodology or the need for rapid results, the use of bacterial surrogates to evaluate the fate of ARGs might be necessary. Correlation analysis (Pearson's correlations) highlighted that E. coli should not be used to evaluate the variation in incoming concentrations of ARGs but could be considered as a surrogate to evaluate the removal of specific ARGs such as bla CTX-M , ermB, and tetM ( supplementary information Figures  S14-S16 ), commonly associated to Enterobacteriaceae and Lactobacillae. These taxa follow E. coli removal patterns during wastewater treatment ( Barrios-Hernández et al., 2020 ;Ferreira Da Silva et al., 2007 ;Ottoson et al., 2006 ).

Conclusions
A one-year sampling campaign of three full-scale WWTPs highlighted that warmer seasonal temperature marginally decreased the concentrations of resistance genes in the influent but increased those of E. coli. However, seasonal temperature variation had an impaired effect on concentrations of antimicrobials in the influent. Instead, rainfall played a major role by diluting the concentrations of antimicrobials as well as fecal indicators such as E. coli , but not of resistance genes. Rainfall increasing the typical hydraulic load of each WWTPs significantly reduced the efficiency of wastewater treatment removal of genes and E. coli, in agreement with previous findings across The Netherlands . Increasing effluent ś turbidity was also related to slightly poorer removal. In addition, we concluded that the occurrence of resistant determinants in the biosolids line followed the occurrence patterns in the influent and that incP-1 plasmids are highly abundant in biosolids. Finally, full-scale activated sludge and granular sludge technologies displayed comparable performance in the ability to remove antimicrobials, resistant determinants, and the fecal indicator E. coli.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. tireless effort during the sample collection and to Gonçalo Macedo for his assistance throughout the data analysis. We finally acknowledge Erwin Tuinhof and the analytical team at Wetsus for the analysis of the antibiotic residues.
This work was performed in the cooperation framework of Wetsus, European Centre of Excellence for Sustainable Water Technology (www.wetsus.eu). Wetsus is cofunded by the Dutch Ministry of Economic Affairs and Ministry of Infrastructure and Environment, the European Union Regional Development Fund, the Province of Fryslân, and the Northern Netherlands Provinces.
Besides, this research has received funding from the European Union's Horizon2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 665874 and was cofunded by STOWA . The authors would also like to thank the members of the research theme Source Separated Sanitation for the shared knowledge and financial support.

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.watres.2020.116752 .

Quality of the extractions
To assess the efficiency and quality of the DNA extraction, all the samples were spiked with an internal standard consisting of 1 x 10 7 copies of synthetic gene fragments (gBlocks, IDT technologies, IA EE.UU) of the synthetic blue fluorescence protein ( bfp ) gene prior to extraction. DNA extracts were quantified using a Quantus TM fluorometer (Promega, NL) according to the supplier's instructions. DNA quality was assessed by gel electrophoresis (agarose at 1.5% m/v) and by measuring absorbance ratio at 260/280 nm and 260/230 nm using a Nanodrop spectrophotometer (ThermoScientific, UK)

qPCR: oligonucleotides, probes, and reaction and conditions
Preparation of qPCR reagent mix and reaction conditions was performed as indicated in Pallares-Vega et al., 2019 , except for the korB assay, for which an increase of the primer concentration was used (400 nm), and for the inclusion of the bfp assay, that follows the average reaction conditions with annealing temperature at 60 °C Table A1 Oligonucleotides and probes used for gene detection by qPCR reactions. In primes/probes with degenerate code, Y stands for pyrimidine bases (C or T), R stands for purine (A or G), S for strong bases (C or G), and V for A, C, G (IUPAC nomenclature).

Sample preparation for antimicrobials and disinfectants measurements in biosolids
For the sludge line samples, aliquots stored at -20 °C were defrosted overnight at 5 °C. AGS samples were homogenized by beadbeating at 4500 rpm for 30 s in the Precellys homogenizer (Bertin Technologies SAS, FR) with the help of 4 mm glass beads (Merk, NL). AS and DS samples were homogenized by vigorous manual agitation. A total of 0.15 g of homogenized AGS and 0.75 g of homogenized AS or DS were used for the analysis.
The conditioning of the sludge matrix was achieved by mixing the sample with 1.5 mL of buffer (ammonium formate/formic acid (50:50, v/v), pH 2). To this mix, the following was added: 0.1 mL of an internal standard with isotopically labeled compounds ( Table B3 ), 1.5 mL of modifier (consisting of in 100 mL: 50 mL of oxalic acid at 1 mol L -1 , 15 mL of ammonia at 5 mol L -1 , 5 mL of formic acid at 99% (v/v) and 35 mL of ultrapure deionized water), 3 mL of methanol at 99% v/v and 1 mL of organic modifier (acetonitrile/methanol (50:50 v/v) + 1 % formic acid). In order to calculate the recovery of each compound in each sample, a parallel vial was prepared, including 0.4 mL of a standard containing a mix of all the tested compounds. For both samples and recovery surrogates, the volume was completed up to 15 mL with ultrapure deionized water. The mix was vortexed for 30 min at speed 8 (1700 rpm) and subjected to sonication for 15 min in a bath sonicator (Bandelin electronic, DE). The vials were centrifuged at 3475 x g for 10 min, and the supernatants used directly for analysis by LC-MS/MS. All samples were injected in an Agilent 6420 Triple Quadrupole LC-MS/MS system with an electrospray ion source. All the compounds were detected in the positive mode after separation in a ZORBAX Eclipse Plus C18 RRHD L = 50 x d = 2.1 mm column with 1.8 μm particle size. Detailed information about the mobile phases and data analyses can be found in Tables B1 , B2 , B3 , B4 .
The limit of detection (LOD) and limit of quantification (LOQ) of the method were determined for each compound as the lowest detectable amount of compound with a signal-to-noise ratio of 3 and 10, respectively. The recovery rates for each compound in each Table B1 Composition of mobile phases and parameters for liquid chromatography (LC) separation of the antimicrobials and disinfectant residues.

Mobile phase
Mobile Phase A Positive electrospray ionization: 2,5 L ultrapure deionized water + 5 mL formic acid (99% v/v), 0,5 mol L -1 ammonia 5 mol L -1 + 0,1 mL Oxalic acid 1 M. Negative electrospray ionization: 2,5 L ultrapure deionized water + 5 mL ammonia (5 mol L -1 ) + 1 mL Formic acid (99% v/v) + 0,1 mL Oxalic acid 1 mol L -1 . Mobile Phase B Positive electrospray ionization: Acetonitrile + 0,1% Formic acid. Negative electrospray ionization: Acetonitrile. sample were calculated from the spiked samples. The samples values were then recalculated by multiplying each result by the corresponding recovery, only if this was among 50-150%. When the recovery value was below or above the aforementioned threshold, the sample was excluded. For the disinfectant residues BAC 12 and BAC 14 , the recoveries were not applied, as the concentrations already present in the samples were 10-30 times higher than the spiked concentrations. Thus, the results displayed are the concentration in the sample without further recalculations.
Only for the statistical analysis, the final effluent samples having bla CTX-M values below the LOQ were replaced by the LOQ value (6 out of 36 values). All other genes were above the LOQ in all samples before the removal was calculated.
For the comparison analysis of single gene occurrences in influent, or the comparison of single genes and E. coli removal performance across all three WWTPs or across the two parallel lines of WWTP3 and its final effluent, an analysis of the variance was used (when normality was met) followed by Tukey post-hoc analysis. The comparisons were made using the log concentrations of genes mL -1 or the log transformed removal values, respectively. If the distribution did not meet the normality, the group comparison was performed with non-parametric test (Kruskal-Wallis).
The influence of abiotic factors in the incoming and removal of antimicrobials genes and E. coli as well as the role of intermediate steps in the removal of genes, were analyzed by linear models and linear mixed models. The summary of the models is displayed in Table C1 , and the construction of the models is described below: To evaluate the contribution of either the overall treatment or the intermediate steps (primary treatment in WWTP2 or A-stage in the AB line of WWTP3) to the removal of genes determinants and E. coli , linear mixed models were used with observations clustered by sampling time-point. The log-transformed concentrations (log10 of gene copies of ARGs, MGEs, and E. coli log10 CFU counts per mL) from each sample type were used as the response variable. The explanatory variables tested in the mixed model were "sample.type", fixed term, representing the location of the sample within the WWTP, and the "gene.type" (only in the gene models), "WWTP" and "sample.code.month" as independent random terms (random intercept modelled). The latter allowed the model to account for paired measurements from the same month of influent (IN) and final effluent (FE), primary treatment (PT) or A stage (AST). The result of these models are coefficients describing the gene reduction per location, and differences in gene concentrations per WWTP, across all resistance genes or E. coli and sampling timepoint.
For single plants, similar linear mixed models (or linear models as for E. coli ) were also used to investigate the influence of additional explanatory factors (sampling parameters) on the removal of either ARGs and MGEs or E. coli CFU counts. In this case, the response variable was the removal value, calculated as the log10 of the ratio of the concentration of genes or CFU counts in the influent versus the final effluent. The explanatory variables (fixed terms) were the WWTP, the average temperature, the turbidity (as a surrogate for TSS presence in effluent), and the hydraulic load factor (HLF). The random terms were the "gene type" (only in the gene model) and the "sample.code.month" that allowed grouping all the genes from the same WWTP and sampling time point. For E. coli , only one explained variable (CFU counts) was available, and thus neither "gene type" nor "sample.code.month" random terms were applicable. The inclusion of WWTP as fixed term lead to singular fit problem. Thus, for E.coli , a linear model was used instead.
The role of explanatory factors for the concentrations of resistance genes and E. coli in influent was also investigated through linear mixed models. The response variable was the log10transformed influent concentration per mL of either ARGs and MGEs gene copies or CFU counts. In addition, the "load.pe" or absolute daily amount of resistance genes and CFU per population equivalent was used as a response variable. The load was obtained from multiplying influent concentrations with the flow on the measurement day, thereby correcting for increased treatment volumes. Last, the load was normalized per population equivalent (pe) "load.pe", (where 1 pe stands for 150 g of total oxygen demand), and used as a response variable (log10 transformed). The average atmospheric temperature and HLF from the day of sampling were used as the fixed term explanatory variables. Again, "gene type" (only in the gene models) and "sample.code.month" were used as the random terms whilst HLF and temperature were used as fixed effects.
The same explanatory factors were also modelled against the concentration of genes in the biosolids, by using a similar approach than for the influent models, but with the log10 transformed concentration of genes (ARGs and MGEs) or CFUs per g of TS as a response variable. Due to high percentage of missing values, qnrS and bla CTX-M genes were not included in the model. Finally, a linear mixed model was used to address the effect of both temperature and HLF in the incoming concentrations of antimicrobials commonly found in influent (azithromycin, clarithromycin, ciprofloxacin, sulfamethoxazole, sulpyridine, trimethoprim, doxycycline, tetracycline, and clindamycin). Missing values were replaced by LOD/ √ 2 of each compound. The response variable was the concentration of antimicrobial in μg L -1, and the explanatory variables were the HLF and temperature in the fixed terms, and "antimicrobial.type" to group for each type of antibiotic and "sample.code.month" in the random terms to account for samples from the same location and sampling time-point. After analysis, the model was considered to not being well-fitted, according to the optical inspection of residuals, and thus, results of this model are not further included in the results and discussion section.
In the linear mixed models with several explanatory variables, the relevance of the explanatory variables was determined through stepwise backward model reduction, and the quality of the models was assessed by visual inspection of the normality of the residuals.
The datasets used in this study are available in Mendeley Data repository doi:10.17632/53fk4cht32.1