Bacterial diversity and health hazards associated with resuspended road matter transported by runoff into storm and combined sewers Open Access

Impermeabilization caused by the construction of roads, sidewalks, buildings, and other infrastructures disrupts the water cycle by increasing runoff (volumes and flow), erosion, and soil leaching. Furthermore, waters flowing over impermeabilized surfaces can lead to the resuspension and dissemination of biological and chemical contaminants that will then be transferred into combined or storm sewers. Runoff among combined sewers will be transferred into a downstream wastewater treatment plant (WWTP), which will reduce their pollutant loads, e.g. Bouchali et al. (2023). However, these combined sewers can saturate during heavy rain events, and their content can sometimes be released directly into streams to avoid street and household overflows (e.g. Pozzi et al. 2024). An alternative in the management of this runoff is thus to implement more local infiltration systems like aquifer recharge systems (e.g. Aigle et al. 2021), road trenches, and rain gardens. Such infiltration systems can involve the use of large stormwater detention and infiltration basins which can handle large volumes of runoff coming from artificialized watersheds. However, the proper functioning of these infiltration systems can be altered by recurrent contaminations which can result in their clogging if not subjected to maintenance routines (Pucher & Langergraber 2019). These aquifer recharge systems were also found to transfer petrol derivatives, nonylphenols, pharmaceuticals, and polar pollutants into their connected natural water compartments (e.g. Wiest et al. 2018; Pinasseau et al. 2020).

Similarly, microbiological contaminants resuspended and spread by runoff can also impact the quality of the receiving downstream systems. So far, these impacts have been mainly documented through the use of fecal indicator bacteria (FIB) such as intestinal Enterococci, fecal coliforms, and E. coli, and quantitative polymerase chain reaction (qPCR) microbial source tracking (MST) assays of fecal bacteria (e.g. Hou et al. 2022). These indicators highlighted the negative impact of runoff on the fecal content of the receiving aquatic systems (Goto & Yan 2011; Crim et al. 2012; Hou et al. 2022), and their relative pressures on these systems were found related to population density, impervious cover (DiDonato et al. 2009; Sidhu et al. 2012) and the type of water management practices (Frenzel & Couvillion 2002). These reports showed the importance of point-source such as discharges from storms (Steele et al. 2018) and combined sewers (Garcia-Armisen & Servais 2007; Hou et al. 2022; Pozzi et al. 2024)), and diffuse pollutions linked to runoff and soil leaching (Garcia-Armisen & Servais 2007), on the spread of such fecal indicator bacteria.

However, microbiological contaminants of these systems coming from runoff are not limited to fecally-borne bacteria and can include several bacterial pathogens and non-pathogenic species from other sources. The impact of these exogenous bacteria on the connected downstream environments remains under study (e.g. Ahmed et al. 2019). Among the few reports on this topic, a study showed that urban runoff bacterial taxa could represent as much as 10% of the bacterial diversity found in aquifer waters fed by an infiltration system, as inferred from a 16S ribosomal ribonucleic acids (rRNA) gene meta-barcoding analysis (Colin et al. 2020). However, the relative contributions of road bacterial taxa over those coming from rainfall on the diversity of microbial taxa found in the investigated downstream compartments had not been defined.

Here, we aimed to further characterizing bacterial community mixing processes over a peri-urban catchment by (i) performing a meta-barcoding deoxyribonucleic acid (DNA) characterization of road bacterial taxa found among runoff communities of a residential peri-urban area with typical European typologies, and (ii) evaluating the contributions of these road bacterial taxa and communities in the buildup of microbial assemblages in the connected combined and storm sewers. The selected area for this study was part of the Yzeron watershed (Lyon conurbation, France) which is one of the main sites of the long-term field observatory dedicated to urban water management (named OTHU; https://www.graie.org/othu/) in France. The hydrology, urbanization, ecology, and microbiology of the Yzeron catchment have been investigated for more than 30 years, facilitating the definition of morphotypes and land use typologies over the catchment. The following hypotheses were tested in this study: (i) DNA meta-barcoding analytical schemes and qPCR assays can be used to differentiate sources of bacterial taxa that can seed road surfaces and lead to the emergence of their core bacterial communities, (ii) DNA meta-barcoding profiles can be used to estimate the relative contributions of these road bacterial taxa in the buildup of bacterial communities in the downstream storm or combined sewers, and (iii) functional inferences from the meta-barcoding taxonomic allocations can identify properties associated with the colonization of road surfaces.

To characterize these road bacterial taxa found in peri-urban runoff, both, 16S rRNA gene and tpm-based meta-barcoding amplicon DNA libraries were generated. The tpm (encoding a thiopurine methyltransferase) gene meta-barcoding approach allows deeper taxonomic allocations down to the species and sub-species levels for more than 60 genera including the Pseudomonas and Aeromonas (Aigle et al. 2021), while segments of the 16S rRNA gene target can allow accurate allocations at the genus level for the Bacteria (and sometimes at the levels of complexes of species through a full sequencing of the 16S rRNA gene). Interestingly, this latter dataset could also be used to infer sources of fecal bacteria that seeded surfaces of a peri-urban zone through bio-informatic tools like FORest ENteric Source IdentifiCation (FORENSIC) (Roguet et al. 2020). These assessments of sources of fecal contaminations were further completed by MST qPCR assays applied to the DNA extracts used to yield the meta-barcoding datasets. These qPCR assays targeted fecal bacteria with host-specific distribution patterns (e.g. Simpson et al. 2002; Marti et al. 2017b; Monteiro et al. 2021). The DNA amplicon gene libraries were further used to infer the relative contributions of road taxa in the buildup of the downstream combined and storm sewage bacterial communities through fast expectation-maximization for microbial source tracking (FEAST) analyses (Shenhav et al. 2019).

Socio-urbanistic surveys were conducted in these areas over three seasons (spring, summer, and fall 2017) within 50 m of the runoff collection area. Each survey corresponded to sixteen counting sessions over day (n = 13) and night (n = 3) sessions of 30 min. This led to the identification of 22 socio-urbanistic variables (technical objects and traces) indicative of peculiar behaviors that may impact the microbial communities colonizing the sampling areas (Table S1b). At sampling time, at least 1 L of water was collected per sampling point. Runoff from G01, G02, G03, and G04 were collected upstream (10 cm) of the drainage grids (Figure 1 and Table S1). Rainfall was collected and termed MW_G06S1 and MW_G06S3 or MW_GaS1 and MW_GaS2 (Table S2). They were collected with 70° sterilized alcohol buckets (dried prior to use) positioned over two distinct open sky areas.

Electrical conductivity, turbidity, pH, oxygen, and water temperature were measured using a multi-parametric sensor probe (Horiba, Piscataway, USA) (Table S2). Total organic carbon, dissolved organic carbon, total suspended matter, total phosphorus, and nitrogen analyses were measured by normalized methods applied by the CARSO laboratory (Lyon, France). Total heterotrophic bacteria were quantified using 1/10 diluted tryptic soy agar (TSA) medium amended with 1.2% agar and inoculated with 100 μL of water sample (1/100 and 1/1,000 dilutions). The number of intestinal enterococci and Escherichia coli were estimated using the IDEXX methods with the Enterolert type E and Colilert kits (IDEXX, Westbrook, USA). The most probable number (MPN) methods were used to detect and quantify Aeromonas caviae and Pseudomonas aeruginosa according to Navratil et al. (2020).

Sampled waters were filtered using a 0.2 μm pore size polycarbonate filter (Merckmillipore, Burlington, USA), and DNA was extracted from the filtered cells using the Qiagen DNA easy Tissue Kit (QUIAGEN, Hilden, Germany). DNA extraction blanks were performed. PCR assays were performed to assess the occurrences of some bacterial pathogens and fecal contaminants in these DNA extracts. These qPCR assays were performed in triplicate using a CFX96 qPCR device (Bio-Rad, Marnes-la-Coquette, France) and the Brilliant II SYBR Green low ROX qPCR mix (Agilent, Vénissieux, France). Total bacterial counts were performed using primers targeting the universal 16S rRNA gene (Marti et al. 2017a). Total Bacteroidales 16S rRNA genes were quantified as described in Layton et al. (2006). Origins of fecal contaminations were inferred from bacterial DNA targets of humans (Seurinck et al. 2005), dogs (Kildare et al. 2007), ruminants (Mieszkin et al. 2010), and bovine (Shanks et al. 2008) intestinal tracts. Sewage-specific bacteria were detected according to Feng & McLellan (2019). Integrons (often encoding ARG) were searched using primers targeting integrase genes of class 1, 2, and 3 according to Marti et al. (2017a). Listeria monocytogenes and Staphylococcus aureus qPCR assays were performed as described in Bouchali et al. (2024). qPCR datasets were normalized by dividing the values by the number of 16S rRNA gene copies estimated per sample.

PCR products of the 16S rRNA gene V5–V7 regions were generated by PCR using primers 799F (barcode + ACCMGGATTAGATACCCKG) and 1193R (CRTCCMCACCTTCCTC) and the PCR temperature cycles described by Bouchali et al. (2022). The tpm gene was amplified using primers ILMN-PTCF2, ILMN-PTCF2m, ILMN-PTCR2, and ILMN-PTCR2m (Aigle et al. 2021) according to Bouchali et al. (2022). Illumina V3 Miseq DNA sequencing of 16S rRNA gene amplicons was performed by MrDNA Services (Shallowwater, Texas, USA). The tpm amplicons were sequenced by Biofidal (Vaulx-en-Velin, France). Raw sequence barcodes and primers were removed from DNA libraries using Trimgalore v0.6.5 software (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). The sequences supplied by MrDNA and Biofidal were in two formats: merged or forward/reverse reads, respectively. 16S rRNA gene sequences (forward and reverse) from the combined sewage were merged using the Mothur package (Schloss et al. 2009). The two datasets (MrDNA and Biofidal) were processed using the dada2 package for R (Callahan et al. 2015). The DADA2 pipeline tutorial (1.16) was used to process the raw 16S rRNA gene and tpm sequences (https://benjjneb.github.io/dada2/tutorial.html). 16S rRNA gene reads were discarded from these datasets if they contained ambiguous bases (maxN = 0) or were less than 340 bp and greater than 400 bp (Table S3). tpm reads were removed from these datasets if they contained ambiguous bases, and were smaller than 200 (forward reads) or 180 bp (reverse reads) (Table S4). Chimeric sequences were removed with the removeBimeraDenovo function of the dada2 package using the consensus method. They accounted for 12.3 and 7.8% of the 16S rRNA gene and tpm sequences respectively. The Decontam v1.16.0 package for R was used to remove contaminating sequences generated from the extraction kit, filters, and over the sequencing process. Contaminant Amplicon sequence variants (ASVs) were identified by increased prevalence in negative controls (with a probability threshold of 0.5) and removed from both datasets if observed (this represented 53 ASVs for the 16S rRNA gene library (Table S5) and 5 for the tpm one (Table S6)). The decontam approach was also used to remove reads from rainfall samples among the runoff (Tables S5 and S6). This led to the removal of 2003 ASVs from the 16S rRNA gene library from the runoff, and 218 ASVs for the tpm library. Overall, 12337 16S rRNA ASVs (Table S7) and 2222 tpm ASVs (Table S8) were conserved for further analyses. 16S rRNA gene reads and tpm reads were taxonomically affiliated using the silva_nr99_v138.1 database (Table S7) and the BACtpm3.0 database (Table S8) with a minimum bootstrap confidence of 80% per allocation. Chord diagrams showing the relative abundance of 16S rRNA or tpm- based taxonomic allocations were made using the ‘circlize’ v0.4.6 package for R.

The FORENSIC platform (Roguet et al. 2020) is a database of 16S rRNA gene V4/V6 regions reference sequences from the Bacteroidales and Clostridiales segregated according to their host-distribution patterns. Bacteroidales ASVs from our peri-urban datasets were retrieved and compared with the FORENSIC database. A total of 535 Bacteroidales ASVs (30380 reads) were recovered. It is to be noted that the investigated runoff amplicon DNA libraries did not show high numbers of Clostridiales ASVs (only 31 over a total of 1,195 reads). As a number of more than 100 ASVs per sample was recommended to run the FORENSIC package, this bacterial group was not further considered. To format the peri-urban 16S rRNA gene sequences, multiple alignments were generated with reference FORENSIC V6 sequences through ClustalX 2.1. Library sequences were trimmed to fit an exact coverage of the sequences from the FORENSIC database. These formatted Bacteroidales sequences were then uploaded and tested on the FORENSIC online platform (https://forensic.sfs.uwm.edu/) to infer best matches with host-specific DNA sequences.

Statistical analyses of the datasets were carried out in R 4.1.0 (2021-05-18, R Core Team, 2021). Spearman rank-order correlations were tested with the corr.test() and Holm adjustment for multiple tests. Spearman correlations were computed using the Hmisc v4.4.0 R package. Samples were analyzed by sample types (storm sewage, rainfall, road surface resuspended matter, or combined sewage) or by sampling zones (resuspended road matter from sites G01 to G04). Diversity indices were computed from the ASV contingency table with the estimate_richness() function of the phyloseq package (McMurdie & Holmes 2013), and differences between indices from the samples were analysed by non-parametric Kruskal–Wallis tests followed by the Dunn test with holm correction (dunnTest() of the FSA v0.9.5 package (Ogle et al. 2021)). Principal component analyses were computed with prcomp() of the factoextra v1.0.7 package on data normalized by the Hellinger method using the decostand() function in the vegan v2.6-4 package (nutriment concentrations assessed as ‘ < 1’ were replaced by 0.5).

ASV distribution pattern analyses were illustrated by non-metric multidimensional scaling (NMDS) ordinations based on Bray–Curtis dissimilarities between samples (ordinate() function of phyloseq package with ‘NMDS’) on relative proportions of taxa (obtained by dividing their respective counts over the total number of reads observed per sample). Significant differences in the dissimilarities between groups of samples were tested by the adonis() function of the vegan package (Oksanen et al. 2024). Permutation tests of multivariate homogeneity of group dispersions were used to determine whether groups varied more (were more spread out) than others using betadisper() of the vegan v2.6-4 package). Indicator ASVs (genera or species) of resuspended runoff road taxa (RM = ASV road runoff libraries being deleted from their rain ASVs) were found using the multipatt() function of the indicspecies package (Cáceres & Legendre 2009) using Hellinger transformed abundance matrices and by comparison with the rainfall, contingency tables using the ‘IndVal.g’ function to correct for group sizes (Tables S7 and S8). The significance of the indicator taxa was based on permutations (n = 9,999). Bacterial community mixings of RM (road-resuspended matter) ASVs into the WN and SN bacterial communities were assessed with FEAST (Shenhav et al. 2019). Analyses were carried out three times and averaged values of the contributions were then used. Differences in the contribution of RM communities to sink WN and SN communities were tested by Wilcoxon signed-rank tests.

The Functional Annotation of Prokaryotic Taxa (FAPROTAX) database was used to make functional inferences from the tpm and 16S rRNA gene contingency tables. Redundant functional groups were removed from the functional contingency table (Tables S12 and S13). The vegan packages (metaMDS() function) were used to construct NMDS ordination on the FAPROTAX functional contingency table (data were transformed into proportions). Bray–Curtis distances were calculated on the relative proportion data. Data significance was assessed with pairwise adonis tests with holm corrections (pairwiseAdonis v0.4.1 package).

Runoff from sub-catchments typical of a small French town was collected for this study. The selected sites were found to be differentiated by human activities (Table S1(a) and (b)) such as the main road G02 site showing the highest number of gasoline engines. The G04 site corresponding to a shopping mall also showed a high number of gasoline engines but was most markedly differentiated from the other sites by higher exposure to various wastes such as plastic bags, sanitary wastes, and cigarette butts. Site G01 was dominated by activities induced by the occurrences of café, restaurants, and a drug store. The G03 site was close to a rural footpath and attracted more joggers and dog walkers than the other sites.

Physico-chemical monitoring of the runoff samples of the above sites did not show large variations (Table S2). The highest pH was 8.1 and the lowest one was 6.8. Turbidity varied from 4.9 to 396 NTU. Only conductivity (52–110 mS/cm on average) and dissolved organic carbon (3.–8.25 on average) tended to increase between the sampling periods. Conversely, the oxido-reduction potential (347–198 mV) tended to decrease between the sampling periods. Sources of exogenous bacteria impacting these peri-urban surfaces were inferred from classical and qPCR assays tracking fecal emitters (using species-specific bacterial DNA signatures) and pathogens. The following bacteria were tracked in the runoff, and their correlations with the monitored physicochemical parameters were assessed (Table S2 and Figure S2): (i) E. coli and intestinal enterococci in order to infer the global impact of fecal contaminations on the quality of the peri-urban runoff and road surfaces, (ii) Bacteroidales respectively specific of human (HF183 target), dog, ruminant/bovine, and sewage, and (iii) P. aeruginosa, A. caviae, L. monocytogenes, S. aureus bacterial pathogens, and integron DNA shuttles, named int1, int2, and int3, contributing at spreading antibiotic resistances among proteobacteria. Overall, total bacterial counts on general culture media, and those inferred from 16S rRNA gene copies, showed about 4 log differences (about 6 × 10⁶ CFU per 100 mL or 1 × 10¹⁰ 16S rRNA gene copies per 100 mL on average). These counts were found higher in the combined sewage (about 6 × 10⁸ CFU per 100 mL or 2 × 10¹¹ 16S rRNA gene copies per 100 mL; KW test with post-hoc Dunn tests with Holm's correction p < 0.05), and were lower but still detectable in the rainfall (Table S2). E. coli and intestinal Enterocci represented about 0.3–0.4% of the culturable bacterial cells in the road runoff (without significant differences between sites), and about 0.5–2.25% of those of the combined sewage (significantly higher E. coli counts in the sewer than the runoff; KW test with posthoc Dunn test with Holm's correction p < 0.05), and suggested significant fecal contamination of runoff flowing over the investigated peri-urban town. MST qPCR assays were carried out to identify some of the fecal emitters involved in these contaminations (Table S2). Runoff from three sites (G01, G03, and G04) showed recurrent contaminations by Bacteroides from human, dog, sewage, and ruminant fecal matter (Figure 1). G02 did not show the HF183 DNA targets for human Bacteroides but showed significant amounts of sewage-specific DNA targets. This site is located on the main road going through the town. Bovine Bacteroidales were detected at G01 and G02, and ruminant Bacteroidales were detected at all sites (Table S2). The significant occurrences of fecal bacteria at the sampling points were indicative of probable contamination by bacterial pathogens. This was confirmed by a significant detection of P. aeruginosa cells (from 7 to 10³ MPN cells/100 mL) and L. monocytogenes DNA targets (1 × 10³–7 × 10⁵ gene copies/100 mL) among all runoff DNA samples. A. caviae (7–2 × 10³ MPN/100 mL), and S. aureus (1 × 10² to 7 × 10³ gene copies/100 mL) were also detected but not among all runoff (Figure 1 and Table S2). S. aureus and L. monocytogenes DNA targets were found at similar levels between the runoff and the sewage. Integrons 1, 2, and 3 DNA targets were also detected from runoff collected at most sampling sites (Table S2). Counts of ruminant Bacteroidales, P. aeruginosa, and A. caviae were correlated (Spearman correlation tests; Figure S2) to turbidity, suspended matter, and total nitrogen concentrations. L. monocytogenes DNA counts were correlated with turbidity and the suspended matter, and the bovine Bacteroidales DNA counts were correlated to total-P concentrations (Figure S2). The relative counts of ruminant Bacteroidales gene copies, P. aeruginosa, and L. monocytogenes were found correlated. This was also the case for specific sewage Bacteroides, total Bacteroidales, and A. caviae numbers. Dog Bacteroidales relative counts were also correlated to the culturable E. coli and S. aureus gene copy numbers (Figure S2). These distribution patterns clearly indicated multiple sources of contaminations in the investigated area with a mixing of bacterial communities coming from multiple animals but also humans and their associated wastes. Further characterizations of these runoff microbial communities were performed below through large-scale DNA meta-barcoding analyses.

DNA meta-barcoding datasets of the 16S rRNA gene and tpm PCR amplicons were generated for this study (Tables S3 and S4). These libraries included runoff DNA reads, and reads from the combined (WN) and storm (SN) sewers, and rainfall (MW). ASVs shared with the blank samples were retrieved from these libraries. FORENSIC analysis of the 535 Bacteroidales ASVs (30,380 reads) found in these generated 16S rRNA gene libraries was undertaken to confirm some of the qPCR MST assays performed on runoff DNA samples (Table 1). Among these Bacteroidales ASVs, sewage ASVs could be detected by FORENSIC with high significance (Table 1) among the combined sewage, thus confirming the qPCR datasets. High-confidence sewage-like ASVs were then detected among the runoff bacterial 16S rRNA gene libraries collected from the sampled peri-urban roads and from the storm sewer. The storm sewer samples were the most contaminated, and the G01 town center and G02 road runoff contributed significantly to these contaminations. Nevertheless, the highest fecal contaminations recorded among these road runoff were related to ruminant fecal taxa, confirming the qPCR assays presented in Figure 1 and Table S2. These analyses also indicated a probable contamination of rainfall but at a low confidence because of the low number of Bacteroidales 16S rRNA gene reads and ASVs in these samples. Still, the qPCR MST assays for human and dog contaminations also indicated a low but significant fecal contamination of these rain waters (Table S2).

After having performed these FORENSIC analyses, ASVs from rainfall were retrieved from the runoff libraries to generate the road-resuspended matter (RM) amplicon libraries. This was done by using the decontam approach on both the 16Sr RNA gene (Table S5) and tpm libraries (Table S6). This resulted in a final dataset of 2,125,891 16S rRNA gene reads and 999,641 tpm reads (Tables S3 and S4) grouped into 12,337 16S rRNA ASVs and 2,222 tpm ASVs (Tables S7 and S8). Diversity indices were computed from these ASV contingency tables (Figure S3, and Tables S3 and S4), and significant variations were observed between the numbers of road surface ASVs (termed road-resuspended matter or RM) and ASVs found in the combined sewage (termed WN). The 16S rRNA-based ASV richness in terms of ASV numbers was significantly lower in the combined sewage (1,505 ASVs) and was higher, but not significantly, in the storm sewage than in the RM of the surface runoff (6,630 vs. 5,024 16S rRNA ASVs) (KW test with post-hoc Dunn test with Holm's correction p < 0.05) (Table S3). This was further supported by a Shannon index of 6.74 computed for the 16S rRNA gene road RM dataset while it was 5.36 for the combined sewage and 686 for samples coming from the storm sewage (KW test with post-hoc Dunn tests with Holm's correction, p < 0.05) (Table S3). Similarly, using the tpm dataset, richness was found highest among the storm sewage samples (index of 6.34) but found relatively similar to the RM road samples (index of 3.62) and those from the combined sewage (Table S4). These trends were confirmed with the Shannon index.

A total of 45 bacterial phyla (5 with a relative abundance greater than 1% in the raw contingency of Table S9) were inferred from the 16S rRNA V5–V7 gene libraries generated in this work. Combined sewage samples showed the lowest number of phyla (n = 16), followed by rainfall and resuspended street surface taxa of the runoff (both n = 33). As expected from the diversity indices (Figure S3), samples from the storm sewage (SN) had the highest number of phyla (n = 42) (Table S9). Combined sewage samples were highest in Firmicutes relative abundances (mean relative abundance per sample of 12%), while the other samples were dominated by Actinobacteriota (9.4 and 5.1% mean relative abundances in road SM and the storm sewage samples). About 18 genera were found to account for more than 1% of the 16S rRNA gene reads (Table S10).

Storm sewage (n = 786) showed the highest number of genera, followed by the road RM samples (n = 704), and those from the combined sewage (n = 364). Flavobacterium and Massilia dominated among the RM and SM taxa. High relative numbers of Sphingomonas and Spirosoma in RM (6.0 and 5.6% mean relative abundance), and of Malikia and Duganella in the storm sewage (4.7 and 4.1% mean relative abundances) were observed (Table S10). In the combined sewage, the Aeromonas, Bacteroides, Prevotella, and Acinetobacter dominated (27, 11, 5.7, and 5.6% mean relative abundances per sample, respectively). Regarding the tpm taxonomic allocations, 15 genera and 123 bacterial species could be inferred (Table S11) (19 with a mean relative abundance over 1%). In road RM, tpm sequences allocated to a Pseudomonadales bacterium coded RIFCSPHIGHO2 and affiliated to Pseudomonas peli by molecular phylogeny, and sequences allocated to Pseudomonas syringae, and Pseudomonas lurida were dominant (23, 11, and 11% mean relative abundances, respectively). In storm sewage, the P. peli DNA reads were also found to dominate with reads from Pseudomonas sp. GV071 and Pseudomonas oleovorans (with 31, 15, and 8% mean relative abundances, respectively). Aeromonas media, Pseudomonas sp. SL4 and A. caviae tpm sequences dominated among the combined sewage tpm libraries (30, 27, and 26% mean relative abundances, respectively).

Bacterial genera and species contingency tables inferred from the 16S rRNA and tpm ASV tables revealed no significant differences in the Bray–Curtis dissimilarities between the SN and road RM profiles, but significant differences were confirmed with those computed from the combined sewage (WN) dataset (pairwise Adonis; p < 0.05; Figures 2(b) and S4(b)). These taxonomic contingency tables per sample were transformed by FAPROTAX into functional contingency tables (Tables S12 and S13), after a cleaning step to eliminate redundant functional groups. These functional groups were used as extrinsic variables, and projected by envfit onto the genus- or species-based NMDS (Figures 2(b) and S4(b)). Several functional groups were found to be strongly correlated with the inferred NMDS ordinations of the samples (Table S14). Genera and tpm-inferred species from combined sewage (WN) were differentiated by their ability to reduce nitrate. Fifteen and 23% of the genera and species had, respectively, this capacity, on average, compared to only 0.83 and 0.79%, and 4.4 and 8.3% for the road RM and SN inferred taxa. WN ordinations also matched the relations with the occurrences of inferred taxa related to animal parasites or symbionts (4.2% of genera on average against 0.85 and 0.45% in the RM and SN). WN and SN were further found to harbor acetogens likely involved in the transformation of carbon dioxide (CO₂) and hydrogen (H₂) into acetic acid (acetate) under anaerobic conditions (acetogenesis) (Figure 2(b)). The relative abundances of methylotrophes, or taxa involved in sulfur respiration, chlorate reduction, and nitrate reduction were found higher in the storm sewage (SN) than road RM according to the genus-based NMDS ordinations (Figure 2(b)). Similar numbers of inferred bacterial plant pathogens (about 2%), and tpm-species involved in cellulolysis or lignolysis (Figure S4B), but higher hydrocarbon degraders (according to the 16S rRNA-inferred genera) were found in the road RM than SN samples (0.13% on average rather than 0.03% in the SN) (Figure 2(b) and Table S12).

Distribution biases of inferred genera (from the 16S rRNA gene target) or species (from the tpm target) were investigated through the indicspecies bioinformatics package. This approach allows to assess of the positive association of a taxa to a particular category of samples through specificity and sensitivity tests (Cáceres & Legendre 2009). This led to the identification of 20 genera and one species, P. syringae, as being potential indicators of road RM while comparing the dataset with the ones from rainfall (Table 2). These bioindicators had high specificity (highly associated with RM and mostly not found in rainfall) and high sensitivity (number of times found in the samples of a category). Five genera were additionally found as being indicative of particular road sites: (i) Seohaeicola reads were associated with G01 in the town center, (ii) Neorickettsia with G03 in a family housing estate, (iii) three genera, Defluviitaleaceae, Fibrobacter and Polyangium which were specific of the G02 site located on the main road going through the town. No specific indicator taxa for G04 (parking of the shopping mall) were identified. All these indicator taxa except one (Neorickettsia) were recorded in the storm sewage (SN). Two indicator taxa (Actinotalea and Fibrobacter) were recorded in the combined sewage but at a lower relative abundance than in the RN samples (Table 2).

Several studies investigated chemical contamination of runoff and street sediments and demonstrated occurrences of nutrients (total-nitrogen (N), total-phosphorus (P)), and enrichment in pollutants like metallic trace elements (Cu, Sb, Pb, and Zn) and polycyclic aromatic hydrocarbons when these matrices are exposed to gasoline engines (e.g. Revitt et al. 2014). These nutrients and pollutants are sources of energy, carbon, and nitrogen for micro-organisms, and likely contribute to the buildup of original microbial communities among street sediments. However, the microbiology of road microbial communities, and the knowledge of taxa getting resuspended in runoff, remain poorly documented (Ahmed et al. 2019). The nature of these road taxa needs to be explored in terms of sources feeding these microbial communities but also of functional units showing the ability to thrive on pollutants under these particularly extreme conditions (desiccation, high temperature in the summertime). Furthermore, the ability of these taxa to get in runoff and spread among the downstream networks like storm and combined sewers needs to be addressed to evaluate their relative impacts on these microbial assemblages. Here, road runoff from a peri-urban town was thus sampled over three seasons, and their microbiological contents were investigated by classical and DNA-based approaches. The selected areas were typical of European small-town typologies with occurrences of café, restaurants, and drug stores, and the proximity of shopping malls. The road runoffs of these areas showed similar quantities of total suspended solids, total-N, and total-P as those reported for other low-density European contexts (e.g. Lundy et al. 2012).

However, fecal contaminations were found to have significantly contributed to the buildup of the microbial communities of the sampled areas. This confirmed previous reports such as those of Lundy et al. (2012) and Ahmed et al. (2019). E. coli and intestinal Enterocci represented about 0.3–0.4% of the culturable bacterial cells found in the runoff collected from these peri-urban zones. Total counts of intestinal Enterocci and E. coli were found slightly higher by about 1 log than those reported in runoff from low urban density areas of Brisbane (Australia) (Sidhu et al. 2012) but similar to those reported for Lisbon city (Monteiro et al. 2021). Interestingly, as reported in these latter studies, these occurrences of FIB were associated with significant qPCR counts of DNA markers indicative of human fecal contaminations, and of mobile genetic elements carrying antibiotic-resistance genes (ARG). Integrons 1, 2, and 3 DNA markers were detected among all sampled runoff as previously reported for other urban and peri-urban contexts (e.g. Hou et al. 2022). Furthermore, as observed in Lisbon, runoff samples from the French peri-urban area were also significantly contaminated by dog fecal matter. Additional fecal emitters, not always tracked in the other studies such as ruminants further confirmed the high impact of animal fecal sources in the buildup of these road microbial communities. Some of these patterns were confirmed by FORENSIC bio-informatic analyses of the ASV reads of the 16S rRNA gene libraries. High scores of ruminant-specific Bacteroidales DNA reads were found by this approach in the road-resuspended matter (RM) 16S rRNA gene libraries. Combined and storm sewage libraries also showed by FORENSIC significant scores of sewage Bacteroidales DNA signatures, and the RM libraries were confirmed to harbor sewage-related taxa. These fecal contaminations were further found associated with bacterial pathogens such as P. aeruginosa, A. caviae, L. monocytogenes, and S. aureus. Dog Bacteroidales relative counts were found to correlate to culturable E. coli counts and S. aureus gene copy numbers. Interestingly, the MPN PCR counts of P. aeruginosa, and qPCR counts of L. monocytogenes, and of ruminant Bacteroidales indicator DNA, were found correlated to the runoff turbidity values and loads of suspended matter. These correlations supported the hypothesis of remobilization of these bacterial species by runoff shear forces operating over street surfaces. It is worth noting that, respectively, DOC and total-P values below 7 and 0.07 mg/L were previously shown to lead to a decay of FIB in surface waters (Surbeck et al. 2010). Total-P values above 0.07 mg/L were not recorded for the investigated runoff but about 50% had DOC above 7 mg/L indicative of conditions that can sustain microbial growth.

To go further into the identification of road bacterial taxa found among runoff, DNA meta-barcoding approaches based on the bacterial 16S rRNA and tpm genes were implemented. These meta-barcoding DNA libraries were submitted to a comparative analysis with the ASVs generated from rainfall (MW) to remove these latter ASVs from the runoff libraries and generate the road-resuspended matter libraries (RM). The genetic diversities of these RM gene libraries were then compared with those of the downstream storm (SN) and combined (WN) sewers. The 16S rRNA gene-based diversity indices showed a lower richness in the combined sewage (1,505 ASVs) than those of the storm sewage and RM datasets (6,630 vs. 5,024 16S rRNA gene ASVs). Similar trends were observed with the tpm gene libraries. This can be explained by strong selective forces operating in a sewer, and directing the evolution of microbial communities toward enrichments of highly specialized taxa such as sulfate-reducers and hydrogen-producing acetogens. These divergences between the WN, and RM/SN compartments in terms of bacterial genetic community profilings were confirmed by NMDS ordinations of the dissimilarities between the ASV contingency tables of the samples.

Taxonomic assignation of the 16S rRNA gene-based ASVs was performed, and core bacterial phyla were found in similar relative numbers between the investigated compartments except for the Firmicutes which were higher in the combined sewage. The high content of Firmicutes building up among combined sewage was previously reported by other authors (e.g. Mathioudakis & Aivasidis 2009; Marti et al. 2017b). Most significant differences in bacterial taxa between the compartments were recorded at the genus level. Road RM was found dominated by DNA reads from the Flavobacterium, Massilia, Sphingomonas, Pseudomonas, and Spirosoma, while combined sewage (WN) was dominated by DNA reads from the Aeromonas, Bacteroides, Prevotella, Acinetobacter, and Flavobacterium. Pseudomonas DNA reads represented between 2 and 5% of the RM, WN and SN bacterial communities. Interestingly, Flavobacterium ASVs were also recorded among storm sewage (SN) but their DNA read numbers were found to be higher than in the other compartments. These trends confirmed the strong association of Flavobacterium with natural water systems (Marti et al. 2017b; Wu et al. 2024). Among SN, these flavobacteria were found associated with other specialized taxa like Malikia (found in activated sludge), Duganella (leaf colonizers), Rheinheimera (broad water inhabitants), and Limnohabitans (found in freshwater systems), with brief descriptions available at https://lpsn.dsmz.de/genus. In other studies, Flavobacterium relative read numbers in combined sewage were also found to be quite high, with a relative proportion reaching about 15%, and other freshwater taxa like Limnohabitans were also found to be significant (Marti et al. 2017b). Here, this latter genus represented only about 0.3% of the WN dataset but showed an occurrence at about 3% in the SN samples. This suggests that strong rain events had probably impacted the WN bacterial community investigated by Marti et al. (2017b). This is further supported by the absence of taxa like Acidovorax (soil inhabitants and plant pathogens), and Aquiflexum (isolated from various water systems) in our study while these were at significant levels in the latter study. Nevertheless, as reported by Marti et al. (2017b), significant DNA read counts of Aeromonas (about 27% of the dataset), Pseudomonas (about 2.2%), and Acinetobacter (about 5.6%) were recorded in WN. A deeper look at the classification of these taxa at the species level (using the tpm datasets) showed A. caviae and A. media to represent about 60% of the Aeromonas total read counts of the WN samples, but only 0.1% of those of the resuspended road matter (RM). This is in line with a recent report that showed A. caviae to be a reliable indicator species of spillovers from combined sewers (Pozzi et al. 2024). Distinct Pseudomonas species were recorded between the RM, WN, and SN samples. RM was dominated by P. syringae, a phytopathogen, while SN and WN showed higher numbers of P. oleovorans known to be able to oxidize C5–C12 n-alkanes (Van Beilen et al. 1994). It is to be noted that P. syringae can colonize multiple plants going from legumes like peas, and crucifers to cherry trees. The occurrences of this species can thus be explained by an accumulation of infected plant debris on road surfaces.

To go forward with the understanding of the bacterial distribution biases observed in the study, the Indicspecies and FAPROTAX computer packages were applied to the datasets. Indicspecies confirmed P. syringae as being a reliable indicator of road-resuspended matter. Additionally, this tool identified 20 genera showing an association with road matter such as Bartonella and Cryptosporangium (a filamentous actinobacteria). The observation of Bartonella DNA reads in RM could be explained by their association with fleas, ticks, and lice among animals that can be seen in peri-urban contexts (Taber et al. 2022). In fact, the MST datasets had shown recurrent contamination of the investigated roads by dog fecal matter. Bartonella-infected flea fecal matter could thus have been released by dogs (and potentially cats, but this was not analysed in this study), and got accumulated in the road deposits getting resuspended in the analysed runoff. The ‘functional groups’ contingency table inferred from the taxonomic classifications further confirmed such associations of intracellular pathogens with road-resuspended matter (RM). Bacterial taxa involved in hydrocarbon degradation were also found to be in higher numbers in the RM and SN compartments than the WN one. This is obviously related to hydrocarbon contaminations from gasoline engines polluting the road sediments and deposits. Interestingly, genera and species from the combined sewer (WN) showed a greater number of taxa able to reduce nitrate. Previous reports have shown that denitrification can take place in sewer systems (Mathioudakis & Aivasidis 2009). WN and SN were also found to harbor acetogens likely involved in the transformation of CO₂ and hydrogen (H₂) into acetic acid (acetate) under anaerobic conditions (acetogenesis). Acetogenesis is a key step in the anaerobic digestion process of organic matter (e.g. Abdelgadir et al. 2014) and was reported to take place in sewers (Jin et al. 2015). These functional inferences and differentiations confirmed species reassortments between compartments. Road surface matter favored the growth of aerobic heterotrophs, and hydrocarbon degraders, while storm and combined sewage favored taxa involved in anaerobic digestion processes of organic matter.

These functional specializations between compartments questioned the ability of road surface bacterial taxa to get established among the connected downstream storm and combined sewers. To define the relative ability of road surface bacterial taxa to colonize the SN and WN compartments, FEAST analyses were undertaken. This confirmed the clear differentiation of the RM and WN bacterial communities. Less than 3% of the RM ASVs were found among WN confirming a clear functional switch toward anaerobic organic matter transformation processes in this compartment. However, up to 53% of the SN bacterial communities were found related to the ASVs detected among RM. This proximity between these communities was previously recorded through Bray–Curtis dissimilarity analyses between the ASV contingency tables of both the 16S rRNA and tpm gene libraries. The RM and SN 16S rRNA gene profiles were found to be closely related and apart from those of the WN sewage. The main differences between the RM and SN communities appeared to be related to sulfate/sulfite respiration which appeared to be harbored by more taxa in the SN bacterial community. This was likely related to the growing conditions in the storm sewer favoring anaerobic digestion processes of sulfate and acetogenesis.

The hygienic quality and microbial diversity of debris and deposits leading to the formation of sediments accumulating over road surfaces have been rarely assessed. Here, significant and recurrent bacterial taxa accumulating and growing in these sediments, and getting resuspended by runoff shear forces, were identified over a peri-urban catchment through qPCR assays and DNA meta-barcoding analyses. Road surface microbial communities were clearly found to be seeded by fecal bacteria from various hosts such as humans, ruminants, and dogs. These fecal contaminations were further found associated with the occurrences of human pathogens like Aeromonas caviae (intestinal tract infections) and Bartonella spp (vector-borne infections transmitted by fleas, and characterized by crusted erythematous papules). These road matter bacterial assemblages were also found enriched in pollutant degraders with high counts of Sphingomonas and Alkanindiges DNA reads. These assemblages clearly showed signs of functional adaptation to the growing conditions prevailing in road sediments. Low numbers of these road surface taxa (less than 3% of the ASVs) were found to be competitive in the downstream combined sewage, but strong proximity (up to 53%) was observed with the bacterial taxa of storm sewage (SN). However, the SN communities showed signs of a functional drift toward anaerobic sulfate-reducing activities.

This study underscores the importance of assessing the hygienic quality of road surface sediments. The effectiveness of road cleaning practices – such as sweeping followed by washing, push broom cleaning, and high-pressure water cleaning (Van Kampen et al. 2020) – should be evaluated in terms of their ability to reduce potential human exposure not only to fecal pathogens, but also to vector-borne agents transmitted by fleas, ticks, and lice. Public awareness should be increased regarding the influence of individual hygienic behaviors on urban surface cleanliness and the transmission of infectious agents. Our findings demonstrate substantial microbiological contamination in runoff, reinforcing the need for thorough quality assessments before any reuse for gardening, cleaning, or similar purposes.

We thank the OTHU (funded, in part, by the French water agency Rhône-Méditerranée- Corse, and the Greater Lyon Urban Community) for its technical support. We also thank the School of Integrated Watershed Sciences H2O'LYON (ANR 17-EURE-0018), and the Lyon Urban School (LUS) (ANR-17-CONV-0004) for their support in the elaboration of this research initiative.

This work was funded by the French national research program for environmental and occupational health of ANSES under the terms of project ‘Iouqmer’ EST 2016/1/120, the MITI CNRS project named Urbamic, the Greater Lyon Urban Community through the Lyon OTHU urban drainage field observatory of ZABR (Zone atelier basin du Rhône), and ANR projects Chypster (ANR-21-CE34-0013) and CityBioHazards (ANR-22-CE34-0016-01). Part of this work was also performed under the terms of the CNRS International Research Programme (IRP 2021–2025) with India on the topic of river-aquifer exchanges. A. Dominguez-Lage was awarded a 50% PhD fellowship from H2O'LYON (ANR 17-EURE-0018).

A.D.L. conceptualized the work; rendered support in data curation; formal analysis; investigated the project; developed the methodology; validated the work; visualized the study; wrote the original draft; wrote and reviewed and edited the article. A.M.P., R.B., B.Y., C.M., J.Y.T., and W.G. rendered support in data curation; and investigated the project. B.C. contributed in data acquisition; rendered support in data curation and formal analysis; developed the methodology, validated the study, conceptualized the work; rendered support in funding acquisition; support in project administration; contributed in resources; supervised the work; validated the work; wrote the original draft; wrote and reviewed and edited the article.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.