We selected sites around Clear Lake, an estuary between Houston and Galveston (Figure 1), based on the availability of long-term water quality data for these locations and ease of access. Water samples were collected on the afternoon of August 24th, 2017, 1 day before Hurricane Harvey landed in the Corpus Christi area. These samples, designated “pre” in this work, were collected at baseflow conditions (Supplementary Figure 1). A second set of samples was collected on August 30th. These samples, designated “HH” throughout, correspond to Hurricane Harvey samples and were collect hours after flow of a major tributary into Clear Lake peaked (Supplementary Figure 1). We added a second sampling site (N) when collecting the HH set to collect floodwater received by segment 1101C (Figure 1). Starting on September 8th we sampled six times to generate a “post” sample set; all post-HH samples were collected during baseflow conditions (Supplementary Figure 1). Samples collected in August and September were designated as summer season. Samples collected in October were designated as fall season. We also generated a mock sample by mixing raw sewage, collected as described previously (Amaral-Zettler et al., 2008), and surface water collected from station H (Figure 1) in March of 2018. The sewage and water were mixed at a ratio of 1 part sewage with 9 parts surface water.

At each sampling station, we collected surface water samples with a bucket lowered from an overpass or dock. Temperature and dissolved oxygen were measured in situ at 3–5 cm beneath the surface with a YSI model 55 dissolved oxygen (DO) probe (YSI Inc., Young Spring, OH). Water samples were split in the field for FIB (E. coli and Enterococci), metagenomic and nutrient analysis. For FIB analysis, unfiltered samples were transported on wet ice and stored at 4°C. Incubations for quantification of FIB were initiated within 24 h of sampling. FIB samples from pre-HH samples were discarded because we were locked out of our laboratory for several days and sampling holding times were exceeded.

Water samples collected before HH were stored on wet ice, returned to the laboratory and filtered to collect microbial samples for metagenomic analysis and to archive nutrient samples within 24 h of collection. Following HH, all samples were filtered in the field immediately upon collection. For metagenomic analysis, water samples were pulled through a Sterivex SVGPL10RC 0.2 μm cartridge (EMD Millipore, Billerica MA) until refusal (no flow at 15 psi), with a hand vacuum pump, as described previously (LaMontagne and Holden, 2003). The volume filtered, which ranged from 75 and 300 ml, was measured with a graduated cylinder. Sample filtrates and cartridges were transported on wet ice, temporarily stored at −20°C, and archived at −80°C. Two technical replicates were generated by collecting duplicate samples from station J on the eve of the storm and from the sewage-spiked samples described above.

Dissolved inorganic nutrient analysis for ammonium, orthophosphate, and nitrate were done by colorimetric analysis in microplates, as described previously (Ringuet et al., 2011). E. coli and Enterococci were enumerated using Colilert and Enterolert in the Quanti-Tray/2000 format following manufacturer recommendations (IDEXX, Westbrook, ME). Microbial community DNA for metagenomic analysis was recovered from the Sterivex cartridges as described previously (Amaral-Zettler et al., 2008) and assessed for molecular weight by agarose gel electrophoresis. These crude extracts were subsequently purified by passage through a OneStep PCR Inhibitor Removal column (Zymo, D6030) and purity was assessed by UV-spectra. Metagenomic analysis followed protocols outlined in the Earth Microbiome Project (Caporaso et al., 2011). Briefly, the V4 region of the 16S rRNA gene was amplified to generate an amplicon library. This library was multiplexed using Illumina designed indices, pooled with equal amounts, and sequenced on an Illumina MiSeq instrument as described by Caporaso et al. (2012).

Water quality data were analyzed and figures were generated with custom scripts presented in Supplementary Files 2, 3. These scripts included functions from packages from ggplot2 (Wickham, 2009). This data set included water quality data collected as described above and public data previously collected by the Texas Commission for Environmental Quality (TCEQ) and cooperating organizations. Data from the TCEQ archive was limited to samples collected between January 1st, 2011 and May 6th, 2021. For this time period, the mean for each of these five segments was calculated by grouping by segment, month and year. This data range was then merged with water quality data generated from samples collected in this study to generate Supplementary File 4.

MiSeq data were processed to determine alpha diversity of the microbial community using functions from DADA2 v. 1.20.0 (Callahan et al., 2016), with custom scripts presented in Supplementary File 5. Briefly, reads were filtered, trimmed, denoised, and merged to yield sequences from 251 to 253 nucleotides long. Chimeras were then removed with the function removeBimeraDenovo in DADA2 and putative non-chimeric sequences were assigned taxonomy and aligned with the functions IdTaxa and AlignSeqs in Decipher v 2.20.0 (Wright et al., 2012). Amplicon sequence variants (ASVs) and taxonomic identifications were merged to create a phylogseq-class object—available as Supplementary File 6—with functions in phyloseq v 1.36.0 (McMurdie and Holmes, 2013). ASVs with uncertain taxonomic identification at the phylum level were then removed before fitting the alignments into a phylogentic tree with functions in phangom v 2.7.1 (Schliep, 2010). Technical replicates (two samples collected at the same time and place but processed independently) were then merged and meta-data (volume filtered, environmental conditions, FIB counts, DIN, etc.) and reference sequences were combined to create a phylogseq-class object—available as Supplementary File 7—with functions in phyloseq and Biostrings (v 2.60.2 Pagès et al., 2021). Reference sequences were also exported in fasta format and compared to public sequences with the Seqmatch application (RDP Taxonomy 18) available through the Ribosomal Database Project (Cole et al., 2013). Default settings were used in Seqmatch. MiSeq data is available in the NCBI SRA under accession number/Bioproject ID: PRJNA795782.

Alpha diversity (richness and Shannon indices) was estimated with the plot_richness function after sewage-spiked samples were removed. Analysis of variance, calculated with a core function in R version 4.1.1, was used to test for significance of differences between samples collected before HH made land fall (sampled August 25th), vs. samples collected immediately after the storm (August 30th) and in September and October. Significance of differences was assessed with a Tukey test using the function HSD.test in the R package agricoloe v 1.3.5 (de Mendiburu, 2020). Coverage of the library of reads used for diversity analysis was visualized with the function rarecure in vegan v 2.5.7 (Oksanen et al., 2020).

The relationship between microbial community structure and nutrient levels was determined with correspondence analysis with custom scripts presented in R markdown in Supplementary File 8. This workflow started with a phyloseq object (S7). A prevalence threshold of 10% was set to remove rare taxa. Counts of the remaining 1,616 ASVs were transformed with Hellinger option prior to non-metric multidimensional scaling analysis (NMDS) with functions in vegan v 2.5.7 (Oksanen et al., 2020). NMDS was first run with all 44 samples, prior to removal of sewage-spiked samples. Goodness of fit of NMDS ordination was visualized with a Shepard plot generated prior to fitting meta-data to the ordination with functions in vegan. The resulting ordination plots were visualized with functions available in R package ggordiplots v 0.4.0 (Quensen, 2020).

Functional composition was predicted from the 1,616 ASVs used in NMDS analysis (above) with PICRUSt2 v 2.3.0-b (Douglas et al., 2020). To prepare the data, an ASV abundance table and fasta files were exported from a phyloseq object (S09) to a biom file (S10) and a sequence file (S11) using R package biomformat v 1.20.0 (McMurdie and Paulson, 2021). This pipeline, including bash scripts used in PICRUCSt2 analysis, are presented in Supplementary File 12.

Differential abundance of pathways (Supplementary Files 13, 14) predicted from ASVs and ASVs themselves was assessed using functions in R package ANCOM-BC v 1.2.2 (Lin and Peddada, 2020), following scripts presented in Supplementary File 15. Pathways that were differentially abundant were plotted with functions available in R package Heatplus v 3.0.0 (Ploner, 2021), following scripts presented in Supplementary File 16. Pathway functions and expected taxonomic range associated with them were identified with the web application MetaCyc v 25.5 (Caspi et al., 2013).

Hurricane Harvey lowered the salinity for Clear Lake. On the eve of HH (August 25th, 2017), surface salinity at stations C, K, and J, which correspond to water body segments 2425 and 2425B (Figure 1), ranged from 9 to 12 practical salinity units (PSU, Figure 2A). These pre-HH salinity levels are within the 95% confidence interval of the 10-year average for salinity for records for these two water body segments (2425 and 2425B, Supplementary Figure 1). Immediately after the storm, salinity dropped to < 1 PSU at all stations sampled herein and then gradually increased to pre-storm levels over the next 2 months (Figure 2A).

Oxygen levels in Clear Lake and tributaries to that system did not show a strong response to HH. Across all segments, oxygen levels averaged 6.8 mg/L before HH and 5.2 mg/L after. This temporal difference was not significant (p = 0.131); however, spatial differences between segments were significant (p = 0.0003). Oxygen levels averaged 6.9–7.1 mg/L for segments 2425 and 2425B, respectively, and less than 6 mg/L for stations 1101C, 1101, and 1113B (Figure 2B). Lowest oxygen concentrations were observed at station H in waterbody segment 1113B, where DIN concentrations are relatively high (see below).

Hurricane Harvey lowered the concentration of dissolved inorganic nutrients (DIN). On the eve of the storm, nitrate/nitrite (NO_x) ranged from 1 to 59 μM (Supplementary Figure 3). These pre-HH NO_x levels are in the range for records for the last 10 years for these segments (Supplementary Figure 3). Immediately after the storm, NO_x ranged from 1 to 9 μM and varied significantly between segments (p < 0.001) and type (pre-storm, HH, and post-storm). Highest levels of DIN were observed for samples collected from segment 1113B, which is approximately 50 m downstream of the outfall pipe of a wastewater treatment plant.

NO_x showed a non-conservative mixing relationship with salinity (Figure 3). That is the system is a sink NO_x. High concentrations (>10 μM) were associated with samples that showed salinities of 3 PSU or less. In contrast, samples with higher salinities (>3 PSU), typically showed NO_x levels of 3 μM or less, which suggests a freshwater source. Ammonium and phosphate levels showed a similar pattern with salinity as NO_x. High concentrations of ammonium (Supplementary Figure 4) and phosphate (Supplementary Figure 5) were associated with low salinities. DIN/P ratios were generally below 16 (Supplementary Figure 6). These ratios were on average highest (7–8) for segments 1101C and 1113B, respectively, and < 2 for segments 2425 and 2425B.

E. coli levels ranged from 488 to 1,733 MPN/100 ml for the six stations sampled on September 1st, 2017 (Figure 4), which was 72 h after HH pasted over the study area. The geometric mean (GM) for this set of samples was 1,018 MPN/100 ml. These values exceeded the statistical threshold value (STV) for single samples and GM recommended for water contact by the EPA (USEPA, 2012) and by the TCEQ for these particular water body segments (TCEQ, 2013). After 1 week, the E. coli levels had decreased to < 100 MPN/100 ml and remained relatively low until the end of October, when levels spiked again.

Enterococci levels ranged from 63 to 3,050 MPN/100 ml for samples collected in the fall of 2017; however, because of logistical issues, Enterococci levels were not measured until September 18th. For this period (post-HH), Enterococci levels did not differ between segments sampled (Supplementary Figure 7), and 22 of 24 samples exceeded the STV for single samples recommended for recreational water contact by the EPA (USEPA, 2012); GM (495 MPN/100 ml) exceeded, by an order of magnitude, the GM recommended for recreational water contact by the EPA (USEPA, 2012). FIB counts of samples taken in the post-HH period, were significantly higher (p < 0.001) than counts for the same segments collected over the last decade, where the GM of Enterococci was 74 MPN/100 ml.

Alpha diversity of the bacterial and archaeal community did not differ significantly between samples collected before and after HH (Figure 5). Average Shannon diversity indices ranged from 4.95 to 4.89 for samples after the event and averaged 4.44 for samples collected immediately before the storm. Diversity was relatively lower for samples collected before HH at stations 1101 and 1113B but only one sample was collected at that time point (Supplementary Figure 8). Average richness ranged from 761 to 633 for samples collected after the event and 488 for samples collected before (Supplementary Figure 9). Rarefaction analysis suggested the sequence library appeared to have the depth to describe alpha diversity (Supplementary Figure 10). After quality control, which included removing sequences that did not classify at the phylum level, average depth of the library was 112,178 reads. In other words, all 44 samples reached an asymptote.

Beta diversity of the bacterial and archaeal community structure, as assessed by NMDS, differed significantly between samples collected before and immediately after HH (Figure 6). The good fit (r² = 0.999) of a stressplot (Supplementary Figure 11), and low stress value (0.036), suggest this model is an excellent fit (Dexter et al., 2018). NMDS showed two clear clusters. Three samples collected before HH in segments (2425 and 2425B), that showed brackish salinities (8–12 PSU), formed a coherent cluster, with similarity to samples collected in the same segments in fall. Samples collected immediately after HH also formed a coherent cluster, with similarity to a sample spiked with sewage. One sample collected in March 2018 in segment 1113B clustered with the HH samples.

Recovery of Clear Lake progressed from the summer through the fall in segments (2425 and 2425B). These segments typically have brackish conditions. In the summer following HH the microbial community within these stations within Clear Lake looked similar, in terms of NMDS, to communities sampled from freshwater tributaries to the estuary (Figure 6). This recovery of the estuary’s microbiome appeared driven by salinity (Supplementary Figure 12). Salinity appeared strongly (p = 0.001) associated with pre-HH samples. NO_x and oxygen appeared strongly (p = 0.016 and 0.020, respectively) associated with post-HH samples. Phosphate also appeared associated with post-HH samples but the significance was weak (P = 0.086).

Bacterial community structure of the Clear Lake system shifted from a system dominated by Cyanobacteria before HH to a system dominated by Proteobacteria and Bacteroidetes immediately after (Figure 7). SAR324 clade (Marine group) was relatively abundant before HH and in the fall in segments we defined as brackish (2425 and 2425B). A total of 59 phyla were detected in 7,491 ASVs generated from 44 samples. Almost all of these ASVs (7,410/7,491) classified as bacteria. After removing taxa with relatively low (<10%) prevalence, almost all of the ASVs (1,534/1,617) were found to be differentially abundant, at a significance threshold of P < 0.05, in a model that tested the factors: salinity, NOx, and sample type (pre, HH, and post). In other words, these factors predicted the abundance of 95% ASVs.

Sample type predicted the majority of abundances. For example, the abundance of 819 (51%) ASVs differed between pre-HH and HH samples and the abundance of 1,007 (62%) ASVs differed between pre-HH and post-HH samples. Of the 10 most abundant ASVs observed in samples collected immediately after HH, nine classified as γ-Proteobacteria. Most of these (7/9) classified within the family Comamonadaceae and showed similarity to bacteria typically observed in freshwater systems. For example, ASV18 showed similarity to Limnohabitans curvus MWH-C1a, which was isolated from a lake (Hahn et al., 2010). The other two highly abundant ASVs (ASV64 and ASV23) showed similarity to γ-Proteobacteria isolated from rhizosphere soil (Jung et al., 2007) and freshwater systems (Hahn, 2003), respectively. ASV6, the most abundant ASV in libraries generated from floodwater samples, accounted for 5–17% of the reads generated in those six libraries. This ASV showed similarity to Aquirufa strains isolated from lakes (Hahn, 2006; Lee et al., 2018). The ASV showed the greatest differential abundance between pre-HH and HH samples (ASV103) showed high similarity to two uncultured bacteria (KP686762 and KP686755) generated from floodwater collected in North Carolina immediately after Hurricane Irene (Balmonte et al., 2016).

PICRUSt2 analysis predicted the abundance of 418 metabolic pathways from 1,616 ASVs generated from 44 samples. Differential abundance analysis suggested that 76 of these pathways were significantly different between samples. Cluster analysis, based on the relative abundance of these 76 pathways, suggested that floodwater samples formed a coherent group (Figure 8). That is, with one exception (sample eH), floodwater samples were relatively similar to each other in terms of predicted pathways. The outlying sample was also similar to samples collected immediately after HH in terms of numerically abundant phyla. In particular, eH and floodwater samples showed relatively high proportions of Proteobacteria and Bacteroidetes.

Comparison of pathways predicted from samples collected immediately before and after HH, identified 29 differentially abundant pathways (Supplementary Figure 12); 14 of these were significantly higher in samples collected before HH and 15 were significantly higher after the storm. Half (7/14) of the pathways associated with samples collected before HH were biosynthesis pathways. These include PWY-5347, which produces methionine and PWY-5840, which produces menaquinol-7. In contrast, only 3 of 14 pathways that were differentially abundant in samples collected immediately after HH were biosynthesis pathways, and two of these biosynthetic pathways are associated with virulence. PWY0-1338 confers resistance to the antibiotic polymyxin and PWY-6143 produces pseudaminic acid, which is associated with pathogenic Gram negative bacteria (Schirm et al., 2003). The vast majority (11/15) of pathways that were more abundant in samples collected immediately after HH were degradation pathways. These include pathways ORNDEG-PWY, ARGDEG-PWY, and ORNARGDEG-PWY, which are associated with degradation of L-arginine, putrescine, 4-aminobutanoate, and L-ornithine (Caspi et al., 2013).

Salinity and NO_x levels appeared associated with the abundance of 30 pathways. Of these 12 of were associated positively with salinity and 13 were associated negatively (Supplementary Figure 13). All but one of the pathways positively associated with salinity were biosynthesis pathways. These included five pathways (PWY-6165, 6349, −6350, −6654, −6167) associated with archaea and PWY-622, which is associated with starch biosynthesis by photoautotrophs (Caspi et al., 2013). In contrast, 6 of 13 pathways negatively associated with salinity were degradation pathways. These included two pathways (PWY-5427, −6956) associated with naphthalene degradation by bacteria and PWY-5088, which is associated with glutamate degradation by members of the Firmicutes phylum (Caspi et al., 2013). NO_x concentrations appeared associated with the abundance of five pathways (Supplementary Figure 14). The two positively associated pathways were degradation pathways; both are associated with mandelate degradation by Proteobacteria. Pathways negatively associated with NO_x concentrations include PWY-6174, which is associated with the mevalonate pathway in archaea, and PWY-5183, which is associated with toluene degradation by Proteobacteria (Caspi et al., 2013).