Six sampling sites were chosen around the Houston Clear Lake area (Figure 1) primarily based upon accessibility after initial rainfall and the subsequent receding of floodwaters. Sampling began as soon as sites were accessible and continued every 1–2 weeks for 2 months. Water samples were collected via a grab-sampling technique and stored inside sterile, plastic bottles (IDEXX, Westbrook, ME) and transported, on ice, to the lab for further analysis within 6 h of collection as described previously (LaMontagne et al., 2022).

Relative quantities of total Bacteroidales and human Bacteroidales were measured via quantitative polymerase chain reaction (qPCR). Approximately 50–100 ml of water sample was filtered onto 0.22 μM Supor membrane filters of 47 mm diameter and stored at-80°C. DNA was extracted from the filters using the DNeasy PowerWater Kit (Qiagen, Hilden, Germany) and quantified using a nanodrop spectrophotometer (Thermo Scientific, Waltham, MA). The qPCR reaction conditions and thermal cycler programs used for both Total Bacteroidales and Human Bacteroidales assays are listed in Supplementary material (Bernhard and Field, 2000; Shanks et al., 2009; Shanks et al., 2016). For the total Bacteroidales assay, a gBlock standard based on the 16S rRNA gene of Bacteroides fragilis was developed (sequence in Supplementary Table S4) and used as standard for the qPCR assay (Integrated DNA Technologies, Newark, NJ). This standard was diluted to create a standard curve between concentrations ranging from 1.32 × 10⁸ to 1.32 × 10¹ copies ml⁻¹. For the Human Bacteroidales assay, a gBlock (Shanks et al., 2016; sequence in Supplementary material) diluted between 3.2 × 10⁸–3.2 × 10¹ copies ml⁻¹ was used to create a standard curve. Eppendorf LoBind® Tubes were used for all assays to ensure maximal recovery of nucleic acid molecules (Eppendorf, Framingham, MA). Gene copies were calculated and reported as copies per 100 ml of water filtered.

Absolute quantification of the ARGs intI1, aint1, and sul1 and the plasmid pBI143 was performed using an Applied Biosystems QuantStudio™ Absolute Q™ Digital PCR System per the manufacturer’s recommendations. Purified DNA samples, which were archived previously (LaMontagne et al., 2022, #4256), were analyzed for all four target genes in a single 4-plex reaction. Each 10 μl-PCR mixture consisted of 2.0 μl of 5X Absolute Q™ DNA Digital PCR Master Mix, 1 μM of each primer, 80 nM of each probe, 1 μl of DNA template, and a final concentration of 0.2 mg/ml of BSA. After centrifugation at 10,000 X g for 1 min, 9 μl of sample reactions were loaded into individual wells of an Absolute Q™ MAP16 Digital PCR plate and covered by 15 μl of isolation buffer. All multiplex reactions were run with the following thermocycler conditions: 10 min of initial denaturation at 98°C, followed by 45 cycles of 95°C for 30 s and 60°C for 1 min. Data analysis was performed on-instrument using provided Absolute Q™ Digital PCR Software (version 6.3.0-RC11). Primers and probes, described previously (Heuer and Smalla, 2007; Barraud et al., 2010; Quintela-Baluja et al., 2021; Fogarty et al., 2024) are shown in Supplementary Table S1. All labeled PCR probes (6,000 pmoles each) were purchased from Thermo Fisher Scientific. All primers and synthetic gBLOCK positive control sequences were purchased from IDT. Gene copies were calculated and reported as copies per 100 ml of water filtered.

As an additional source tracking approach, we further analyzed 16S rRNA sequence data that was previously published for these same water samples. The samples were initially sequenced using a MiSeq Illumina platform as described by LaMontagne et al. (2022) and results deposited in the NCBI SRA under accession number/Bioproject ID: PRJNA795782. We used SourceTracker2 to reanalyze these data with a focus on human source contributions. Raw sequencing reads were collected and analyzed with the QIIME 2 software package version 2022.8 (https://qiime2.org; Bolyen et al., 2019) and plugins associated with this version. After importing raw reads to the Qiime2, paired-end sequencing reads were demultiplexed using the demux plugin.¹ Quality control, filtering chimeric sequences and feature table construction were done using the q2-dada2 plugin (Callahan et al., 2016) with trimming parameters based on the demux output visualization files. Taxonomy classification was done against the SILVA 138 database (Quast et al., 2012; Yilmaz et al., 2014) using the “feature-classifier” plugin (Bokulich et al., 2018) with the “classify-sklearn.” Taxonomy-based filtering was conducted to remove mitochondria, chloroplast, Archaea and Eukaryota using “taxa-filter-table and taxa-filter-seqs” command from both taxonomic and rep-seqs tables. A taxonomic bar-plot was generated using the filtered table and assigned taxonomy file and visualize through the QIIME 2 view website. The resulting filtered taxonomy file was used to extract sample specific taxonomic classification. SourceTracker2 software with QIIME Gibbs Sampler plugin was used to determine source-sink proportions and taxonomic classification based on ASV abundance (Knights et al., 2011).

Multivariate regression analysis was conducted on previously published (LaMontagne et al., 2022) water quality data [dissolved inorganic nutrients (DIN), salinity, temperature, oxygen and pH] and the metagenomic and dPCR data produced as described above. Values for missing a set of DIN measurements, which were not collected at station H (Figure 1) immediately before Hurricane Harvey made landfall, were imputed from DIN values pulled from a public archive water quality measurements for that water body² and processed with scripts presented in Supplementary file S2. The R (R Core Team, 2021) package mice (v 3.16.0) was used to impute the missing values with scripts presented in Supplementary file S3. The input Supplementary file S4 in this analysis contained 29 data points from the public archive and 6 data points from LaMontagne et al. (2022).

Metagenomic data processed previously (LaMontagne et al., 2022) with functions available in DADA2 v1.32.0 (Callahan et al., 2016), and custom scripts available here https://github.com/MgLaMontagne/dada2-optimize, was subsetted to the samples collected on August 25^th, August 30^th and September 8^th 2017. This sequence data, and the corresponding metadata (water quality and dPCR), was analyzed with scripts presented in Supplementary file S5. Analysis of variance (ANOVA) was conducted with base functions in R and further analyzed with functions available in agricolae v1.3-7 (de Mendiburu and Yaseen, 2020) and scripts also presented in Supplementary file S5. Metadata is available as Supplementary file S6. Multivariate analysis was conducted on vectors produced by non-metric multi-dimensional scaling (NMDS) with functions available in the R packages phyloseq v1.48.0 (McMurdie and Holmes, 2013) and vegan v2.6–6.1 (Oksanen et al., 2024) with scripts available in Supplementary file S7.

Total Bacteroidales copy numbers peaked in the immediate aftermath of Hurricane Harvey. Levels of this marker, across all sites, were highest on August 31 and September 8 (Figure 2A). At the first sampling after the rainfall ended, total Bacteroidales copies ranged between 5.7 × 10⁴–2.4 × 10⁵ copies 100 ml⁻¹ across the six sites. The following week (September 8), total Bacteroidales copy numbers ranged between 4.2 × 10³ and 1.5 × 10⁶ copies 100 ml⁻¹ across the six sites (Figure 2A). The average human Bacteriodales values ranged between 1.6 × 10³ and 4.0 × 10³ copies 100 ml⁻¹ and peaked on September 8 (Figure 2B). Following the initial peak, total Bacteroidales remained somewhat high following the hurricane with some additional elevated values observed in October. In contrast, human Bacteroidales values leveled off following the peak observed during the September 8 sampling. As previously reported in LaMontagne et al. (2022), E. coli numbers in these water samples peaked on August 31, averaging 1,087 MPN ml⁻¹ and then dropped to baseline within a week (Supplementary Figure S8).

Human source contributions, as assessed by metagenomic analysis, were greatest immediately after Hurricane Harvey (August 31) and then decreased precipitously and remained relatively low for all subsequent sampling dates (Figure 2C). To further confirm the presence of enteric-like bacteria, we extracted the sample-specific ASVs and determined taxonomy based on the Silva database (Quast et al., 2012; Yilmaz et al., 2014). Based on the bacterial abundance counts >5, bacteria similar to that found in the sewage source samples (Gammaproteobacteria and Bacteroidia) dominated all water samples collected immediately after flooding (August 31; LaMontagne et al., 2022).

Copy numbers of pBI143, intI1 and aint1 were significantly higher in water samples collected immediately after Hurricane Harvey than before or after the storm (Supplementary Table S3). The highest copy number for pBI143 (9.4 × 10⁴ copies/100 ml, Figure 3) was observed in samples collected immediately after Hurricane Harvey at station H (Figure 3). Lowest copies (6 copies/100 ml) for this plasmid were observed before Hurricane Harvey at station C, which is in Clear Lake (Figure 3). The highest copy number for intI1 (1.2 × 10⁵ copies/100 ml, Figure 3) was observed in samples collected immediately after Hurricane Harvey at station N. Lowest copies (7,008 copies/100 ml) for this ARG were observed before Hurricane Harvey at station K, which is also in Clear Lake (Figure 1). The highest copy number for aint1 (9,131/100 ml) was observed in samples collected immediately after Hurricane Harvey at station H (Supplementary Figure S1). Lowest copies (620 copies/100 ml) for this ARG were observed before Hurricane Harvey at station K. Copy numbers of sul1 did not differ significantly between sampling dates (Supplementary Figure S1).

Copy numbers of the plasmid pBI143 and ARGs (intI1, sul1 and aint1) all differed significantly (Supplementary Table S3) between stations. These differences between locations did not show a significant interaction with the date (pre/Hurricane Harvey/post, Supplementary Table S3). The highest copy number for pBI143, sul1 and aint1 were observed for samples collected at station H, where the geometric means across all dates for these markers were 29,775, 70,594 and 4,881 copies/100 ml, respectively. Copy numbers for pBI143 and sul1 at stations C and K (within Clear Lake) were at least an order of magnitude lower. The highest copy number for intI1 was observed for samples collected at station N, where the geometric mean across all dates at that station was 85,348/100 ml. These counts were significantly higher than the geometric mean for intI1 for samples collected at stations C and K.

NMDS analysis of metagenomic data showed a coherent cluster containing the six samples collected immediately after Hurricane Harvey (Figure 4). This cluster appeared driven by contamination with human waste, as assessed by copies of the plasmid pBI143 and the ARG intI1. That is these two parameters, which were measured by dPCR, corresponded to the structure of the microbial community, as assessed by sequencing of 16S amplicons. Salinity, temperature, total dissolved solids and oxygen levels also show significant relationships with bacterial community structure; however, dissolved inorganic nutrients and pH did not.