Wastewater bypass is a major temporary point-source of antibiotic resistance genes and multi-resistance risk factors in a Swiss river

Untreated combined sewage (bypass) is often discharged by wastewater treatment plants to receiving rivers during stormwater events, where it may contribute to increased levels of antibiotic resistance genes (ARGs) and multi-resistance risk factors (multi-resistant bacteria and multi-resistance genomic determinants (MGDs)) in the receiving water. Other contamination sources, such as soil runoff and resuspended river sediment could also play a role during stormwater events. Here we report on stormwater event-based sampling campaigns to determine temporal dynamics of ARGs and multi-resistance risk factors in bypass, treated effluent, and the receiving river, as well as complimentary data on catchment soils and surface sediments. Both indicator ARGs (qPCR) and resistome (ARG profiles revealed by metagenomics) indicated bypass as the main contributor to the increased levels of ARGs in the river during stormwater events. Furthermore, we showed for the first time that the risk of exposure to bypass-borne multi-resistance risk factors increase under stormwater events and that many of these MGDs were plasmid associated and thus potentially mobile. In addition, elevated resistance risk factors persisted for some time (up to 22 h) in the receiving water after stormwater events, likely due to inputs from distributed overflows in the catchment. This indicates temporal dynamics should be considered when interpreting the risks of exposure to resistance from event-based contamination. We propose that reducing bypass from wastewater treatment plants may be an important intervention option for reducing dissemination of antibiotic resistance.


Introduction
Antibiotic resistance has been considered one of the biggest challenges to public health, and its global increase has been recognized as an impending public health crisis by intergovernmental entities. Antibiotic resistance genes (ARGs) and bacteria are transmitted not only via direct human-to-human interaction, but also via multi-sectoral routes, such as interdependent routes among humans, food animals, and environments (McEwen and Collignon, 2018). As examples, ARGs and bacteria are enriched in reservoirs such as human patients, hospitals, and livestock, sometimes driven by antibiotic overuse (McEwen and Collignon, 2018). ARGs and bacteria can subsequently be transmitted to environments via sewage, hospital wastewaters, farm runoff, manure, etc. Humans and animals may be exposed to these environmental ARGs and bacteria (McEwen and Collignon, 2018).
Wastewater treatment plants (WWTPs) are a known route through which sewage-borne resistance genes are discharged into the environment (Bürgmann et al., 2018;Rizzo et al., 2013). Even though many studies showed that the level (genes per volume, or absolute abundance) of ARGs decreases during wastewater treatment processes (Ju et al., 2019;Marano et al., 2020;Rodriguez-Mozaz et al., 2015), profound levels of ARGs remain in treated effluents (Lee et al., 2021;Rodriguez-Mozaz et al., 2015). For instance, the levels of common ARGs and bacteria in treated effluents are significantly higher than upstream waters by up to 2 order of magnitude (Lee et al., 2021;Rodriguez-Mozaz et al., 2015;Sabri et al., 2018). Accordingly, the levels of ARGs and bacteria at receiving waters increase after receiving effluents (Lee et al., 2021;Rodriguez-Mozaz et al., 2015), even though the exact degree of increase depends on the proportion of effluent discharge to river discharge (Lee et al., 2021;Ort and Siegrist, 2009).
The abovementioned studies have generally been conducted under conditions of normal WWTP operation. However, many sewer systems or WWTPs will regularly experience high flow events caused by high intensity precipitation. If such flows exceed the capacity of the WWTP or sewer system, it may lead to wastewater bypass (the release of untreated combined sewage into receiving waters) (Weyrauch et al., 2010). During stormwater events, the level of ARGs in wastewater-receiving waters can increase due to the input of wastewater bypass, but high volume flows from precipitation may reduce or compensate this effect. Resistance levels may also increase during stormwater events due to inputs from other non-point sources. Soil runoff from river catchment, and sediment resuspension are considered as potential non-point sources under stormwater events (Tsihrintzis and Hamid, 1997). Storm drains and urban runoff (Baral et al., 2018) may represent temporary point sources at various points in a catchment. However, bypass is expected to contain the highest abundance of resistance determinants (ARGs and resistant bacteria) and multi-resistance risk factors (multi-resistant bacteria and multi-resistance genomic determinants (MGDs)) because untreated sewage contains resistance determinants (Ju et al., 2019;Rodriguez-Mozaz et al., 2015), and multi-resistant bacteria at considerably higher levels than treated effluent (Czekalski et al., 2012).
The potential combined effects of different contributors on receiving waters during stormwater events were investigated in the Stroubles Creek, Virginia, USA by monitoring surface water qualities (Garner et al., 2017). Two other studies monitored combined-sewage overflows and the receiving surface waters noted elevated levels of ARGs in the receiving waters during stormwater events in the Hudson River, Raritan Bay, and Passaic River, USA (Eramo et al., 2017), and North Shore Channel, USA (Chaudhary et al., 2018). In a study performed in the Antelope Creek in Lincoln, Nebraska, USA (Baral et al., 2018) that does not receive wastewater, the authors found that discharges from storm drains had the strongest impact on the riverine resistome (total ARGs identified using shot-gun sequencing) compared to impacts from (potential) sewage leakage, street sweepings, soils, and sediments. While wastewater-receiving river waters have been studied quite intensively so far, there have thus been comparatively few studies systemically examining the combined effects among potential contamination sources. For instance, the contribution of 'bypass' has not been systematically compared with other potential sources (i.e., treated effluent, surface soil runoff, and sediment resuspension), thus whether those other sources also contribute significantly to the increase of riverine ARGs or not remains unanswered. Furthermore, the potential existence of bypass-borne 'multi-resistance' risk factors in the river has not been known.
In this study we investigated the effect of stormwater events leading to bypass discharge on a highly wastewater-impacted river, the River Murg near the WWTP Münchwilen, Switzerland. The main goals of this study were (1) to assess the impact of stormwater-related disturbance and resilience of the resistome of the Murg river by monitoring temporal dynamics of the resistome during stormwater events that lead to combined-sewage bypass, (2) to identify the key source(s) which contribute the most to the increase of riverine resistance levels, and (3) to assess the impact of stormwater events in terms of multi-resistance risk factors. We used various molecular biological, microbiological, and ecological parameters to identify the key factor among different pollution sources. Those parameters include absolute and relative abundance of well-known anthropogenic antibiotic resistance markers (sul1, and intI1) by quantitative PCR (qPCR) and alpha-and betadiversity analysis of the resistome based on environmental shotgun metagenomics. This approach was motivated by the hypothesis that different sources would be distinguishable by the different (relative) abundance and/or the composition and diversity of resistance indicators. To analyze multi-resistance risk factors, we used metagenomeassembled contigs and heterotrophic cultivationpresumptive multiresistant bacteria and MGDs, which are defined as contigs where ≥ 2 ARGs conferring resistance to different antibiotic classes are co-located in this study. The underlying hypothesis was that untreated combinedsewage (bypass) is more likely to contain high levels of multiresistance risk factors, and thus their abundance might increase in receiving waters during stormwater events.

Site description and samplings
Samples were obtained from the River Murg near Münchwilen, Thurgau, in Switzerland which is one of the most wastewater-impacted rivers in Switzerland. This site was well studied under dry-flow condition in 2018 (Lee et al., 2021), and has the following characteristics: (1) Wastewater effluents comprises a high proportion of the downstream river discharge (33.0-38.0% under base-flow condition), (2) Under baseflow conditions, no other known point-source of wastewater inputs exist other than studied WWTP, (3) All potential combined sewer overflows into the river are from the sewer system that connects to the studied WWTP. Other than the sampled bypass which is discharged near the WWTP (Fig. 1-left), there are additional upstream points where combined-sewage may be discharged: The nearest one is located 800 m, and the farthest one 11.0 km upstream of the WWTP (H. Zbinden, personal communication to authors and (unpublished) sewer planning document, June 29, 2021), (4) The WWTP receives combined sewage (community sewage + storm drain) during stormwater events, and the bypass (untreated combined sewage) is discharged after primary sedimentation to the receiving river when the treatment and rainwater storage capacity of the WWTP is exceeded during heavy rainfall, (5) the river catchment is widely utilized as agricultural area (pasture and meadow for livestock farms) (BFS, 2013). Coordinates and key information on the sampling sites as well as information on the WWTP Münchwilen are summarized in Dataset S1.
On June 25 (DRY1) and August 27 (DRY2) in 2019, water samples were obtained under base-flow condition (no precipitation in the previous 24 h) over 24 h in hourly 1 L batches using autosamplers (ISCO, USA) installed at three sampling points: River Murg 200 m upstream (US) and 500 m downstream (DS) of the wastewater discharge point, and from effluent (EF) (prior to discharge). Water samples were cooled with ice and cooling packs inside the autosamplers. Surface soils and sediments (< 5 cm) from the river catchment were obtained in a separate campaign under base-flow condition on May 28, 2019 (DRY0) as shown in Fig. 1-right: Two soil samples each were obtained from forests (F1, 2), meadows (M1, 2), and pastures (P1, 2). River sediments were sampled once at the DS and US locations, and at two sites (M1, 2) further upstream (0.45 and 1.0 km). Both surface soil and sediment samples were obtained by subsampling and pooling 5 subsamples from each location into sterile containers on site.
Event-based sampling campaigns were carried out over the summer of 2019. Weather forecasts (provided by MeteoSwiss) for the Münchwilen region were monitored from July to September 2020. When heavy rainfall was predicted for the catchment area of the Murg, autosamplers were deployed in the same way as described for base-flow sampling for US, DS, and EF sites. An additional autosampler was installed at the wastewater bypass line (BP) of the WWTP prior to discharge ( Fig. 1left). Two event-based sampling campaigns were eventually performed on August 12-13 (13:00-12:00) and September 23 (01:00-24:00), 2019 when there was a total of 10.5 and 20.0 mm of precipitation during 48 h (24 h on the sampling date plus the previous 24 h, refer to Fig. S1) (MeteoSwiss, 2019). During these events, BP could only be sampled when there was sufficient bypass flow while the autosampler was operating. Only 2 samples (at 13:00 Aug. 12, and 12:00 Aug. 13) could be obtained for EF in RAIN1 because of a malfunction of the autosampler between14:00 (Aug. 12)-11:00 (Aug. 13).
All samples were cooled at 4 • C in the dark while being transported to the laboratory within 32 h (from the starting time of auto-samplings).
For RAIN1 and RAIN2 sampling campaigns, equal volumes of hourly water samples (each 1.0 L) were pooled from time intervals before, during, between, and/or after stormwater events for 2-7 consecutive hours as shown in Table S1. Those pooled samples were stored at 4 • C in the dark, and further analysis (cultivation and biomass filtration) was performed the next day (after < 16 h).
The information on BP flow condition and sample pooling was encoded into the sample labeling scheme, which denotes samples as BP or no BP (nBP) depending on BP flow at the time of sampling. The sampling times of pooled samples is also given. For example, "BP-2(+ 8-9 h)" denotes the second sample taken during active bypass flow, and is a composite sample taken after 8 and 9 h from the start. All sample designations are listed in Table S1.

Heterotrophic plate counts
Colony counts of presumptive clarithromycin and tetracycline resistant bacteria (CLR/TET) was shown to be an useful indicator for anthropogenic resistance inputs in our previous studies (Czekalski et al., 2012;Lee et al., 2021). The CLR/TET colonies were cultivated for water samples using biomass concentrated on 0.2 µm pore size cellulose nitrate filters followed by incubation on R2A agar plates in the presence of clarithromycin (4.0 mg/L) and tetracycline (16.0 mg/L) as outlined in our previous publications (Czekalski et al., 2012;Lee et al., 2021). Samples were diluted before filtration using 0.85% NaCl according to the previously optimized ranges (10 mL loading volumes of 10 − 3 -10 − 1 diluted samples) (Lee et al., 2021). Contamination controls using a blank solution (sterile 0.85% NaCl) were performed for each sampling campaign, and no growth of colonies was confirmed. Technical triplicates were incubated for each sample, and standard errors of triplicates are shown as error bars.

Biomass filtration, DNA extraction, and quantitative PCR
To obtain concentrated suspended biomass for DNA extraction, water samples were filtered through a sterilized 0.2 µm pore size cellulose-nitrate filter (Sartorius, Germany) using autoclaved Nalgene™ filter units (Thermo Fisher Scientific, USA). The maximum filtration volume was up to 1.0 L for river waters, 0.5 L for EF, and 0.1 L for bypass, but the exact volume varied by sample because filters were clogged at different volumes (Dataset S2). After biomass filtration, the filters were stored at -20 • C until processing for DNA extraction. Soil and sediment samples were frozen immediately after arrival at the laboratory and stored at -20 • C until processed further for DNA extraction.
DNA extraction was performed using DNeasy PowerWater Kit (Qiagen, Germany) for water samples, and DNeasy PowerMax Soil Kit (Qiagen, Germany) for soils and sediments according to the manufacturer's instruction. After extraction, DNA quality indicators (260/280 and 260/230 absorbance ratios) were checked using NanoDrop One spectrophotometer (Thermo Fisher Scientific, USA), and concentrations were analyzed using both NanoDrop and Qubit™ dsDNA BR Assay Kit (Thermo Fisher Scientific, USA) (Dataset S2). The extracted DNA samples were stored at -20 • C until analyzed. qPCR targeting two resistance indicators (sul1 and intI1) and the 16S rRNA gene were performed as outlined in our previous studies (Czekalski et al., 2014;Ju et al., 2019;Lee et al., 2021). The two resistance indicator genes (sul1 and intI1) used in our study are two of the most widely used genetic markers for tracking anthropogenic sources of resistance (Berendonk et al., 2015;Gillings et al., 2015). The analysis was performed using LightCycler 480 Probes Master kit (Roche, Switzerland) on a LightCycler 480 II (Roche, Switzerland). Primers are given in Table S2, and key information for qPCR validation is given in Table S3. Standards were run in quintuplicates, samples in triplicate. Measurements with Cp values above negative controls but below the limit of detection (LOD; average Cp of lowest valid standard), or which had a standard deviation of Cp values of triplicates > 0.5 were not quantified, and labeled as 'Detected but not quantifiable (D.N.Q.)'. Samples without valid Cp or with average Cp values ≥ smallest Cp value from the non-template controls were indicated as 'Not-detected (N.D.)'. Standard errors of triplicate measurements are displayed as error bars in the figures.
To calculate absolute abundances from qPCR results, we calculated the total number of gene copies per biomass filter and normalized by filtration volumes for water samples (copies per volume). For soils and sediments, the total copies measured in wet mass were normalized to dry mass (copies per g dried soils or sediments). The dry mass of soils and sediments was determined according to Standard Methods (APHA-AW-WA-WPCF, 1981).
The selected DNA samples from our study were sequenced using the Illumina Novaseq6000 with a paired-end (2 × 150) strategy by Novogene Europe (Cambridge, UK). The reads containing adapters and low quality reads (N > 10% and quality score ≤ 5) were removed by Novogene, and the read qualities were double checked by the authors using FastQC v0.11.4 (Andrews, 2010).
For resistome analysis, we used a read-based annotation approach because assembly efficiencies (percentage aligned sequences calculated using Eq. (1)) were relatively low for soil and sediment samples. As a result, resistance gene profiles could not be properly represented using the de novo assembly based annotation approach used in our previous publication (Lee et al., 2021) (Fig. S4).
Where n denotes the total number of assembled contigs in a sample and average coverage indicates the average sequencing depths per contig calculated according to Albertsen et al. (2013).
Read-based ARG annotation was performed using DeepARG shortreads pipeline v1.0.2 with default parameters. This pipeline was developed to process short-reads and to find and annotate resistance genes using a deep learning algorithm after incorporating several ARG databases publicly available (Arango-Argoty et al., 2018). In short, this pipeline finds and quantifies ARG-like and 16S rRNA gene-like reads, and normalizes the reads assigned to ARGs to 16S rRNA gene-like reads (Arango-Argoty et al., 2018). Both our samples and downloaded metagenomes (from SRA) were analyzed using this approach.
In water samples only (where assembly efficiencies were relatively high, > 51.0%), a de novo assembly based approach was used for the limited purpose of analyzing MGDs. We followed the work-flow described in our previous publications (Ju et al., 2019;Lee et al., 2021). A short summary of the bioinformatics work-flow is also suggested in the supporting information (SI). Identification of ARGs was based on two published databases -CARD v3.1.0 protein homolog model (Alcock et al., 2020) for ARGs, and INTEGRALL v1.2 (Moura et al., 2009) for intI1.

Analysis of multi-resistance genomic determinants
In order to list-up MGDs, all contigs containing ≥ 2 ARGs that confer resistance to different classes of antibiotics were sub-selected (Dataset S6). For defining resistance classes, we strictly followed the classification of antibiotics shown in CARD v3.1.0. The occurrence of MGDs was quantitatively assessed in terms of the relative abundance defined as follows in this study (Eq. (2)): Where, Cov MGD indicates the average coverage of a MGD; Cov Contig denotes the average coverage of a contig; m, n indicate the total number of MGDs and contigs, respectively, in a sample.
The ARGs physically co-located within MGDs were further visualized in a directed network using R (igraph). The network analysis was performed using an adjacency matrix produced from two vertices (two colocated ARGs) and an edge (co-occurrence frequency). The edge information was produced according to the following equation (Eq. (3) A detailed example of the derivation process for edges is given in the SI.
We tried to assign taxonomy to MGDs using Kraken2 (Wood et al., 2019), Kaiju v1.7.2 (Menzel et al., 2016) and BLAST against NCBI-nt database (Sayers et al., 2019) as outlined in our previous publication (Lee et al., 2021). In short, we considered the assigned taxonomy as a valid information only if all three approaches yielded a consensus classification at genus level.

Statistical analysis and visualization
All statistics and graphs were produced using R. The Shannon-index was used as a parameter for alpha-diversity, and calculated using Vegan (Oksanen et al., 2013). The dissimilarity of resistome among metagenome samples (beta-diversity analysis) were analyzed using non-metric multi-dimensional scaling (NMDS) in Vegan. Pairwise t-test (p-adjustment method: Benjamini-Hochberg) was performed to test for significant differences among potential contamination sources (wastewaters, surface sediments, soils) in terms of alpha-diversity index.

Absolute abundance of resistance indicators (sul1, intI1, and CLR/ TET) in waters
Two sampling campaigns performed under dry weather (baseflow) conditions were used to establish reference values for resistance indicator levels in US, DS, and EF under normal operation, as shown by horizontal lines in Fig. 2. The results from these campaigns also showed the expected elevated absolute and relative abundance of resistance indicators in EF and the expected general increase of the abundance in DS compared to US as a result of EF discharge. We observed little intraday variability (Figs. S3 and S5).
Analysis of the resistance indicators in water samples from stormwater event-based sampling revealed that BP contained high levels of resistance, even higher than EF (p < 0.01; Kruskal-Wallis test; Fig. S2). Specifically, the levels of sul1, intI1, and CLR/TET in BP were higher than in EF by up to 2.1 orders of magnitude (Fig. 2). Accordingly, the levels of those indicators in up-and downstream water of the receiving river (US and DS) also increased when bypass flow was active. The levels of indicators in both US and DS during bypass events (BP in RAIN1-2; peri-in Fig. S3) were significantly higher than the reference (MOR -NIG in DRY1-2; pre-in Fig. S3) levels (p < 0.01; Kruskal-Wallis and the posthoc Dunn's test) in most cases. One exception occurred for intI1 in US where the difference between pre-and peri-groups was not statistically significant (p = 0.14). It is important to note that we learned from the operators after the campaign that the US sampling point also receives combined-sewage overflows from discharge points further upstream.
The levels of sul1, intI1, and CLR/TET in US and DS were highest during bypass flow events (BP-1 in RAIN1, and BP-3 in RAIN2). The peak levels in bypass-affected DS samples exceeded the DS levels at dry conditions by up to 2.4 orders of magnitude, and even exceeded the levels in EF. This indicates that aquatic riverine resistance at receiving points was profoundly influenced by bypass from WWTP and upstream points of discharge during bypass events. In contrast, levels in EF were not affected or increased only slightly during or after the bypass event (Fig. 2). The levels of CLR/TET in EF at RAIN2 was higher than for baseflow conditions (DRY) (p = 0.0008; Kruskal-Wallis test, n = 8 for each of DRY-EF and RAIN2-EF). However, this was also the case for the EF sample at nBP(+ 1-3 h) and DS levels from this campaign, thus probably reflecting a change in the effluent not related to the bypass event (Fig. 2).
We observed that elevated levels of resistance indicators at US and DS temporarily persisted for some time after BP flow had stopped (Fig. 2). The levels returned to the pre-disturbance status (the levels at Absolute abundance (gene copies or colony forming unit per volume) of resistance indicators in water samples during bypass events: resistance genes (sul1, intI1), and multi-resistant bacterial counts (CLR/TET, presumptive clarithromycin and tetracycline resistant bacteria). (a) Results from the RAIN1 sampling campaign, and (b) RAIN2 sampling campaign. Error bars indicate standard errors for technical triplicates. The averaged abundance during the DRY sampling campaigns (DRY1 + 2) for EF (effluents), DS (downstream), and US (upstream) samples are indicated by colored horizontal lines (dark gray for EF; light gray for DS; blue for US). The averages ± standard errors were plotted in dotted horizontal lines. BP indicates bypass samples, which are only obtained while bypass flow is active during or after precipitation events. NA indicates not available for BP and not analyzed for EF (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article). DRY) only after 17-22 h in RAIN1. The temporal persistence of resistance was even more profound in the campaign RAIN2 when there were heavier rains. The levels of sul1, intI1, and CLR/TET at US and DS were still higher by up to 1.0 order of magnitude compared to the predisturbance status (the levels at DRY) 13 h after the last bypass event (BP-3(+ 10-11 h)).
For most samples, the DS levels of resistance indicators were above the US level (except for intI1 and sul1 levels during RAIN2, nBP(+ 15-17 h) sample) showing a persistent effect of the WWTP (EF and BP). Notably however, the US levels increased as well especially during RAIN2. While the existence of further combined-sewage discharge points upstream may explain this observation, we could not a priori rule out the possibility of other contributions (sediment resuspension and surface soil runoffs). In order to quantitatively assess the potential contributions of different sources, we analyzed relative abundances of resistance indicator genes in the water samples and potential sources (bypass, effluent, sediment and soil) by normalizing to 16S rRNA gene levels.

Relative abundance of sul1 and intI1 in waters, soils, and sediments
Relative abundance (gene copies per 16S rRNA gene copy) of sul1 and intI1 at US and DS increased dramatically during bypass events (BP-1 in RAIN1, and BP-1 to -3 in RAIN2) (Fig. 3). Similar to the pattern that we observed in absolute abundance analysis, the increased relative abundance also temporarily persisted in receiving waters, and gradually decreased over time (Fig. 3). The abundance of sul1 and intI1 was below the LOD (D.N.Q. or N.D.) for most soil samples. Where it could be determined, the relative abundance of these indicators in soil was far below the values for water samples (Fig. 3). The sediment samples showed higher values compared to soils, especially for the sediment obtained from DS. However, all relative abundances determined for sediment were lower than in river water with exception of sul1 in US in some pre-and post-bypass samples (nBP(+ 3-5 h) to nBP(+ 24 h) in RAIN1, and nBP(+ 1-3 h) in RAIN2) (Fig. 3) (M01-43 in Dataset S3). To broaden the dataset for soil, 19 published soil metagenomes were downloaded from NCBI-SRA. We analyzed the ARG content using a readbased approach (see Dataset S3 for key information on raw reads). Alpha-and beta-diversity analysis of resistomes was performed using relative abundance data of ARG identified at resistance subtype-level ANOVA and the post-hoc pairwise t-test of Shannon diversity index values of the retrieved resistomes from the three potential contamination sources in the River Murg catchment (wastewater, sediments, and soils-CH) revealed that the Shannon-index of wastewater resistomes was significantly higher than of soil or sediment (Fig. 4a). The Shannonindex values of soil resistomes from other studies (soils-CN and soils-USA) were similar to or lower than for soils-CH. The two highest Shannon index values in US and DS were observed during and right after bypass discharges (BP-3 and nBP(+ 12-14 h) in RAIN2 for US; BP-3 in RAIN2 and BP-1 in RAIN1 for DS), which were higher than sediments and soils.
Beta-diversity analysis of ARG subtypes using NMDS based on Bray-Curtis distance (stress = 7.49) showed that dissimilarities between resistomes of water and other compartments were profound (Fig. 4b).
Water samples clustered together mostly in the 2nd and 3rd quadrants on the ordination plot, and soils were mostly located in the 1st and 4th quadrants of the plot. The dissimilarities between sediments and water resistomes (especially for US samples) were low compared to those between the soils and water, but still clearly distinct from highly stormwater-influenced river water samples (e.g., BP-3 in RAIN2 for US; BP-1 in RAIN1 and BP-3 in RAIN2 for DS). The resistomes of high stormwater-influenced water samples were however most similar to the bypass samples.

Riverine resistome during stormwater events
Resistance classes with high explanatory power for the resistome beta-diversity (and thus differing between resistomes in different habitats or with different influence of bypass) were selected as follows: (1)  (a)). Statistically significant differences (t-test) among potential contamination sources (wastewaters, sediments, and surface soils) are indicated by red brackets. (b) Non-metric multi-dimensional scaling (NMDS) analysis of resistome structure (resistance gene relative abundance data) in different compartments. BP, EF, US, DS indicates bypass, effluents, upstream, and downstream waters, respectively. Soils-CH are from this study. Soils-CN (Sinkiang/Inner Mongolia, China), Soils-USA (compost unamended; Virginia, USA), and Soils-C (compost amended; Virginia, USA) are from other studies (Zheng et al., 2021;Chen et al., 2019) (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article). resistome profiles were aggregated by class of antibiotic resistance in Dataset S5 and analyzed using NMDS (stress = 5.67) (Fig. S6), (2) the correlation between resistance classes and the ordination, and (3) the resistance classes showing high statistical significance (p-value < 0.0001) of the correlation were selected. A total of 16 resistance types were selected with these criteria, conferring resistance to macrolidelincosamide-streptogramin (MLS), aminoglycoside, bacitracin, betalactam, diaminopyrimidine, fluoroquinolone, fosmidomycin, glycopeptide, multidrug, nitroimidazole, peptide, pleuromutilin, rifamycin, sulfonamide, tetracenomycin-C, and tetracycline.
The relative abundances of those 16 resistance classes of resistance in water, soil, and sediment samples are shown in Fig. 5. Bypass and effluent samples contained in high abundance MLS, aminoglycoside, beta-lactam, diaminopyrimidine, fluoroquinolone, sulfonamide, and tetracycline resistance genes. Those 7 resistance classes showed high relative abundance also in high stormwater-influenced US and DS samples (BP-1 in RAIN1, and BP-1 to -3 in RAIN2). In soils and sediments, the following 7 types of ARGs were prevalent: fosmidomycin, glycopeptide, multidrug, nitroimidazole, pleuromutilin, rifamycin, and tetracenomycin-C resistance genes. In spite of being selected as resistance classes with high explanatory value, bacitracin resistance genes were highly prevalent in waters in general but also occurred in soil and sediment, and the "multidrug" resistance class was abundant in all sample types but especially abundant in the soil and sediment. Both are therefore shown separately in Fig. 5 and were not considered further.
The relative abundance of ARGs in the classes not selected by our criteria are shown in Fig. S7.

Analysis of multi-resistance genomic determinants
The results shown above provide strong evidence that resistance levels increased profoundly in receiving river waters during stormwater events, and that the majority of the increase in resistance determinants in river waters originates from bypass according to the results presented in Sections 3.1-3.3.
Based on these findings, we hypothesized that the risk of exposure to various 'multi-resistance' risk factors would also increase. While this was already partly supported by our CLR/TET analysis in Fig. 2, we aimed to obtain a more comprehensive picture by screening de novo assembled contigs for evidence of multi-resistance. Contigs containing Fig. 5. Resistome profiles by class of antibiotic resistance in river waters (upstream and downstream waters), wastewaters (bypass and effluents), soils and sediments. Total 16 types of resistance classes were selected using biplot analysisthe 14 classes of antibiotic resistance which significantly correlated with the ordination (Fig. S6) (p < 0.0001) are shown in (a); the others (bacitracin and multidrug resistance genes) were displayed separately in (b) because those had much higher values of relative abundance than the previous 14 classes. The order of stacked bars is same as the order of resistance classes in the legend. Unit: Relative abundance per Mille (ARG reads per 1000 16S rRNA gene reads).
≥ 2 ARGs were classified as MGDs, and their relative abundances in each sample were calculated according to Eq. (2). The temporal dynamics of MGDs in terms of relative abundance (Fig. 6) followed a similar pattern to that of other resistance indicators (Figs. 2 and 3). Highest abundances were found in all BP samples. High abundances of MGDs were also observed in the river during BP1 (RAIN1) and BP-1-3 (RAIN2). Elevated relative abundance of MGDs persisted temporarily, but decreased over time (Fig. 6). Notably the relative abundance of MGDs in EF was low compared to BP and to BP-affected river water.
To provide a better understanding of the nature of the resistances contained in the MGDs, we visualized the co-localization of ARGs on MGDs with a network analysis. The directed network shows the frequency with which each gene co-occurs on MGDs with each other gene (Fig. 7a). According to references, many of those genes are expected to be associated with plasmids. For instance, various aac(6 ′ )-Ib, aadA, catB, cmlA, dfrA, sul1, and intI1 homologues were found to be located in bacterial plasmids isolated from wild animals (Dolejska and Papagiannitsis, 2018), activated sludges and treated wastewaters (Tennstedt et al., 2003), sediments, and various water environments (river waters, drinking and wastewaters) (Ma et al., 2017). Among genes associated with animal-, human-and environmental-origin plasmids, the following subtypes also occurred in MGDs of our study: aac(6 ′ )-Ib, aadA, aadA5, catB2, catB3, cmlA5, dfrA14, OXA-2, OXA-10, OXA-129, sul1, and intI1. Assuming that the ARGs that are directly linked to the subtypes listed above are also associated with plasmids, a large portion of the ARGs (35 out of 72 subtypes; Fig. 7b) identified on MGDs are likely plasmid-located. Those potentially plasmid-associated ARGs confer resistance to aminoglycoside, beta-lactam, diaminopyrimidine, fluoroquinolone, phenicol, and sulfonamide antibiotic classes. The ARGs that are previously reported to be associated with gene cassettes of intI1 (Tennstedt et al., 2003) are indicated by larger nodes in Fig. 7b.
The MGDs retrieved from water samples were then subjected to contig-based taxonomy assignment in order to infer potential bacterial hosts (Dataset S7). After annotating taxonomy using two read-based taxonomy assignment tools (Kaiju and Kraken2), a total of 12 out of 126 MGDs showed a consensus at genus level. Those 12 MGDs were additionally subjected to BLAST analysis against NCBI-nt DB. Finally, a total of 5 MGDs showed a consensus among all three approaches (Table 1), and we considered the assigned taxonomy as potential hosts to those MGDs. Four of the five assigned hosts were potential pathogens, including Pseudomonas aeruginosa, Escherichia coli, and Citrobacter sp. On the other hand, Tolumonas auensis was described as an isolate from freshwater sediment capable of toluene production, and is not a known human pathogen (Fischer-Romero et al., 1996).

Wastewater bypass as a major contamination source for resistance during stormwater events
Various lines of evidence (i.e., relative abundance of resistance indicators and alpha-and beta-diversity analysis of various compartments of river systems) revealed that wastewater bypass was a major source of resistance during stormwater events in the River Murg. This indicates WWTPs and sewer systems are important intervention points for tackling dissemination of ARGs in the aquatic environment during stormwater events. Without proper interventions, the risk of public exposure to those genes will remain high.
The risk of exposure to MLS, aminoglycoside, beta-lactam, diaminopyrimidine, fluoroquinolone, sulfonamide, and tetracycline resistance genes through contact with river water is considerably elevated during such events (Fig. 5). Furthermore, bypass events occur frequently in the Murg catchment near Münchwilen. For instance, they occurred on a total of 118 calendar dates in 2019, the year of our study. The bypass events occurred especially frequently during August ~ September, occurring on 54 out of 61 days. The elevated exposure risk also persists for several hours after the stormwater event (discussed in detail below). Considering that beta-lactam antibiotics are among the most commonly prescribed antibiotics in Switzerland for both in-and out-patients (FOPH and FSVO, 2020), the discharge and potential for public exposure to bacteria with those resistance groups could be potentially problematic. While the actual risks of exposure and likelihood of spreading resistance determinants to the population through this route remain unknown, our findings may provide some justification to recommend caution against exposure to river water during and after strong precipitation events in rivers receiving high levels of bypass. Obtaining a quantitative overview of the frequency and magnitude of bypass discharge would be an important asset in this sense.
Considering that both relative and absolute abundances of resistance indicators were relatively stable over time for EF and considerably lower than those for bypass (i.e. incoming combined sewage) reducing combined sewer overflows or bypass by increasing the treatment or Fig. 6. Relative abundance of total multi-resistance genomic determinants (calculated according to Eq. (2)) in water samples. US, EF, DS, and BP indicate upstream, effluent, downstream, and bypass water samples. NA indicates not available for BP and not analyzed for EF, and the asterisk (*) denotes analyzed but not-detected. retention capacities of WWTP and retention basins could be a way to reduce the amount of resistance factors released to receiving waters. Scaling-up existing treatment facilities could be an option, or at least, this aspect could be considered in the early stage of WWTP installation or renovation, i.e., making capacities of new WWTPs be high enough to handle large quantities of incoming combined sewage.

A temporal persistence of bypass-born resistance at receiving waters after stormwater events
Our observation that bypass-borne resistance at US and DS sites persisted for an unexpectedly long time (for 22 h until it fell to predisturbance levels in RAIN1) indicates that the temporal dynamics of resistance deserves further investigation and should be considered when interpreting the fate of event-based resistance inputs in the river and associated risks. The downstream transport of event-based pollutant inputs has been widely studied (Jamieson et al., 2005;Nevers and Boehm, 2010;Parsaie and Haghiabi, 2017), and hydraulic processes (advection, and dispersion) are regarded as main drivers. Our previous study showed that EF is fully mixed cross-sectionally after 500 m downstream distance at the Münchwilen study site (Lee et al., 2021). Thus, the prolonged increase of bypass-borne resistance genes at both DS and US sampling sites appears to be likely due to the advective transport and longitudinal dispersion of upstream inputs from combined sewer overflows. Several such potential discharge points exist in the catchment, however data on their contribution during the observed events is not available.
Given our current data and study design, it is not possible to quantitatively explain the temporal dynamics of resistance genes from bypass and combined-sewer overflows in the receiving waters. We are missing quantitative information e.g., on the discharge, flow duration, and resistance indicator loads in the various sewer overflows connected to the studied river. It is even less possible to predict the situation in other watersheds. The impact of bypass can be expected to depend on additional site-specific properties, such as the proportion of bypass discharges to river discharge, characteristics of incoming combined sewages, river hydraulics, etc. Therefore, further studies will be required to better understand the factors that influence the total loadings, peak levels and duration of event-based inputs of resistance factors from bypass and sewer overflows and to derive larger scale assessments of the The entire networks visualizing total ARGs associated with MGDs, (b) Sub-networks visualizing the ARGs potentially associated with plasmids. The classification of antibiotic resistance is according to CARD v3.1.0. The edges (co-occurrence frequencies calculated according to Eq. (3)) with > 0.5 (50%) are shown as thick red lines, and those with ≤ 0.5 are shown as thin light-red lines. The direction of the arrow indicates co-occurrence of the gene of origin with the target gene. The beta-lactam_1-6 indicate the classes conferring resistance to carbapenem, cephalosporin, penam, monobactam, cephamycin, and penem antibiotics, respectively. The large nodes in (b) are potentially associated with gene cassettes of intI1 (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Table 1
The multi-resistance genomic determinants (MGDs) with successfully assigned genus level taxonomies. Only MGDs with consensus assignment from all three assignment methods used (Kaiju, Kraken2, BLASTn) were shown. P ident indicates the percentage of identification (%), and Q hsp denotes the query coverage for highscoring segment pair (%). contribution of bypass to AMR levels in receiving waters. Catchmentwide measurements of bypass properties (discharge, and representative values of resistance levels in BP, etc.) and inputs during stormwater events combined with hydraulic modeling could provide a clearer picture.

The risk of exposure not only to resistance, but also to 'multiresistance' increases at bypass-receiving waters
To the best of our knowledge, we showed for the first time that WWTP bypass of untreated sewage leads to high relative abundance of potentially problematic multi-resistance factors in bypass-receiving rivers, thus increasing the risk of exposure to multi-resistant bacteria during and shortly after stormwater events.
Untreated sewage contains a resistome much more closely related to the resistome of the human gut compared to effluent (Ju et al., 2019). Thus, the risk of these inputs could be even higher. Many MGDs containing the ARGs shown in Fig. 7b are potentially associated with plasmids, and could be transmitted to previously susceptible cells via horizontal gene transfer (HGT) in the presence of antibiotic-mediated selection. The plasmids containing those MGDs could also persist or evolve in environmental microbial communities, for instance, in the form of attached growths near wastewater-receiving points where continuous anthropogenic disturbance such as EF, and also event-based disturbance (stormwater events) exist.
Recent studies provide evidence that ARGs could persist even in the absence of antibiotic-driven selection. Plasmids are a key to mediate these processes. For instance, under laboratory condition, it has been shown that ARG-encoding non-mobile plasmids persist in E.coli over long timescales under non-selective conditions (Wein et al., 2019). Another study showed that conjugal plasmids encoding ARGs were transmitted to donor cells via HGT at high rates in E.coli even without antibiotic-mediated selection (Lopatkin et al., 2017). The ubiquity of HGT in highly dense population (e.g., biofilms), positive selection coupled with other compensations, and/or population dynamics could be contributing factors (Lopatkin et al., 2017). These findings indicate that there is a risk of environmental persistence of disseminated MGDs even in locations where antibiotic levels are not high, such as downstream river locations, and reservoirs (i.e., lakes).
Contig-based taxonomy assignment results in Table 1 suggest that three MGDs (k121_1708054, k121_852659, and k121_1425685) retrieved from bypass (M13) and bypass-receiving DS water (M14) are hosted by potential opportunistic pathogens (E.coli, and Citrobacter sp.). The risk of potential exposure to multi-resistant E.coli, and Citrobacter sp. in bypass-receiving water under stormwater events therefore has to be considered. Even though these organisms are often found in the human gut and do not normally lead to serious infections in healthy individuals, the risk of infection still exists once exposed to a high dose of pathogenic strains (> infectious dose), such as diarrhea in the case of E.coli (Hunter, 2003), and urinary tract infection for Citrobacter sp. (Abbott, 2011). Even though their relative abundances in our analysis do not appear high (0.0038-0.0065% of contigs in each sample), care should be taken to minimize human exposure to these pathogens. Using our current approach, it was not possible to assign taxonomy to many other MGDs, so a comprehensive overview of MGD-host relationships is not provided. One of the ways to increase success rates for contig-based taxonomy assignment in the future could be to obtain longer contigs so that those sequences could have a higher likelihood to include taxonomic markers (e.g., housekeeping genes). While it would still be difficult to assign taxonomy to highly mobile fragments (i.e., plasmid sequences) that are usually shared by many taxa, it might be possible to assign taxonomy at least to chromosomal associated sequences in this way. Practically speaking, applying long-read sequencing, or combining short-and long-read sequencings could help to obtain longer sequenced or de novo assembled fragments, thus increase the success rate of taxonomic assignment.

Conclusions
• WWTPs are potentially important intervention points for preventing or minimizing discharges of bypass-borne ARGs. Future interventions could be made by increasing the proportion of treated wastewater discharges to incoming combined-sewages during stormwater events. • Temporal persistence of bypass-borne ARGs in the receiving water should be considered when interpreting the fate of aquatic ARGs during and after stormwater events. Transport of upstream combined-sewage inputs could be the reasonfuture study involving hydraulic aspects is required. • The risk of exposure to multi-resistance risk factors increased profoundly in the bypass-receiving river due to bypass inputs during stormwater events. • Large portion of bypass-borne MGDs were expected to be associated with plasmidsproper interventions for tackling discharges of bypass-borne MGDs are required to prevent potential persistence and evolution of those factors in the environment.

Data availability
The raw sequencing data are deposited in NCBI-Sequence Read