Determinants of presence and removal of antibiotic resistance genes during WWTP treatment: A cross-sectional study

Wastewater treatment plants (WWTPs), linking human fecal residues and the environment, are considered as hotspots for the spread of antimicrobial resistance (AMR). In order to evaluate the role of WWTPs and underlying operational parameters for the removal of AMR, the presence and removal ef- ﬁ ciency of a selected set of 6 antimicrobial resistance genes (ARGs) and 2 mobile genetic elements (MGEs) was evaluated by means of qPCR in in ﬂ uent and ef ﬂ uent samples from 62 Dutch WWTPs. The role of possible factors impacting the concentrations of ARGs and MGEs in the in ﬂ uent and their removal was identi ﬁ ed through statistical analysis. ARGs and the class I integron-integrase gene ( intI1 ) were, on average, removed to a similar extent (1.76 log reduction) or better ( þ 0.30 e 1.90 logs) than the total bacteria (measured as 16S rRNA gene). In contrast, broad-host-range plasmids (IncP-1) had a signi ﬁ cantly increased (p < 0.001) relative abundance after treatment. The presence of healthcare institutions in the area served did only slightly increase the concentrations of ARGs or MGEs in in ﬂ uent. From the extended panel of operational parameters, rainfall, increasing the hydraulic load of the plant, most signi ﬁ cantly (p < 0.05) affected the treatment ef ﬁ ciency by decreasing it on average (cid:2) 0.38 logs per time the ﬂ ow exceeded the average daily ﬂ ow. Our results suggest that overall, WWTP treatments do not favor the proliferation of the assessed resistance genes but might increase the relative abundance of broad-host- range plasmids of the IncP-1 type.


Introduction
Antimicrobial resistance (AMR) is a growing problem worldwide. Although it is a natural and ancient phenomenon (D'Costa et al., 2011), its occurrence in natural environments has been accelerated by anthropogenic activities. One of the essential vectors for the dissemination of human-related AMRs into the environment is wastewater, as it collects fecal residues including antibioticresistant bacteria (ARB) and their genes (ARGs) (Baquero et al., 2008).
In Europe, wastewater is treated in wastewater treatment plants (WWTPs), and their effluents are commonly discharged into natural water bodies. The main goal of sewage treatment is to remove organic components (measured as chemical (COD) and biological (BOD) oxygen demand), phosphorus and nitrogen nutrients (P, N) as well as suspended solids, but not bacteria, or their genes (Council of the European Communities, 1991). The core biological secondary treatment units of WWTPs, involving activated sludge, are composed of open microbiomes comprising complex networks of microbial populations (Weissbrodt et al., 2014). The microbial communities of WWTPs are considered as hotspots for horizontal gene transfer (HGT) and putative proliferation of AMR (Berendonk et al., 2015;Rizzo et al., 2013).
Therefore, in the last years, several studies in different countries have evaluated the fate of ARGs in full-scale WWTPs (Czekalski et al., 2012;Makowska et al., 2016;Rafraf et al., 2016;Rodriguez-Mozaz et al., 2014). Whereas most studies have detected lower absolute concentrations of ARGs after treatment, inconsistent results have been found on changes in the concentration of ARGs relative to 16S rRNA ("relative abundance") (Lee et al., 2017;Makowska et al., 2016;Rafraf et al., 2016). Many factors might be the cause of this disparity, including changes in community composition along with the treatment, the presence of AMR selective agents in the wastewater, as well as the sampling design. Moreover, up to date, studies have rarely investigated the relationship between the efficiency of ARG removal and the process design and operational conditions of WWTPs. One approach has been to evaluate possible relations between the presence or removal of ARGs and water quality parameters as temperature, total organic carbon (TOC), BOD, nutrients and TSS (Ben et al., 2017;Di Cesare et al., 2016;Laht et al., 2014;Novo et al., 2013). Even though some correlations were found (e.g., between temperature and sulfonamide ARB, or between TOC and ermB, tetA and qnrS ARGs), these results did not universally apply to all investigated ARGs and are based on a low number of investigated WWTPs. Regarding the operation of the plant, most of the recent studies focus on assessing or comparing the efficiency of advanced treatments or disinfection technologies such as biological post filtration or disinfection by chlorination, UV or peracetic acid (Di Cesare et al., 2016;Laht et al., 2014).
In order to limit the emissions of ARGs to surface water, information on the efficiency of WWTP treatment and the role of plant processes is needed. Therefore, the aim of our study was 1) to evaluate the presence and removal of ARGs and MGEs in Dutch WWTPs and investigate changes in the relative gene abundance, and to elucidate 2) the influence of catchment area factors on loads of ARGs and MGEs in raw wastewater as well as 3) the role of WWTP process configurations for their removal. A selected panel of six relevant ARGs (Berendonk et al., 2015) and two mobile genetic elements (MGEs) related with AMR, intl1 and korB (IncP-1 plasmids), were analyzed across an extended number of 62 treatment plants. Moreover, information regarding the presence of possible explanatory variables in the catchment area (i.e., healthcare institutions), the amount of treated water during the sampling day (hydraulic load), and the design configuration and operational parameters of the WWTPs were gathered. The correlation of process and/or catchment parameters with ARG abundance and removals was studied with linear mixed models.

Characteristics of the selected wastewater treatment plants and sampling procedure
Influent and effluent samples were collected in 62 Dutch WWTPs distributed across the country (Fig. S1 in the supplementary information). The selected plants comprised varied sizes (4,800e1,060,500 Population Equivalents (PE) measured as 150g TOD) and different process configurations (Table S1 in the supplementary information). No installation with advanced treatment was included. The sampling was performed at a single time point per plant within the period dating from April to November of 2016. Detailed information about the sampling points and sampling procedures can be found elsewhere (Schmitt, 2017). Briefly, 24-h flow-proportional composite samples were gathered from the influent and the effluent of the WWTPs. Samples were processed within 24 h upon collection, as follows: a total of 30 mL of influent and 250 mL of effluent were filtered through 0.45-mm estercellulose membranes (Merk Millipore, DE) and filters were frozen at À20 C until extraction.

DNA extraction and quantification of ARGs and MGEs by quantitative polymerase chain reaction (qPCR)
DNA filters were extracted with the DNeasy kit Power Water (Qiagen, NL) following the manufacturer's instructions. DNA quantification and purification were determined by fluorometry using Qubit ® (Thermofisher, US).
The diluted DNA was subjected to quantitative polymerase chain reaction (qPCR) analysis of the selected panel of genes. The 16S rRNA gene was selected as a proxy to quantify total bacteria. qPCR targets were chosen from a proposed target panel to track ARGs (Berendonk et al., 2015). The chosen ARGs confer resistance to the antibiotics with the highest consumption in The Netherlands, namely: macrolides (ermB), tetracyclines (tetM), sulfonamides (sul1 and sul2), fluoroquinolones (qnrS) and extended-spectrum betalactamases (bla ctxM ) (Table 1). Moreover, two genes assessing the presence of MGEs were included in the study: intI1 and korB. The first target allows detecting integrase of class I Integron, a wellknown MGE related to the acquisition and exchange of ARGs  (Gillings et al., 2015). The second, korB gene, targets plasmids belonging to Incompatibility group P-1 (Jechalke et al., 2013) that are another kind of MGEs associated with ARG and serve as models to study HGT events in complex environmental samples (Klümper et al., 2015;Tsutsui et al., 2010). Further information about oligonucleotides, probes, reaction mix, and conditions can be found in the supplementary information and Table S2.

Catchment area and WWTP metadata for statistical analyses
Information about healthcare institutions was obtained from a separate project (Schmitt, 2017). In brief, information on localization of hospitals and polyclinics was obtained from a registry of Dutch health data maintained by the National Institute for Public Health and the Environment (www.volksgezondheidenzorg.info) and amended with a separate list of Dutch hospitals (https://nl. wikipedia.org/wiki/Lijst_van_Nederlandse_ziekenhuizen).
Information on care homes was obtained from a registry of Dutch health care institutions maintained by the Dutch patient federation. (https://www.zorgkaartnederland.nl/overzicht/sectoren/ verpleeghuizen-en-verzorgingshuizen/zorgaanbieders/plaatsen). Both were matched with the areas served by specific WWTPs. The Daily Flow of the WWTPs on the sampling days was obtained from the corresponding waterboards or WWTPs. Detailed information about plant design and performance parameters for the WWTPs in 2016 was obtained through the Dutch Statistical Office (Central Bureau Statistiek eCBS, 2018 https://www.cbs.nl/en-gb). The list of parameters taken into account for statistical analysis is summarized in Table 2.

Statistical analyses
Statistical analyses were performed with R 3.4.4 (R Core Team, 2018) and Rstudio (https://www.rstudio.com/). Correlation analysis (Pearson's correlation) between the removal of ARGs, MGEs, 16S rRNA, and the fecal indicator bacteria Escherichia coli were performed with the package corrplot (Wei and Simko, 2016). The data regarding E. coli removal was obtained from a previous study (Schmitt, 2017).
For the analysis of the relation of the absolute resistance gene concentrations in the WWTP's influent to catchment area factors, the log-transformed absolute concentration of 8 genes per liter of influent was set as the outcome variable. Bla CTX-M was excluded from the statistical analysis as the data set was not complete (this gene was found below the detection limit in >10% of the WWTPs effluent samples). Explanatory variables tested included the presence of hospitals, polyclinics and nursing homes in the catchment area as well as the effect of rainfall events increasing the hydraulic load (amount of water processed) of the WWTP during the sampling day. We named this variable "Hydraulic Load Factor," and it was calculated as the ratio of the Daily Flow (DF) during the sampling day over the average daily flow (derived from the annual flow) and expressed as "times x ADF." We used the annual flow as the dry weather flow (DWF) was not known for all plants. The hydraulic load factors observed in the 62 studied WWTP fall well within the distribution of hydraulic load factors retrieved from daily flows of two exemplary individual WWTP over two years and are therefore representative for Dutch conditions (data not shown). A linear mixed model involving the R packages lme4 and lmerTest (to determine p values through Satterthwaite approximation) (Bates et al., 2016;Kuznetsova et al., 2017) was used with the genes and the Hydraulic Load Factor and area parameters as fixed factors and a random intercept for the plant identity. Parameters were first tested univariably, including testing for interaction between gene identity and the other factors. Then, the MuMin package (Barton, 2018) was used for the identification of the best minimum adequate models by fitting all possible submodels using maximum likelihood methods, and ranking them by the corrected Akaike information criterion (AICc), retaining all models differing from the best model by less than an AICc of 2. Demonstration of single best models was chosen instead of model averaging (Dormann et al., 2018). The quality of the model was investigated by visually inspecting normality of the residuals.
For the analysis of the influence of plant parameters on the removal efficiency, the outcome variable was the log-transformed removal efficiency per gene per plant (i.e., the log reduction). Variables tested included the following: size of the plant (based on 150 g TOD PE ), presence/absence of primary sedimentation in combination with recirculation of dewatered sludge, type of P removal (none, biological, or biological and chemical), average sludge retention time (SRT), average hydraulic retention time (HRT), anaerobic contact time, and the average concentration of total suspended solids (TSS) in effluent. As denitrification was applied in all plants, it was not included in the statistical analysis. A linear mixed model was used, with the resistance gene and the plant parameters as fixed factors and a random intercept for the plant identity. In the first step, all parameters were tested in bivariate models (gene and plant parameter), and the significance of the interaction between genes and parameters was determined. To adjust for the effect of the Hydraulic Load Factor which was found to be highly significant, trivariate models were run (including the following three explanatory factors: gene, Hydraulic Load Factor, and their interaction, in addition to one additional parameter). Finally, all parameters (gene, Hydraulic Load Factor, and their interaction, the presence of primary settler with and without sludge recirculation, type of P removal, HRT, SRT, and TSS) were tested in a full model. The anaerobic contact time was excluded from this model since data on the design of this parameter was only available for 11 out of the 62 WWTPs. The MuMin package (Barton, 2018) was used for model reduction (i.e., identification of the best submodels of this model) using maximum likelihood methods, Table 1 Panel of genes used in this study. Genes are arranged within three groups of interest: all bacteria, antibiotic resistance genes (ARGs), and mobile genetic elements (MGEs).

Group
Gene Function All bacteria 16S rRNA For normalization to the concentration of bacteria ARGs ermB Resistance to macrolides sul1 Resistance to sulfonamides sul2 Resistance to sulfonamides tetM Resistance to tetracyclines qnrS Resistance to quinolones bla CTXM Resistance to extended spectrum b-lactams

MGEs intI1
Integrase of type 1 integrons korB Broad host range plasmids of incompatibility group incP-1 choosing the models differing from the best model by less than 2 AICc. The model reduction was performed on a subset of 37 plants with complete observations. The identified minimum adequate submodels were then re-run for the largest set of WWTP for which all relevant data was available. The quality of the model was investigated through analysis of normality of the residuals. Collinearity amongst explanatory variables was assessed by variance inflation factors.

Prevalence of ARGs and MGEs in the influent and role of the catchment area
We assessed the occurrence of different ARGs and MGEs in the influent of 62 WWTP. Our results, summarized in Fig. 1 and Table S3 in the supplementary information, showed that the most predominant genes in the influent were the ARGs sul1, ermB and the MGE intI1 (6.54 log copies mL À1 on average). The concentrations of tetM, sul2, qnrS, and MGE korB, were on average 10 times lower, while the lowest concentrations (4.40 log copies mL À1 ) were obtained for the beta-lactamase gene bla ctxM . These findings are in general accordance with the concentrations found in other European studies by Czekalski et al. (2012), Rodriguez-Mozaz et al. (2014) and Di Cesare et al. (2016), but 2e3 logs lower than the ones found in other northern European countries (Laht et al., 2014) ( Table S3). The resistance to sulfonamides (sul1, sul2) and macrolides (ermB) seem to be most wide-spread, although the use of these antibiotics in humans is not so extensive anymore (10% of the total antibiotic consumption in humans in 2016 (Nethmap/Maran, 2018)). Their high prevalence may be the result of the combination of prolonged use of these antibiotics in both humans and animal husbandry, their association with MGEs (Baran et al., 2011;Davies and Davies, 2010) and the presence of residues from these antibiotic families in wastewater (Schmitt et al., 2016).
The possible contribution of healthcare institutions to the ARGs and MGEs levels in influent of WWTPs was also investigated in this study. We found a small effect of the presence of hospitals on concentrations of resistance genes in influent: while the presence of hospitals was included in all best models indicating a possible role of hospitals as point sources of resistance, the increase in gene concentrations in influents due to presence of hospital wastewater was relatively small (<0.10 log unit), and hospital presence was not significant within these models. The effects of the presence of care homes and polyclinics were yet lower (Table 3 and Fig. S2 in the supplementary information). The role of healthcare institutions as an important source of antibiotic resistance emissions (bacteria and genes) into the sewer system has been already demonstrated in previous studies (Buelow et al., 2018;Lorenzo et al., 2018;Rodriguez-Mozaz et al., 2014). However, despite the elevated concentrations of ARGs in healthcare sewage systems as compared to community or industry wastewater, healthcare institutions' discharges are estimated to represent just 1% of the total influent Table 2 Studied variables and operational parameters that possibly affect ARGs and MGEs loads and removal efficiency. Information was gathered from the WWTPs, the waterboards, and the Dutch Statistical Office (CBS, 2018  volume (Kümmerer, 2009;Schmitt, 2017). Thus, the contribution of healthcare institutions wastewater was likely too small to increase the overall concentration of the genes tested in WWTP influents. Similar conclusions were recently drawn by Buelow et al. (2018).
Rainfall occurring on the sampling date, increasing the Hydraulic Load Factor, was found to decrease the concentration of ARGs in the influent slightly. This was likely a consequence of the dilution of human waste in the sewer with the rainfall inflow (Metcalf and Eddy, 2003). Still, the effect of rainfall on the gene prevalence (À0.11 logs/per time ADF was doubled) was not significant by itself (Table 3 and Fig. S3 in the supplementary information). These results could be explained by the homogenizing effect of the 24-h composite samples that would include wastewater from both rainshowers but also the dry period. Such an effect was also observed by Lucas et al. (2014) while monitoring fecal indicator bacteria in two WWTPs in Paris. Besides, the variability in both catchment area characteristics and the rainfall events could also have influenced the magnitude of the effect as was previously observed for fecal indicators and ARGs in influent and combine sewers overflows (Eramo et al., 2017;Lucas et al., 2014).

ARGs and MGE removal efficiencies
In the WWTPs studied, the average removal of ARGs was similar or higher than the average removal of total bacteria measured as 16S rRNA gene (Fig. 1, Table S3), meaning, that the average relative abundance of ARGs after treatment did not increase ( Fig. 2 and Table S4 in the supplementary information). We observed a 1.76 ± 0.40 log reduction for 16S, which implies a decrease of 98.2% on average. The highest removal rates were observed for qnrS, ermB and tetM genes (2.65 ± 0.68, 2.65 ± 0.74, 2.53 ± 0.68 average log reduction, respectively). For bla ctxM , 16% of the effluent samples had concentrations below the detection limit, and the average removal excluding those samples was 2.44 ± 0.56 logs. The sul2 removal was 2.00 ± 0.48 log copies, and closer to the 16S rRNA gene removal was sul1 with 1.82 ± 0.53 log copies reduction. Similar values were obtained for the MGE intI1 (1.80 ± 0.49). These results are comparable or better to those obtained in WWTPs including advance treatments or disinfection processes (Di Cesare et al., 2016;Laht et al., 2014;Wen et al., 2016) as described in Table S3. On the other hand, lower average reduction than the 16S rRNA gene was observed for the korB gene, with 0.89 ± 0.60 log removal (87.2%), meaning that the relative abundance of this gene was significantly increased (p < 0.001) after the treatment (Fig. 2, Table S4). Since broad-host-range plasmids as IncP-1 are known to disseminate into a great variety of environmental bacteria families, their removal might be countered by HGT events (Bellanger et al., 2014;Klümper et al., 2015). In addition, IncP-1 plasmids often include genes that encode for metal resistance and the degradation of xenobiotic compounds, thus, conferring metabolic advantages to bacteria in activated sludge and likely enhancing their dissemination (Dr€ oge et al., 2000). The presence of IncP-1 plasmids in WWTPs has been previously investigated by culturing and molecular-based techniques (Bahl et al., 2009;Moura et al., 2010), but to the best of our knowledge, this is the first time that their occurrence has been quantified, revealing an increase in their relative abundance after wastewater treatment.
Generally, it can be concluded that Dutch WWTPs do not contribute to enhancing antimicrobial resistance. The absolute concentration of the tested ARGs and of intI1 genes decreased on average 98.4e99.8% after treatment. Besides, the average relative abundance of ARGs either decreased or remained identical (although, a slight relative increase was found in some plants for sul1, sul2 and intI1, Table S1). Most of the available studies agree that WWTPs reduce the absolute numbers of both total bacteria and ARGs in wastewater. Yet, the effect of treatment on the relative ARG abundance differs greatly depending on the studied genes.  Table 3 Effect of operational parameters and catchment factors on concentrations of ARGs and MGEs in the influent. Estimates of the effects of the explanatory variables on gene concentrations (in log 10 copies/L) are given (with their p values between brackets and in italice significant estimates are shown in bold) for the 6 best models that were of nearly identical quality as determined by AICc. Gene identity was also included in all models. For interaction terms, the genes for which significant interactions with explanatory variables were found are listed (for korB the increase of concentrations in influent with hospitals was higher than for 16S; for qnrS the increase of concentrations in influents with care homes was lower than for 16S). Acronym: n: number of plants with available information for that parameter.  Significant differences observed after treatment on each gene were tested by a paired Wilcoxon test, and values are expressed above each gene (****): Highly significant, (p < 0.0001); ns: no significant differences observed. Rafraf et al. (2016) found an increase of up to 0.50e2.40 logs in the absolute concentration of some ARGs after treatment. The decrease of absolute and relative abundance of ARGs found in this study might be the result of the combination of the low human use of antibiotics in the Netherlands (possibly limiting selective pressures of these substances in sewage) together with continuous surveillance and upgrading of the Dutch wastewater facilities. Despite the average of 2.30 ± 0.30 log removal of resistance genes, Dutch WWTPs still release approximately 10 6 ARG copies per liter of effluent. The impact of the discharge of ARG-containing effluent on the receiving waterbodies was not evaluated in this study, but it has been addressed elsewhere (LaPara et al., 2015;Marti et al., 2013;Rodriguez-Mozaz et al., 2014;Sabri et al., 2018); in most of the cases, the discharge of WWTPs effluents increases the ARGs content in the receiving aquatic ecosystems. This illustrates that human exposure to ARG emitted from WWTP is possible, e.g., through recreation in surface waters. The exact public health burden of the presence of specific resistance genes in surface water is difficult to quantify, however, recreational exposure has been linked with higher ESBL carriage in surfers (Leonard et al., 2018).

Removal of ARGs as compared to MGEs, 16S rRNA gene, and E. coli
The removal of all ARGs and MGEs was positively correlated with the removal of the 16S rRNA gene (r ¼ 0.68e0.87), Fig. 3. Moreover, a strong and significant correlation was observed between the removal of ermB and tetM (r ¼ 0.96, p < 0.001). This could be explained by their typical co-location on diverse transposon families (Brenciani et al., 2007) that are usually present in Grampositive bacteria, mainly in the order Lactobacillales (Park et al., 2010). These bacteria are common fecal microorganisms present in the wastewater influent and in general, they are partially removed during WWTP treatment (Cai et al., 2014). The removal of another type of fecal bacteria, Escherichia coli, was also significantly correlated with the decrease of the beta-lactam gene bla CTX-M (r ¼ 0.79, p < 0.01) in accordance with the co-location of these genes in Enterobacteriaceae (Bradford, 2001). Sul1 and intI1 proved to be correlated in their persistence (r ¼ 0.92, p < 0.001), in line with the co-location of sul1 on class 1 integrons (intI1). Their resilience to treatment may be due to their association with diverse broad host range plasmids that are horizontally transferred among diverse bacteria (Gillings et al., 2015).

Influence of the design and process parameters of WWTPs in the efficiency of ARGs and MGEs removal
The impact of specific treatment processes on gene removal (as included in Table 2) was studied by statistical analysis (linear mixed models). These analyses were run both with and without the korB gene, as this was the only gene for which the relative abundance increased during wastewater treatment, suggesting a different effect of the treatment dynamics than for the rest of the gene panel.
Our results show that an increase of the Hydraulic Load Factor caused by rainfall events was the dominant variable explaining differences in reduction of ARGs and MGEs between the plants (significant effect in the univariate and most multivariate models, both with and without korB; Fig. 4, Table 4, Table 5, Tables S5 and S6 in the supplementary information). On average, the efficiency was reduced by 0.38 logs per time the ADF was exceeded. Therefore, rainfall reduced the incoming loads of ARGs in influent as well as the efficiency of their removal during treatment, yet only the latter was found to be statistically significant.
WWTPs are usually optimally operated in the so-called Dry Weather Flow (here represented by the ADF). While wastewater volumes up to a maximum flow of 3e6 times the DWF conditions can be processed (EPA, 1995), a higher hydraulic load can disturb the treatment processes. These disturbances are mainly related to the reduction of the WWTP's optimal HRT and differences in the influent composition (Metcalf and Eddy, 2003). In the secondary treatment, shorter HRTs and different influent composition (physical and biochemical) can affect both, biomass growth and the dynamics of biological processes, resulting in poorer treatment performance (Capodaglio, 2007). No significant effects of HRTs on removal rates were found in the univariate analysis, and only a slight effect of HRT was found in multivariate analyses (HRT was included in a few of the best adequate models, but not significant in itself). The increase in removal with 1 h HRT was limited in all these  A linear fit of the removal of resistance genes (16S rRNA gene, ermB, sul1, and qnrS, on log 10 scale) during WWTP treatment versus Hydraulic Load Factor measured in times the Average Daily Flow (ADF). models to 0.01 log gene copies. Yet, shorter HRT occurring during higher hydraulic loads on the sampling day likely contributed to poorer treatment efficiency.
When the increase in the hydraulic load is intense or lasts for a long period, it can cause washing out of activated sludge, leading to an increase in total suspended solids (TSS), including bacteria and ARGs, in the effluent (Capodaglio, 2007;Rouleau et al., 1997). TSS in effluent could also be increased by increasing TSS in the influents caused by rainfall events (Mines et al., 2007). In this study, effluent's TSS were only marginally influencing ARGs removal in multivariate analysis, and not significant in univariate analysis (Table 4). This was likely due to the use of annual average TSS values, therefore not reflecting the sampling day. However, TSS values from the sampling day were not significantly correlated with the concentrations of ARGs in effluents in a recent study (Ben et al., 2017). Those samplings were performed under dry weather conditions though. Therefore, the association between TSS and WWTP's removal efficiency during rainy weather events still remains unclear, and it is advisable to investigate the influence of this parameter in future studies.
Average SRT only slightly affected gene removal in multivariate analysis and was not significant in the univariate analysis. The SRT is optimized in each plant to achieve the best conditions for nutrient removal (Smith et al., 2014). Although higher SRTs have shown to improve the removal of pharmaceuticals (De Sotto et al., 2016), its effect on ARB and ARGs removal is still controversial and restricted to bench scale studies (De Sotto et al., 2016;Neyestani et al., 2017). Higher SRT might favor the grazing of bacteria by protozoa, but on the other hand, it might also favor HGT events (Tsutsui et al., 2010). Moreover, we also evaluated the effect of primary sedimentation processes. These are meant to reduce debris, TSS, and BOD by mechanical and/or settling procedures (Puig et al., 2010). Bacteria associated with such particles might be removed during primary settling. However, such an effect may be masked by primary settlers receiving recirculated water from thickeners and dewatering digested sludge that contains high amounts of ARGs (Gao et al., 2012). In our study, the presence of primary sedimentation with and without recirculation of digested sludge seemed to result in slightly decreased ARG removal. This appeared from the inclusion of this parameter in two best models, retrieved on a subset of the data with complete observations, albeit these parameters not being significant in themselves (Table 4 and  Table S5). However, when this statistical model was repeated on the full dataset, the effect of primary sedimentation with and without recirculation of digested sludge was less pronounced (i.e., the estimates were numerically smaller and the p-values higher) (Table 5  and Table S6). Lastly, the effects of the remaining design and operational WWTPs parameters investigated (Table 2) were not statistically significant, namely the size of the plant and the type of P removal (chemical or biological) (Table 4).
Thus, the increased hydraulic load caused by rainfall remains the single clear parameter of the dataset determining the ARG removal. The simplest model describing resistance gene reduction, therefore, includes the hydraulic load only ( Table S7 in the supplementary information). According to this model, ermB, sul2, qnrS, and tetM are removed significantly more efficiently than 16S rRNA, while korB increases in relative abundance. The effect of the hydraulic load on the removal efficiency differs per gene: qnrS, tetM, and ermB are significantly better removed at higher hydraulic load as compared to 16S rRNA (À0.69 to À0.83 logs/per time ADF was doubled) (Table S7).

The challenge of comparing and anticipating treatment efficiencies
Albeit an increasing number of studies addressing AMR in WWTPs in the last decade, no conventional treatment or operational strategy has been identified that can improve ARG removal. Unlike in laboratory approaches, full-scale studies involve dozens of variables at once. Effects caused by parameters that are targeted in a specific study can be masked by others (environmental, design). As the majority of the available studies include relatively few locations (from 1 to 5 WWTPs), meta-analyses might be used to aggregate data from single studies and more sensitively identify explanatory factors for ARG removal. However, some studies do neither gather nor include crucial metadata about plant design and operational parameters along with the sampling campaign. Preferably, the collected parameters should be specific for the sampling dates rather than representing average values. In any case, cooperation with water authorities in both the sampling design and the evaluation of the results might ease the access to operational process information. Furthermore, the comparison of results between studies is not always possible since the ARGs assessed often differ. This might be helped by a consensus panel of ARGs, such as suggested by Berendonk et al. (2015) which was used in this work. Additionally, not all studies report both reductions in absolute concentrations and reduction of the relative abundance of ARGs normalized to the 16S rRNA gene. The absolute ARGs concentrations in influent or discharged effluents provide valuable information for risk assessments and to estimate the plant performance. On the other hand, relative abundance is relevant to point to possible selective processes within the plant. Integrating each of the aforementioned points into future studies would help to build a comprehensive data frame that might result in a better understanding of the efficiency of ARG removal in wastewater treatment.

Conclusions
From an analysis of the influent concentrations and the removal of ARGs in a large number of WWTPs in the Netherlands, we conclude that: Table 4 Effect of operational parameters on the reduction of ARGs and MGEs (univariate models excluding korB). Models give the effect of single operational parameters, adjusted for the Hydraulic Load Factor and its interaction with the gene identity. Acronyms: n: number of plants with available information for that parameter; beta: model estimate; SE: standard error of the estimate; z: z statistics; p: p-value; NP: no primary clarification; PNR: primary clarification without recirculation of activated sludge; PR: Primary clarification with recirculation of activated sludge; B: biological phosphorus (P) removal; C: Chemical P removal; BC: Biochemical P removal. HRT: Hydraulic Retention Time; SRT: Sludge Retention Time; ACT: Anaerobic Contact Time. The model for the Hydraulic Load Factor includes the Hydraulic Load Factor, the gene, and their interaction only. In all models, ermB, tetM, and qnrS are significantly better removed than 16S. In some models, sul2 is also significantly better removed than 16S. An increase in the Hydraulic Load Factor leads to significantly reduced removal for ermB, tetM, and qnrS as compared to 16S in all univariate models. Rainfall causing an increase in the usual WWTP hydraulic load marginally reduced the amount of incoming ARGs but significantly reduced the WWTP's removal efficiency of the ARGs and intI1. WWTP design parameters as size, presence of primary clarification, type of P removal and operational parameters as HRT, SRT, anaerobic contact time and effluent TSS were not found to affect the removal of the studied ARGs and MGEs significantly. However, the use of average annual data instead of actual data on the sampling day for some of these parameters probably masked their possible effect.

Acknowledgments
We would like to thank all the Dutch waterboards and WWTP operators for their collaboration in the sampling collection. We also warmly acknowledge Dr. Sven Jechalke and Dr. Kornelia Smalla Table 5 Effect of operational parameters on the reduction of ARGs and MGEs (excluding korB) e results of multivariate modeling. The estimates of the effects of explanatory variables on gene removal (on log 10 scale) during WWTP treatment (p values are shown in italic between brackets; significant estimates are shown in bold) are shown for the 11 best models that were of nearly identical quality, as determined by AICc. For each model, both the results based on the subset of 37 plants for which all parameters were known and the results for the larger subset of plants for which these specific model parameters were known are given. The number of plants included in the model is shown under "n". Gene identity was also included in all models, as was the interaction between the Hydraulic Load Factor and the gene identity. Genes for which significant interactions with the Hydraulic Load Factor were found are listed (for ermB, qnrS, and tetM, the decrease in log removal with increasing Hydraulic Load Factor was lower than for 16S rRNA). Acronyms: n: number of plants in the model; p: p-value PNR: primary clarification without recirculation of activated sludge; PR: Primary clarification with recirculation of activated sludge; N: not included in the model. from the Julius Kühn-Institut in Germany for providing the controls for the IncP-1 (korB) qPCR assay. Finally, we would like to thank Adrian Gonzalez Fraile for his help during the setup of the korB assay. This work was partly funded by the Ministry VWS to RIVM. Wetsus is co-funded by the Dutch Ministry of Economic Affairs and Ministry of Infrastructure and Environment, the European Union Regional Development Fund, the Province of Fryslan 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. The authors would also like to thank the members of the research theme Source Separated Sanitation for the shared knowledge and financial support.