Unraveling the riverine antibiotic resistome: The downstream fate of anthropogenic inputs

are one of the main routes by which the


Introduction
Antibiotic resistance is increasingly recognized by international and governmental entities as a growing global public health threat.According to a 2014 report by the Wellcome Trust and the British government, more than 50,0 0 0 cases of antibiotic resistant infec-supply, irrigation, and recreation.Research in this area has greatly intensified over the last decade ( Bürgmann et al., 2018 ;Rizzo et al., 2013 ;Zhang et al., 2009 ) and an increasing number of studies have investigated anthropogenic impacts on receiving rivers.Among the earliest investigations were the studies on the Poudre River in Colorado, United States, which proposed quantitative polymerase chain reaction (qPCR)-based quantification of various antibiotic resistance genes (ARGs), along with phylogenetic analysis (e.g., tetW ), as a framework for tracking anthropogenic inputs.Anthropogenic input of ARGs to the receiving river was well-apparent using this approach ( Storteboom et al., 2010 ).More recently, the advent of shotgun metagenomic sequencing has greatly advanced the resolution in the ability to characterize large-scale impacts of anthropogenic ARG inputs, as was observed in the Han river catchment in Korea ( Lee et al., 2020 ).The authors noted a strong association of fecal contamination as evidence of anthropogenic activities shaping the composition of the downstream antibiotic resistome (collective ARGs in a microbial community).In Switzerland, a study on rivers and lakes identified the occurrence of extended spectrum β lactamase-and carbapenemase-producing Enterobacteriaceae , which presumably originated from anthropogenic activities ( Zurfluh et al., 2013 ).Another recent study revealed that stream microbiota are significantly altered by the input of treated wastewater in natural streams ( Mansfeldt et al., 2020 ).
ARGs and resistant bacteria (ARB) can persist or proliferate in environmental systems by various mechanisms.Horizontal gene transfer may occur, potentially resulting in new combinations of ARGs or the transfer of resistance to environmentally-adapted bacteria that could in turn change the role of the environment as reservoirs of resistance for clinically-relevant bacteria.Furthermore, the possibility of resistance selection under sub-inhibitory concentrations of antibiotics has been reported ( Andersson and Hughes, 2012 ).Recently, first attempts have been made to estimate predicted no effect concentrations (PNECs) for resistance selection in environmental settings ( Bengtsson-Palme and Larsson, 2016 ).However, current PNECs are an estimate extrapolated from data on isolated bacteria, and could vary substantially under in-situ environmental conditions and with environmental bacteria.
The above examples make clear that treated wastewater discharges have a significant impact on the abundance and types of ARB and ARGs in receiving rivers.Thus, it is crucial that we gain a better understanding of the downstream fate of the anthropogenic antibiotic resistome in receiving rivers.In this sense, few previous studies have attempted to investigate the downstream behavior of various indicators of resistance (e.g., ARBs, ARGs, mobile genetic elements commonly associated with ARGs), and no clear picture of such behavior has as yet emerged.For instance, a study performed in two wastewater treatment plants (WWTPs) and their receiving river in China reported that the levels of wastewater-origin ARGs ( tetC, sul1 ) and the class 1 integron integrase gene ( intI1 ), decreased significantly 1.2 ~2.5 km downstream of the wastewater discharge point ( Li et al., 2016 ).On the other hand, a study performed in a Dutch stream showed that the downstream levels of sul1, sul2, ermB, tetW , and intI1 persisted, or even increased for certain genes over a 20 km downstream distance ( Sabri et al., 2018 ).Mass-flow analyses of ARB and ARG are missing.These contradictory results regarding the downstream behavior of resistance determinants could be in principle attributable to various factors -different geo-hydrological conditions, potential inputs from nonpoint (e.g.agricultural) sources, and the possible existence of biological drivers (i.e., horizontal and/or vertical gene transfer).An improved understanding of the fate of the wastewater-origin antibiotic resistomes and underlying causes would therefore require an integrated approach across multiple disciplines.
The purpose of this study was to track wastewater-origin antibiotic resistomes and identify the key mechanisms governing their fate in two of the most substantially wastewater-impacted rivers in Switzerland.It was hypothesized, that short-distance (up to 1~2 km from wastewater discharge point) behavior of wastewaterorigin resistance determinant concentrations would be governed mostly by hydrological effects such as mixing and dilution.Thus, we used conservative chemical tracers to determine dilution effects and further investigated the contribution of dilution on downstream dynamics of resistance determinants.On the other hand, it was expected that over longer distances (more than 1~2 km; up to 13.7 km to the next downstream WWTP) fate of ARGs and ARB would depend also on additional source/sink mechanisms, such as from biological processes (e.g., death or growth of wastewaterorigin ARB, in-situ resistance (co)selection by antibiotics or metals, and horizontal and/or vertical gene transfer), and/or non-point sources (e.g., agricultural runoff) which are expected to diffuse into the system continuously from a large catchment area.Therefore, the potential effects of biological drivers over long downstream distances were investigated after accounting for hydrological effects.To provide a comprehensive assessment of indicators of antibiotic resistance in the environment, we combined various approaches: cultivation of heterotrophic bacteria on media containing antibiotics, quantification of key indicators of anthropogenic sources of antibiotic resistance by qPCR, and broad profiling of the resistome in selected samples using shotgun metagenomic sequencing.Our study contributes to a systematic, interdisciplinary understanding of the mechanisms driving the fate of the wastewater antibiotic resistome in anthropogenically-impacted rivers.

Site description and field work
A list of WWTP effluent-receiving Swiss rivers without known upstream point-source inputs (e.g., other WWTPs) was obtained from a database provided by Eawag, the Swiss Federal Institute of Aquatic Science and Technology (retrieved 2018) ( Eawag, 2014 ).Two sites were selected according to the following criteria: 1) Greatest proportion of effluent discharge to river discharge, 2) Least number of side streams (for minimum dilution effect from side streams), and 3) Longest distance until the receiving river reaches another downstream WWTP.The selected sites were the river Suze in Villeret (VIL) in canton Bern, and the river Murg in Münchwilen (MUE) in canton Thurgau (maps with all sampling points, see supplementary Fig. S1 and S2).At the sampled sections, both are shallow (generally < 30 cm depth under low flow conditions, maximum depth 1 m) rivers of Strahler order number 3 and 6, and a mean annual runoff of 2.03 and 1.61 m 3 /s, respectively.The river beds are mostly gravel.To avoid elevated flow conditions we sampled only under dry weather conditions at the time of sampling and during at least the previous 36 hours.
At VIL, we studied a 23.7 km stretch of the Suze that we sampled from 10 km upstream (US5) of the effluent (EF) discharge point of WWTP Villeret and at 8 downstream sites located from 0.5 km (D1) to 13.7 km flow distance downstream (D8) before the Suze reaches the discharge point of another WWTP.Four sampling campaigns were performed in 2018 on July 09 (VI L1), July 19 (VIL2), July 30 (VIL3), and November 05 (VIL4).Different combinations of locations were sampled in each sampling campaign as described in detail in the SI (pp.4-5).Daily discharge measures from two gauging stations, one near US, and the other near D8 were obtained and are given in Dataset S1.
At MUE, we studied a 7.0 km stretch of the Murg that we sampled from 0.2 km upstream (US) of WWTP Münchwilen and at 8 downstream sites located from 0.5 (D1) to 6.8 km flow distance downstream (D8 before the Murg reaches the discharge point of another WWTP.Three sampling campaigns were per-formed in 2018 on July 26 (MUE1), August 03 (MUE2), and August 06 (MUE3).Discharge data was obtained from a gauging station near D4 (Dataset S1).
Samplings were performed according to other projects performed in Swiss rivers and WWTPs ( Mansfeldt et al., 2020 ;Ju et al., 2019 ).At each river sampling location grab samples (5L in sterilized water containers) were obtained by combining water from just below the surface at three points along a river transect: in the middle and roughly equidistant from the banks to each side of the middle point.EF samples were obtained from the final effluent of the WWTPs prior to discharge.Temperature ( °C), conductivity ( μS/cm), pH, and dissolved oxygen (DO, mg/L), were measured on site in an aliquot of the sample using a portable multiparameter probe (Multi 3630 IDS, WTW, Germany) at the time of sampling.To make sure EF was fully mixed with receiving water at D1, the conductivity values across the cross-section were measured, and no significant deviation was observed ( < 0.5%).All samples were cooled at 4 °C in the dark while transported to our laboratory on the same day.Samples were processed on the same and next day within 36 hours.For organic micropollutant analysis, water samples were obtained separately, and stored in precombusted glass bottles on site, cooled at 4 °C during transportation, and frozen at -20 °C in the dark upon arrival at the laboratory until analyzed.Sediment samples were obtained from 5 select locations (US, D1, D2, D5, and D8) for select campaigns (VIL1~3 and MUE1~3), and frozen at -20 °C upon arrival at our laboratory.
To better constrain flow velocities in the rivers, salt tracer experiment using NaCl and flow-velocity measurement were performed in separate sampling campaigns in August 2019 (August 23, 2019 for VIL, and August 27, 2019 for MUE) under comparable flow conditions.The results were summarized in Dataset S4 (e.g., flowvelocity and hydraulic residence time).
Further details on sites, field sampling procedures and hydrological experiments are given in the Method section of the SI (pp.4-5).The exact sampling locations (GPS-coordinates) are given in Table S1.

DNA extraction and quantitative PCR
Two aliquots of each water sample were filtered through two 0.2 μm pore size membrane filters, using 0.5 L for EF and 1.0 L for river water samples.Replicate filters were then processed separately.DNA was extracted from the filters using DNeasy Power-Water Kit (Qiagen, Germany) following the manufacturer's instructions.For sediment samples, DNA was extracted from about 20 g of wet sediment using the DNeasy PowerMax Soil Kit (Qiagen, Germany).Extraction blanks confirmed absence of DNA contamination (see SI, p.8).The concentration and qualities of extracted DNA were measured using a NanoDrop One spectrophotometer (Thermo Fisher Scientific, USA) (Dataset S2).
Presence and abundance of key indicator genes for anthropogenic ARG inputs ( sul1, tetW, ermB, bla CTX and integron integrase class 1 gene intI1 ) ( Berendonk et al., 2015 ;Gillings et al., 2015 ;Ju et al., 2019 ), were determined by qPCR as described previously ( Czekalski et al., 2012 ;Czekalski et al., 2014 ).The detailed qPCR protocols are reported in the SI.Absence of contamination from filtration and extraction procedures was confirmed using an experimental control by qPCR analysis as shown in the SI.

Metagenome and 16S rRNA gene amplicon sequencing analysis
Shotgun metagenomics and 16S rRNA gene amplicon sequencing analysis were performed using Illumina platforms for samples from three selected sampling campaigns.Samples were selected for sequencing according to the following rationale: For VIL the samples were selected only from campaign VIL1 (i.e., 6 samples: US, EF, D1, D2, D5, D8) as the samples from other campaigns (VIL2-3) showed similar patterns of resistance determinant profiles downstream.For MUE, 6 samples from MUE2 and MUE3 campaigns (i.e., US, EF, D1, D3, D5, D8 from each sampling) were selected as the far downstream behaviors of certain ARG (e.g., sul1 ) and intI1 were significantly different from each other in those campaigns.DNA extracts from replicated filters were pooled.All library construction and sequencing was performed by Novogene (Hong Kong).
To analyze 16S rRNA gene amplicon sequencing data, we used the DADA2 pipeline ( Callahan et al., 2016 ), and followed the workflow suggested by the developers.The detailed protocol is described in the SI.

Chemical analysis
Metals, ions (i.e., dissolved cations and anions), nutrients, and dissolved organic carbon were measured as described in Ju et al. (2019) using high-resolution inductively coupled plasma mass spectrometry, ion chromatography, flow-injection analysis, and total organic carbon analyzer, respectively, as described in the SI.Dissolved micropollutants (i.e., pharmaceuticals, antibiotics) were measured as described in Ju et al. (2019) using liquid chromatography triple quad mass spectrometry with electrospray ionization in the SI.Total dried solid (TS) were measured in sediment samples according to standard methods ( APHA-AWWA-WPCF., 1981 ).

Estimating the dilution effect on downstream levels of resistance determinants
Under continuous discharge and after complete horizontal and vertical mixing, the discharged load of a conservative tracer (e.g., sodium) entering the river through EF is expected to be conserved along the river continuum.Under this assumption, any change in the concentration of the conservative tracer would be due to dilution effects by additional water inflows (i.e., groundwater and/or tributary inputs) and additional inputs of the tracer with these inflows.We used sodium and two micropollutants as conservative tracers (i.e., 4/5-methylbenzotriazole, carbamazepine) because these substances had high concentrations in EF compared to the US river and are known to not substantially degrade or adsorb in the river system.The rationale for selecting the conservative tracers is described in more detail in the SI.
Starting with these mass conservation assumptions, under steady state conditions, the dilution parameter ( DP , the ratio between external water inflow and streamflow at the downstream section between any two points A and B along a river stretch) can be estimated from a ratio of tracer concentrations according to Eq. 1 : where A indicates an upstream location; B denotes a downstream location; C indicates the concentration of a tracer; C denotes the average concentration of the tracer in the external inflow between A and B.
The derivation of Eq. 1 is schematized in Fig. 1 , and also described in detail in the SI.Concentrations C A or C B were measured directly for all compounds, C was estimated according to the following equation ( Eq. 2 ) for sodium (the values shown in Dataset S9.3) and assumed to be 0 for 4/5-methylbenzotriazole and carbamazepine.In short, the difference in sodium loadings (mass per time) between the point of EF discharged and the downstream point where gauging stations were located (D8 for VIL; D4 for MUE) was divided by the quantity of additional water inflows: Where, N a D 8 or D 4 denotes the sodium concentration measured at D8 or D4; Q US, EF, D 8 or D 4 indicates the river flow quantity or wastewater effluent discharge (volume per time) at US, EF, D8 or D4.
The value DP should be the same for all conservative tracers.To test our hypothesis that "the short distance dynamics of resistance determinants is largely governed by dilution effects", we calculated DP over short distance (i.e., DP D1 → D2 for VIL, and DP D1 → D3 for MUE), and compared DP values over the same distance for resistance determinants ( sul1, intI1, ermB, tetW, and CLR/TET resistant bacteria) with values for the conservative tracers.Higher DP values for resistance determinants would indicate a lower than expected concentration in the downstream and thus removal.
The expected downstream concentrations of resistance determinants considering dilution as a main driver can be calculated using the DP A → B of conservative tracers according to the following rela-tionship: (3) where C resist−A indicates the concentration of a resistance determinant at an upstream location A; C resist−B denotes the concentration of a resistance determinant at a downstream location B; D P A → B ( for X ) indicates the DP of a conservative tracer (X) between A and B Eq. 3 assumes that resistance determinants behave conservatively over the studied distances and that there are no significant inputs of resistance determinants from the diluting water inflows (i.e., C for resistance determinants ≈ 0 in Eq. 1 ).Therefore, deviations from measured to predicted values can indicate violation of these assumptions.We calculated the predicted concentration by dilution effects for each resistance determinant under these assumptions for all downstream sections of the rivers.

Estimating the river discharge over downstream distance
The river discharge ( Q ) at was estimated for several downstream locations where there were not gauging stations.The estimated Q values were used when calculating loadings of chemical and resistance indicators over downstream distance.The Q EF values were obtained from each WWTP, and Q US values were either obtained from gauging station (for VIL), or calculated as shown in Eq. 8 in the SI (for MUE).
Where, Q D (n ) indicates the river discharge ( Q ) at the downstream location D(n) (2 ≤ n ≤ 8); D P D ( n −1 ) → n denotes the dilution parameter between D(n-1) and D(n).

Upstream water quality and WWTP effluent
In agreement with the criteria for site selection, the levels of intI1 and target ARGs upstream of the WWTP were generally low, except for the MUE2 campaign where we observed elevated upstream levels of ermB , and intl1 ( Fig. 2 ; Fig. S4 and S5 in the SI).Chemical water quality likewise did not suggest significant pollution inputs from either tributaries or upstream locations for either VIL or MUE as most micropollutants were below the limit of quantification (Dataset S10~11).Certain micropollutants (e.g., 4/5methylbenzotriazole, benzotriazole, and diclofenac) were sporadically detected in very low quantities.For cultivable multi-resistant bacteria ( Fig. 3 ), especially CLR/TET resistance, relatively high upstream values were observed in VIL2 and MUE2 US samples.These findings indicate that while there is no indication of significant upstream pollution, some pollution, probably from periodical urban or agricultural activities, may affect the river.There was a settlement upstream of the WWTP and livestock farming activities (i.e., pastures and meadows for livestock) in the catchments, including the upstream sections, in both sites ( BAFU, 2013 ).While we assume surface runoff from the agricultural sites to be minimal as our samplings were performed under dry-weather conditions, it cannot be ruled out that some inputs from agricultural activity occasionally affected the river.Further investigations into the nature of these transient microbial contaminations were not carried out in this project, but future work could employ microbial source tracking or microbial fingerprinting approaches to determine their sources.The effluent from both WWTPs contained considerable levels of pollutants.For instance, effluent concentrations were higher than the upstream levels by approximately 1 order of magnitude for sodium, 1~2 order of magnitude for ARGs and intI1 , and more than 2 orders of magnitude for micropollutants (Datasets S8~11).These results are in line with previous results from a large-scale investigation of micropollutants in Swiss streams ( Stamm et al., 2016 ).

Short range fate of antibiotic resistance determinants in the downstream river
Focusing on the immediate impact of the WWTP effluents (US versus D1 to D3 sites), there were significant impacts of WWTP effluents on the receiving rivers in both VIL and MUE.The estimated proportions of EF in the downstream receiving waters (D1) estimated by conductivity were 10.5 ~35.9 % for VIL1~4, and 33.0 5 38.0 % for MUE1~3 (Dataset S9.2).Accordingly, significant increases of sul1, ermB, tetW and intI1 as quantified by qPCR were observed at D1 compared to US (p < 0.01 paired t-test; Fig. 2 & Fig. S5).However, the measured levels of these antibiotic resistance indicator genes rapidly decreased nearly to upstream levels over 2.5 and 2 km downstream distance (D2 or D3 locations) in VIL and MUE, respectively.
The same dynamic was also observed for multi-resistant bacteria ( Fig. 3 ), especially CLR/TET resistance.SMX/TMP/TET resistance was often below the limit of detection (5.0 CFU/mL), but clearly exceeded it in the D1 samples and was thus also higher there than further downstream (from D2 on).
Several processes may contribute to the observed decrease of resistance determinants, including dilution by additional water inflows via groundwater and/or tributary inputs, biological deterioration (e.g.cell death or dormancy due to exposure to sunlight, lower ambient temperature, predation, etc.), and cell sedimentation.

Dilution effects strongly affect short distance dynamics of effluent resistance determinants
To determine the importance of dilution effects, we compared DP calculated over a short distance (D1 to D2 for VIL; D1 to D3 for MUE) downstream of the WWTP discharge point ( D P D 1 → 2 for VIL; D P D 1 → 3 for MUE) from conservative chemical tracer concentrations (e.g., sodium, 4/5-methylbenzotriazole, carbamazepine) as well as ARG and intI1 levels ( Fig. 4 ).The average D P D 1 → 2 or 3 of the target antibiotic resistance indicator genes levels were always higher than for conservative tracers, indicating possible removal mechanisms at play.However, according to the paired t-test under the null-hypothesis of "No significant differences of dilution parameters between different pairs of bio-and conservative indicators", only the differences between sul1 and tetW versus sodium were significant at p < 0.05 (p-adjusted using Benjamini-Hochberg method) ( Fig. 4 ), confirming non-conservative behavior and additional removal mechanisms.As D P D 1 → 2 or 3 for sodium took up a large portion of the values for sul1 and tetW ( i .e .DP ( N a + ) D 1 → 2 or 3 = 0 .72 ∼ 0 .92 × DP ( sul1 ) D 1 → 2 or 3 and 0 .59 ∼ 0 .96 × DP ( tetW ) D 1 → 2 or 3 ), the dilution effects quantified by sodium nonetheless explain the majority of the concentration decrease for these parameters.This result implies that the observed rapid decrease in the downstream levels of wastewater-origin resistance determinants immediately downstream of the WWTPs was mainly governed by dilution in the studied systems.Dilution effects thus need to be carefully considered in studies of the environmental fate of resistance determinants, and loadings instead of concentrations need to be determined to accurately assess environmental fate.

Additional source/sink effects become apparent over longer downstream distances
We hypothesized that additional source/sink mechanisms affect the downstream behaviors of antibiotic resistance indicator genes over longer distances.To analyze this in more detail the daily loading (copies/day) for the target ARGs and intI1 at the point of discharge (as the sum of upstream and EF loadings), and for short (D2 for VIL; D3 for MUE) and far downstream distances (D8) were calculated by multiplying resistance levels (copies/m 3 ) with the discharge (m 3 /day) at each location and then compared.The discharge was either obtained directly from nearby gauging stations, or estimated under consideration of the EF discharge (m 3 /day) using sodium as an indicator, and according to the equation derived under the mass-conservation assumption ( Eq. 4 ).
The downstream behaviors of the target antibiotic resistance indicator genes varied by indicator and also by sampling campaign.For instance, the load decrease from wastewater discharge (US + EF) to the furthest downstream point (~D8) was pronounced and consistent for ermB and tetW in all the samplings ( Fig. 4 b & c).The average load reduction was 81 ±17 % for ermB , and 70 ±15 % for tetW over 13.7 km distance in VIL1~4; 95 ±5 % for ermB , and 96 ±2 % for tetW over the 6.8 km distance in MUE1~3 (Dataset S13).In contrast, the downstream behavior of the sul1 loadings was inconsistent between sampling campaigns.A pronounced decrease over distance (64-94 % at D8) was seen in VIL1 ~4, but little reduction over distance in MUE1 ~2 (7 and 29% at D8), and a strong increase in MUE3.The downstream fate of intI1 was also variable, for instance, as intI1 loads did not decrease and in some instances even increased .
To further analyze if there is a break point where sul1 and intI1 start to deviate from conservative behavior, we calculated the predicted levels of resistance determinants over the whole study distance considering dilution as a major driver using Eq. 3 (Dataset S8).The predictions based on three different conservative tracers are visualized for sul1 and intI1 in Fig. 2 .In VIL, measured levels are always below predictions, except for intI1 in VIL2.For MUE, in contrast, we see measured values exceeding predicted values in several instances, for intI1 even for most downstream locations.In MUE3 where the pronounced increase of sul1 loading was observed between D5 and D8, the level of sul1 started to exceed predicted values between 5 ~6.8 km distance.The concentration of intI1 increased also very rapidly between 5 ~6.8 km downstream distance in MUE3, which indicates either a pronounced proliferation or a non-point source of sul1 and intI1 in this stretch of the river.We will discuss potential mechanisms (e.g., biological drivers, on-site selection, additional anthropogenic source input, and surface sediment inputs) in section 3.10 .
A number of mechanisms may contribute to the generally observed removal: Sedimentation (especially of cell aggregates or flocs) and cell death by predation, UV light, or various other environmental conditions unfavorable to wastewater bacteria.With the available data we are not able to determine the contribution of various mechanisms.Future studies could investigate the persistence of resistant bacteria or molecular resistance markers in micro-or mesocosm experiments or in a turbulent flow system mimicking natural streams to answer such questions.Modeling transport and sedimentation of wastewater-origin particles using the information on particle size, mass, and flow characteristics could provide information on the importance of sedimentation.

WWTP effluent affects the downstream riverine resistome
To obtain a broader view of the river antibiotic resistome we retrieved the ARG content of metagenomes obtained for selected sampling campaigns (VIL1, MUE2, and MUE3).Overall, 65 ARG subtypes were identified, 49 of them occurred in both upstream and downstream river samples (Dataset S7).The antibiotic resistome in the receiving water closest to the discharge point (D1) was clearly influenced by the input from EF.For instance, a total 28 out of 36 ARG subtypes found in D1 were also detected in EF (B, C, F, G in Fig. 5 a) while 16 of these were not observed US (B, C in Fig. 5 a).The 16 EF-derived resistance genes confer resistance to the following antibiotic classes: 1 x aminoglycoside, 4 x beta-lactam antibi-otics, 1 x chloramphenicol, 2 x macrolide, 1 x quinolone, 4 x tetracycline; 3 subtypes were multidrug resistance genes.Of these 16 genes, 14 were no longer detected in the far downstream (6.8 13.7 km downstream distance in MUE and VIL, designated as D_Far in Fig. 5 a).This is in agreement with the results for qPCR enumeration of ARGs and intI1 reported above and implies that the majority of ARGs that occurred exclusively in EF do not persist at detectable levels for a long distance in rivers where significant amounts of additional water inflows and additional removal mechanisms are expected.

Diversity of the river resistome and microbiome along the river continuum
We analyzed the alpha-, and beta-diversity of river resistomes and microbiomes as another way to observe potential effects of the WWTP discharge and to see if the dynamics in the resistome are strongly correlated with the microbial community, as noted e.g. for changes observed during wastewater treatment ( Ju et al., 2019 ).As expected, Shannon alpha-diversity of ARGs was higher in EF compared to US samples by 20.2 ~225.4 % ( Fig. 5 b).Accordingly, the impact of the EF resistome was observed, especially for VIL1 and MUE3, as an increase in ARG diversity in river water at the D1 sites.The ARG alpha-diversity decreased downstream in all sampling campaigns.However, for the microbial community as represented by amplicon sequence variants (ASVs), alpha-diversity was not consistently higher in EF versus US, and consequently also did not change significantly from US to D1 and did not consistently decrease downstream ( Fig. 5 b).This indicates that the downstream dynamics of the overall microbial community and the antibiotic resistome were decoupled.
Similar conclusions were obtained from beta-diversity analysis.Procrustes analysis between microbial communities and antibiotic resistome ( Fig. 5 c) revealed a strong structural correlation between microbial communities and antibiotic resistome only when the most strongly effluent-affected sites (D1) were consid- ered (p = 0.002; Fig. 5 c).When the D1 samples were excluded, the correlation was barely significant (p = 0.04; Fig. 5 d), indicating that the structural correlation between microbial communities and resistome largely resulted from the impacts of WWTPs on D1.The weak structural correlation in less impacted waters suggests a lack of strong drivers, such as selective pressures or the influx of external ARB.

Resistome analysis confirms effluent effect and abatement
To quantitatively investigate the dynamics of the resistome along the river continuum in more detail, the seven most prevalent ARG subtypes that appeared in more than 9 out of 18 samples were chosen for detailed analysis.We calculated the proportion of each gene in this set based on relative abundance (GPM, gene per million) ( Fig. 6 a).The proportions (%) of each of six genes ( aph(3 ')-I, aph(6)-I, mexT, tetQ, aadA , and sul1 ) to the whole seven genes were lower in US than in EF and D1 ( Fig. 6 a).The bacA gene, in contrast, comprised a larger proportion in US (i.e., up to 83 % in VIL1:US), D5 and D8 (i.e., up to 93 % in VIL1:D8) than in EF and D1.It was therefore excluded from the plots of relative and cumulative abundance of the assembled ARG in the metagenome in Fig. 6 b, and their individual and cumulative environmental level (Gene per liter) in Fig. 6 c.Both relative and absolute abundances showed a similar pattern -the abundances of the selected ARG were higher in EF and D1 compared to US ( Fig. 6 b).The abundances of those six genes decreased along the downstream locations except for sul1 in the far downstream location (D8) in the MUE3 campaign ( Fig. 6 b).This analysis confirmed a quantitative effect of the effluent on the abundance of prevalent resistance genes in the river resistome, and suggests additional candidates for tracking anthropogenic sources of resistance in future studies ( aph(3 ')-I, aph(6)-I, mexT, tetQ, aadA , and sul1 ) that may be useful for tracking resistance determinants from wastewater.Several of these genes have been used as resistance indicators for environmental samples mainly in combination with culture-dependent approaches ( Rizzo et al., 2013 ;Zhang et al., 2009 ), but much less frequently with culture-independent approaches ( Rizzo et al., 2013 ;Sharma et al., 2016 ).

Contig analysis indicates ARG co-location
We hypothesized that the similar dynamics of some ARGs could derive from their co-location in the same host or presence on the same genetic elements, so we analyzed their loci within the contigs.aph(3 ')-I and aph(6)-I were indeed found to be located on the same contig in many samples from VIL1 (EF, D1, D2, and D5) ( Fig. 7 a), which may explain the strongly similar dynamics between aph(3 ')-I and aph( 6)-I GPM in VIL1 with R 2 = 0.98 (p < 0.001) ( Fig. 6 b).Another case of co-location between ARGs was observed between sul1 and aadA .The contigs containing those two genes were found in D1 in VIL1, and EF in MUE3 (Fig. S13c; Fig. 7 b).Unlike aph(3 ')-I and aph(6)-I in VIL1 however, the dynamics of sul1 was not similar to that of aadA especially in D5 and D8 in MUE3 where high abundance of sul1 was observed while no aadA was assembled ( Fig. 6 b).In agreement with this observation, the sul1 -containing contigs retrieved from D5 and D8 in MUE3 did not contain aadA , and this was the only type that was identified in those samples ( Fig. 7 b).This indicates a significant shift of the bacterial populations that yielded sul1 -containing contigs between D1 and D3 in MUE3, with wastewater-derived contigs containing sul1 -aadA pairs not persisting.

bacA, an ARG with high natural prevalence in environmental bacteria
Unlike the other prevalent genes, the proportion of bacA (also known as UppP , undecaprenyl-diphosphate or -pyrophosphate phosphatase) was greatest in US samples, and was also abundant in many downstream locations (except D5 and D8 in MUE3 where sul1 occupied the largest proportion) ( Fig. 6 a).In order to further assess whether bacA was intrinsic in our river water samples, we identified potential hosts by assigning taxonomy to the metagenome-assembled contigs using a combination of methods.The contigs for which all three methods agreed at the genus level are shown in Table 1 .The four genera identified as potential hosts of bacA contigs derived from less-disturbed freshwater samples (US, D2, D5, D8) were Pseudomonas, Acidovorax, Limnohabitans , and Aeromonas .Among them, Pseudomonas, Acidovorax , and Limnohabitans are typical inhabitants of freshwater and soil environments ( Peix et al., 2009 ;Willems, 2014 ).However, considering that the proportions of bacA in the contigs to the total bacA in the sample in terms of reads per kilobase (RPK) were low for river water samples (except for D8 in MUE3), we assume that homologues of Metagenome-assembled bacA -containing contigs to which taxonomy was successfully assigned at genus level.Taxonomy assignment was performed using Kaiju, Kraken2, and the basic local alignment tool for nucleotides (Blastn), and only contigs with consensus from all three approaches at the genus level are shown.For Blastn, the quality criteria were P Ident > 90.0 %, and Q cov > 90 %.P tot_ bacA indicates the proportion of bacA in the contig to the total bacA in the sample in terms of reads per kilobase.P Ident indicates the percentage of identical match.Q Cov indicates the query coverage.bacA could be present in many other environmental bacteria.Thus, our data suggests that bacA is probably unsuited for tracking anthropogenic sources of antibiotic resistance.

Exploring the potential reasons for rapid increase of sul1 in far downstream locations in MUE3
Both qPCR-based, and metegenomic analysis confirmed that sul1 and intI1 increased in the downstream of MUE, and especially strongly in one of the sampling campaigns (i.e., MUE3) between 5.0 -6.8 km downstream distance ( Fig. 2 b & 6 ).
To figure out if there was a biological driver for this unexpected increase of sul1 and intI1 , we first characterized the sul1 -containing contigs.The sul1 gene is known to be highly mobilized and is often associated with intI1 ( Gillings et al., 2008 ;Gillings et al., 2015 ).Indeed, all contigs containing sul1 associated with intI1 retrieved from the river (D3 ~8) were homologs of a single dominant type that appeared to be also plasmid-associated as it contained the plasmid-associated gene parA ( Davis et al., 1992 ).We could unsurprisingly not obtain a meaningful taxonomic assignment for these sequences.It could thus not be demonstrated whether the downstream increase in MUE3 was due to an increase in an EF-derived or an environmental organism or from a local contamination source.Further information could be obtained in future studies by isolation or construction of metagenome-assembled genomes.We therefore turned to chemical indicators to further study the potential for local or non-point sources of contamination as an explanation.
To evaluate non-point-source inputs of pollutants, we chose to evaluate a few micropollutants that may serve as indicators of contamination.Sulfamethazine (also known as sulfadimidine) is used in pig husbandry ( Stoob et al., 2007 ), and mecoprop is a weed control agent used primarily in urban settings in Switzerland ( Wittmer et al., 2010 ).It was assumed that the levels of these pollutants in downstream locations would increase or be persistently high if a pronounced agricultural or urban surface runoff existed, which could accompany resistance genes and bacteria potentially existing in agricultural or urban areas.In VIL, sulfamethazine was below detection (LOD ~0.8 ng/L) in all samples except one US sample, while in MUE there appeared to be a source in WWTP effluent especially during the MUE3 campaign, but the compound was not observed to increase in downstream locations.The concentrations and downstream dynamics of mecoprop varied between campaigns (Fig. S8).For VIL1~3 and MUE1 and MUE3 mecoprop concentrations were low ( < 60 ng/L), while there seemed to be a strong, effluent-associated input for MUE2 and concentrations remained high further downstream ( > 200 ng/L).A slight increase in the downstream range > 5km observed in MUE2 and between 1.0 -2.0 km in MUE3 may be due to fluctuating input of mecoprop from the WWTP effluent.Concentrations did not further increase in the far downstream locations (D8) where the sudden increase of sul1 and intI1 was observed (Fig. S8).Based on these, but also the other analyzed micropopllutants, we found no evidence for significant downstream contamination sources.However, these chemical indicators are not conclusive, as the analyzed compounds were not a comprehensive selection to trace non-point sources in the downstream river section (e.g., from manure or pesticide applications, although these are not very likely under dry-weather conditions).So while we found no evidence for such contamination we can also not conclusively rule them out as an explanation for the marked sul1 and intI1 increase observed for MUE3.
Finally, the potential for in-situ resistance selection in the water was assessed using the concentration of antibiotics and metals in downstream locations in MUE3.Sulfamethoxazole and its derivative (N4-acetylsulfamethoxazole) were the antibiotics with the highest concentration among all the antibiotics analyzed, but downstream concentration (sulfamethoxazole in the range of 33 to 95 ng/L in MUE3) remained far below the published PNEC for resistance selection (e.g., 16,0 0 0 ng/L) ( Bengtsson-Palme and Larsson, 2016 ).The concentration of trimethoprim, which is usually prescribed together with sulfamethoxazole, was also much lower than its PNEC for resistance selection (e.g.500 ng/L) ( Bengtsson-Palme and Larsson, 2016 ).Even though the vast majority of metals analyzed in this study remained below the limit of quantification or below their estimated minimum co-selective concentrations for dissolved metals in water (MCC waterDC ), the concentrations of two metals (i.e., copper and nickel) were higher than their MCC waterDC (1.5 μg/L for Cu, and 0.29 μg/L for Ni) ( Seiler and Berendonk, 2012 ).However, MCC waterDC is a predictive value and actual selective levels could be higher, also their levels in far downstream locations (D8) in MUE3 were not specifically higher than at other locations within the same sampling campaign, nor at the same locations than in other samplings where the increase of sul1 or intI1 was less pronounced (Dataset S12).Furthermore, we did not observe co-localization between sul1 and any other genes potentially conferring Cu, Ni or any other metal resistance based on contig-based co-localization search in MUE ( Fig. 7 b).Overall no convincing evidence for in-situ resistance co-selection by Cu and Ni as an explanation for the downstream increase of sul1 and intI1 was found.We further note that the estimated river retention time per km was relatively short (i.e., 51.4 and 49.5 mins/km for VIL and MUE, respectively, Dataset S4), which makes the likelihood of in-situ resistance selection in the water even less plausible.
As a final possible explanation we considered the possibility of cell migration from other river compartments to the water.According to qPCR enumeration of ARGs and intI1 in surface sediments, we did not observe the increase of sul1 and intI1 levels in D5 and D8 in MUE3 in terms of either absolute and relative abundance (Fig. S6).Furthermore, the relative abundance of sul1 and intI1 in sediment was generally similar to, or lower than the values for water in MUE1~3 (Fig. S7).If sediment resuspension was a major source for aquatic sul1 and intI1 elevation in MUE3 D8, relative abundances of sul1 and intI1 in water samples would be expected to remain unchanged or to drop.While we could not completely exclude the possibility of contribution of sediment resuspension, we assume that there could be other sources (e.g., stream biofilms) where sul1 and intI1 were selectively enriched in terms of both absolute and relative abundance.Considering the downstream levels of both resistance determinants and nutrients remained relatively high in MUE due to high EF inputs and low downstream dilution effects, especially in the third campaign (Fig. S5 & Dataset S9), insystem growth is also a plausible hypothesis.
The reason for and the nature of the striking increase of sul1 and intI1 (but not of tetW, ermB, bla CTX ), during the MUE3 campaign thus remains open and would require further study to resolve.What the observation shows unambiguously, is that unexpected and perhaps not directly anthropogenic contaminationdriven increases of ARGs are possible.As in particular sul1 and intI1 are commonly applied for tracking anthropogenic sources of ARG in the environment ( Berendonk et al., 2015 ;Gillings et al., 2015 ), we caution that monitoring strategies should employ a multi-target strategy to be robust .

Conclusions
• Downstream levels of antibiotic resistance determinants decreased rapidly over 2.0 -2.5 km distance due to dilution effects and decay over longer distance due to other removal mechanisms.This would suggest that public exposure to wastewater-origin antibiotic resistance might be most acute only over short distances (few kilometers) from points of discharge, especially if a pronounced input of additional water inflows exists.• We also observed at least one instance where sul1 and intI1 dynamically increased in the river, without being able to establish any link to a local anthropogenic contamination.Other river compartments where in-system growth of biomass could take place (i.e., stream biofilms) could be included as a monitoring target in future studies.• Metagenomics-based resistome analysis yielded consistent conclusions with qPCR analysis of select targets (e.g., sul1 ) and also identified promising targets for future monitoring of anthro-pogenic sources of antibiotic resistance (e.g., aph(3 ')-I, aph(6)-I, mexT, tetQ , and aadA ).In general metagenomics, qPCR and cultivation-based assays yielded consistent trends.Public health advice could be based on quantifying indicator genes or technically simpler cultivation-based indicators.• A weak structural correlation between resistome and microbiome, and low levels of (co)selective agents revealed a lack of driving forces in less-disturbed river waters (downstream over 3 km distance, plus upstream locations).• We showed that contig-based taxonomic assignment and analysis of the genetic neighborhood of assembled ARG can reveal important, if limited, additional information about shifts in ARG host identities, mobilization, and co-localisation of ARG that would otherwise remain hidden.

Fig. 1 .
Fig. 1.Derivation of dilution parameter (DP) from an upstream point A to the downstream point B using the concentration of a pollute as a marker under the mass-flow assumption.

Fig. 2 .
Fig. 2. Levels (gene copies/mL) of sul1 and intI1 in the upstream near effluent discharge point, and downstream river water quantified by qPCR.Average sul1 and intI1 concentrations in the (a) river Suze near Villeret (VIL) and (b) river Murg near Münchwilen (MUE).The dotted lines are the estimated levels considering only dilution as a major driver according to the Eq. 3 using sodium, carbamazepine (CBZ), and 4/5-methylbenzotriazole (4/5MeBT) as conservative tracers.The point of EF discharged was indicated by a light red vertical line.Symbols indicate the average and tips of error bars are the lower and upper values of biological duplicates.The limit of detection for both sul1 and intI1 is 12.5 copies/mL for all the samples shown here.

Fig. 3 .
Fig. 3. Heterotrophic plate counts (CFU/mL) for clarithromycin and tetracycline multi-resistant (CLR/TET) and sulfamethoxazole, trimethoprim, and TET multi-resistant bacteria (SMX/TMP/TET) from the upstream and downstream river water in (a) Villeret (VIL), and (b) in Münchwilen (MUE).The predicted values were calculated using a selected conservative tracer (i.e., sodium) according to Eq. 3 , and are shown as dotted lines in red, blue, or black.The limit of detection (LOD) was 0.5 CFU/mL for CLR/TET multiresistant bacteria and 5.0 CFU/mL for SMX/TMP/TET multi-resistant bacteria.The LOD for SMX/TMP/TET is shown as a yellow dotted horizontal line.The error bars indicate standard errors among technical triplicates.

Fig. 4 .
Fig. 4. (a) Dilution parameter ( DP ) values over short downstream distance (i.e., D1 to D2 for VIL, D1 to D3 for MUE) among different biological and conservative indicators; (b) ARGs and intI1 loadings at upstream near EF (US), treated wastewater effluents (EF), short downstream (D2 or 3); and long downstream distance (D8) in Villeret (VIL), and (c) in Münchwilen (MUE).The treatment pairs with significant difference in between were asterisked ( * ) in (a).The error bars represent upper and lower values of biological duplicates for each gene in (b ~c).

Fig. 5 .
Fig. 5. Metagenomic analysis of effluent and river antibiotic resistomes at Villeret (VIL) and Muenchwilen (MUE) sites for the selected sampling campaigns (VIL1, MUE2, and MUE3).(a) Venn diagram showing occurrence of antibiotic resistance gene subtypes in the treated wastewater (EF) and in river water upstream (US) and downstream (D1, 0.5 km distance) of the effluent discharge point and in the far downstream (D_Far, 6.8 -13.7 km distance).The presence of ARGs were counted from all three (VIL1, MUE2, MUE3) consolidated-campaigns for each treatment using the presence-absence table shown in Dataset S6.(b) Shannon α-diversity of ASVs (blue) and metagenomeassembled ARG subtypes (red).Procrustes analysis between ASVs (round dot symbols) and resistome (blue arrow tips) where EF was included (c) and EF & D1 were excluded (d).The length of blue arrows indicates the size of Procrustes errors.The error bars represent upper and lower values of biological duplicates for each gene.

Fig. 6 .
Fig. 6.Dynamics of prevalent and widespread metagenome-assembled ARGs along the river continuum.(a) The proportion of each gene among the 7 most frequently occurring and widespread ARGs ( aph(3 ')-I, aph(6)-I, mexT, tetQ, aadA, sul1 , and bacA ).(b) and (c) Stacked bar charts of the abundance of the 6 ARGs that were effluent-associated (omitting bacA ); (b) relative abundance (GPM, gene per million) and (c) absolute abundance (GPL, gene per liter).Sample EF is shaded in red and the other river water samples are shaded in blue.

Fig. 7 .
Fig. 7. Gene arrangement on contigs containing aph and sul1 genes.(a) Contigs containing aph(3 ') and aph(6), retrieved from all samplings (VIL1, MUE2 & 3).(b) Contigs containing sul1 retrieved from MUE3.All annotated genes showed > 90.0 % percent identity (P Ident ) at the protein level to reference proteins, using DIAMOND protein search against NCBI nr protein database.The contig IDs are italicized (e.g.k121_XXXXXX ).P tot_ aph indicates the proportion of average coverage for the aph -containing contig to the sum of average coverages for all the aph -containing contigs identified in the sample.P tot_ sul1 denotes the proportion of average coverage for the sul1 -containing contig to the sum of average coverages for all the sul1 -containing contigs identified in the sample.Only contigs with lengths > 1,0 0 0 bp are shown.