The data reported here were generated as part of a three-phase collection. During the first phase of sampling (No of Samples = 189; June 11 to October 8, 2014), sampling occurred irrespective of ambient or antecedent weather conditions; phases two and three were conducted to reflect dry weather (i.e., sampling was only conducted if <0.20 cm of rain fall during the preceding 48 h). Three waterways were the main focus of sampling (Figure 1): Ley Creek (No. of samples tested for ≥1 microbial targets (N) = 291; Number of sites where these samples were collected (S) = 58], Onondaga Creek (N = 1,700; S = 109), and Harbor Brook (N = 612; S = 38; Figure 2). Additionally, 1 tributary of Ley Creek [Sanders Creek (N = 1; S = 1)], and 4 tributaries of Onondaga Creek [Cold Brook (N = 34; S = 6), Hopper Brook (N = 48; S = 8), Kimber Brook (N = 4; S = 1), and West Branch (N = 126; S = 3); Figure 1] were sampled.

All sampled waterways follow a gradient from rural areas dominated by agriculture and forested to outlets on Onondaga Lake, which are located in the urban core of Syracuse, NY. The watersheds for the three waterways are 35 km² (Harbor Brook), 76 km² (Ley Creek), and 285 km² (Onondaga Creek), representing 5, 10, and 49% of the Onondaga Lake drainage area, respectively (Figure 1). Sample site locations along each waterway were selected to ensure accessibility for samplers, and to ensure both rural and urban sites were represented. A subset of samples (N = 31) were collected for microbial source-tracking marker (MST) analysis between August 2015 and August 2017 from Cold Brook (N = 8), Harbor Brook (N = 10), Hopper Brook (N = 10), and Onondaga Creek (N = 3). The sites used to collect MST data were selected to reflect sites where (i) more than one fecal source appeared to be present, including human and non-human (e.g., agricultural inputs, wetlands or ponds with large waterfowl populations), (ii) where past source-tracking efforts have been unable to identify sources of fecal contamination, and (iii) where fecal coliforms were consistently high during routine monitoring efforts. Samples from these sites were only tested for MST presence if fecal coliforms in a concomitantly collected sample were ≥ 200 cfu/100 mL.

To characterize the impact of land use in the area immediately around each sample collection site the proportion of land within 122, 366, and 1,098 m of each site under pasture-hay, cropland, forest-wetland, and developed (>20% impervious) cover was determined [for an example of the code used see: https://github.com/wellerd2/Calculating-land-use-land-cover-and-landscape-structure-parameters] (Rehmann and Soupir, 2009). Briefly, land use data was extracted from the National Land Cover database (https://www.mrlc.gov/national-land-cover-database-nlcd-2016) for the relevant year, and buffer distance. Buffer distances were selected based on Leafy Green Marketing Agreement (https://lgmatech.com/wp-content/uploads/2020/08/CA-LGMA-Metrics-August-2020_Final_Clean_9-18-20.pdf) recommended buffers from the location of any adjacent land uses that are likely to present a food safety risk (Table 7 from Crop Land and Water Sources Adjacent Land Uses). For example, the recommended distance from cropland to a concentrated animal feeding operation (CAFO) with >1,000 head is ~1,200 feet (366 m). We selected this buffer, and buffers 1/3 smaller (122 m) and larger (1,098 m) to obtain land use data. Additionally, to compare water quality between rural vs. urban areas, sites were characterized as urban if they were within the boundaries of the city of Syracuse (Figure 2). It is also important to note that four of the five sites to the north of the city limits, while in a predominantly developed landscape are in patches of non-developed land cover. However, recognizing that dichotomizing into rural vs. urban oversimplifies that complex landscape around Syracuse, separate analyses using the amount of developed cover around the site were also implemented.

Each site was sampled up to twice a month between May and November, with samples being collected annually from 2008 to 2017 (except 2010 or 2011). Sampling within each waterway was performed in a downstream to upstream manner to avoid sampling the same slug of water. All samples were collected at the same time, varying volumes outlined were collected in different bottles for the various assays including microbiological, physiochemical, and nutrient/sediment analyses.

For microbiological analysis, direct grab samples were collected from the centerline of the waterway by submerging a sterilized 150-mL plastic bottle just below the water surface. Approximately 125 mL of each water sample was transferred to a sterilized and pre-preserved (with sodium thiosulfate) 150-mL coliform bottle. Approximately 10 mL of each water sample was transferred to a 25-mL glass vial for chloride analysis (see below). When samples were collected for MST analysis, a separate 1-L grab sample was also obtained using the direct grab method. Total coliform, fecal coliform, E. coli and/or Enterococcus levels were enumerated using the 125-mL sample <8 h after collection. FIB and total coliform levels were enumerated by the Onondaga County Health Department (a certified laboratory) using NY State Department of Health standard methods (Supplementary Table 1). Separately, samples collected for MST detection were processed <24 h after collection by the Wadsworth Laboratory (NY State Department of Health, Albany, NY; New York State Department of Environmental Conservation, 2017). Briefly, after undergoing anerobic enrichment to select for members of Bacteroidales (Green et al., 2012, 2014, 2019), separate PCR-screens were performed, for host-specific avian (Lu et al., 2008; Green et al., 2012), canid (New York State Department of Environmental Conservation, 2017), human (Bernhard and Field, 2000a; Shanks et al., 2007, 2009; Green et al., 2014), and ruminant fecal markers (Bernhard and Field, 2000b; Supplementary Table 1).

At each sampling event, data were also collected on physicochemical water quality parameters and weather factors. Data on water temperature, pH, dissolved oxygen, specific conductance, and turbidity were collected in situ using a YSI 650 MDS handheld device (YSI Inc., Yellow Springs, OH, USA) equipped with a 6600 or 6820-V2 multi-parameter water quality probe. Chloride was measured in situ using a Hach™ pocket colorimeter. To enumerate nutrient and sediment levels, a separate 7.6 L sample was collected from the same location, and at the same time as the sample used to enumerate FIB levels. The 7.6 L sample was divided into separate 1 and 0.5 L aliquots for nutrient and sediments analysis, respectively (Supplementary Table 1). It is important to note that not all parameters were measured in all samples (e.g., due to equipment malfunctions, availability of staff, time conflicts, feasibility of obtaining samples; S2). Weather data (air temperature, rainfall, relative humidity, and wind speed) for each sampling event were obtained from the Network for Environmental and Weather Applications (NEWA). Data was obtained using the weather station closest to each sampling site location (http://newa.cornell.edu/).

All analyses were performed using R version 3.6.2 (R Foundation for Statistical Computing, Vienna, Austria). For all analyses, all data were log₁₀-transformed. While 2,816 samples were collected, not all analytes were enumerated in all samples; as a result, 5-point summaries and proportion of missingness were calculated for each variable (Supplementary Tables 2, 5, 6).

Bayesian mixed models were developed to (i) characterize spatial and temporal patterns in log₁₀ FIB levels, and human and ruminant MST markers presence, and (ii) identify environmental factors associated with log10 FIB levels or probability of MST detection. Log-linear Bayesian mixed models were developed when the outcome was log10 E. coli, Enterococcus, fecal coliform, or total coliform levels, while Bayesian linear mixed models were developed when the outcome was probability of human or ruminant MST detection. To characterize spatial and temporal patterns, models were developed with fixed effects of elevation, latitude, longitude, rurality (i.e., if the sample was collected inside or outside Syracuse, NY city-limits), and waterway, and of year and month, respectively. For models with a spatial fixed effect, month was included as a random effect, while waterway was included as a random effect for models with temporal fixed effects.

To identify relationships between fecal indicators and environmental factors, Bayesian mixed models were developed with random effects of month and waterway. These models included a single fixed effect for land use (proportion of land within 122, 366, or 1,098 m of the sampling site under crop, pasture-hay, forest-wetland or developed cover), weather (air temperature at sampling and 0–3 d before sampling, rainfall 0–1, 1–2, and 2–3 d before sampling, and wind speed 0–1 d before sampling), or water quality parameters (log10 chloride, conductivity, total dissolved solids, nitrate, total organic carbon, total phosphorous, salinity, total suspended solids, and turbidity levels as well as pH). Since not all analytes were measured in all samples, models were implemented only if there were at least 20 pairwise observations for both the microbial target and environmental factor of interest (Supplementary Tables 5, 6). Due to the large number of samples tested for fecal coliforms (N = 2,658) compared to the other FIBs, fecal coliform models could be implemented for all covariates considered here. Conversely, Enterococcus models could not be implemented with chloride, total dissolved solids, total organic carbon, or total phosphorous as the explanatory variables.

Models were fit using the brms package, uninformative priors, 3 chains, and thinning set to 10 (Bürkner, 2017, 2018). While the number of iterations per chain was set to 5,000 (burn-in of 2,500) for most models, there were convergence issues for a subset of models (Bürkner, 2017, 2018). For those models with convergence issues, the number of iterations per chain and burn-in iterations were increased per the package author's recommendation (Bürkner, 2017, 2018). The median, mean, maximum a posteriori (MAP), and 89% credibility interval for the effect estimates were calculated using the bayestestR package (Makowski et al., 2019a). Bayesian credibility intervals are not interpreted in the same way as frequentist confidence intervals. Instead the 89% credibility interval indicates the highest probability portion of the posterior distribution of possible effect estimates and should be interpreted as: “Given the observed data, the effect estimate has a 89% probability of falling between x and y” (Makowski et al., 2019a).

Unlike frequentist regression models, Bayesian models are not based on assessing statistical significance; instead, Bayesian models offer a probabilistic (more flexible) view of the parameters and their corresponding uncertainty. As such, Bayesian analyses are more robust to small sample size, do not rely on a researcher-specified p-value threshold, and are less prone to Type I errors compared to frequentist approaches. Therefore, rather than concluding that an association is present when the P-value is below a certain threshold, conclusions can be drawn by determining if the parameter is outside of a range of practically negligible effect. This measure, called the ROPE percentage, indicates the magnitude of effect, and is calculated by determining the percent overlap between the 89% credibility interval and the range of practically no effect. Thus, the closer the ROPE percentage is to 0, the more confident we can be that the given factor has a substantial effect on FIB levels or probability of MST detection. Specifically, we use the following cutoffs for ROPE interpretation: >99% ~ negligible effect, >97.5% ~ probably negligible effect, between 2.5 and 97.5% ~ uncertain effect, <2.5% ~ non-negligible effect, <1% ~ significant effect (Makowski et al., 2019a). Associations in Bayesian regression can also be assessed using the probability of direction (PD), which is an index of if a positive or negative effect exists regardless of if that effect is negligible or non-negligible. The PD correlates strongly with frequentist P-values with PD values near 1.0 indicating greater certainty that the effect of the factor is truly positive or negative (i.e., indicates confidence in the direction of the association (Makowski et al., 2019a,b). Specifically, PD values of 0.95, 0.975, 0.995, and 0.9995 correspond to two-sided frequentist P-values of 0.10, 0.05, 0.01, and 0.001, respectively (Makowski et al., 2019a,b). The ROPE and PD measures are independent, and a factor can have a high PD, but low ROPE value, indicating high certainty in the direction of a negligible effect. The inverse would indicate that the parameter has a non-negligible (or substantial) effect, but the direction of that effect is unclear. Overall, PD and ROPE can be used in conjunction to make a statement of the form, “the effect of the given factor has a probability of PD of being negative/positive, and can be considered as significant/nonsignificant (Makowski et al., 2019a).” Measures of (PS) were also calculated, and indicate the probability that the parameter's effect is above a given threshold, representing a negligible effect in the median's direction; this is a unidirectional equivalence test that indicates if the effect is both non-negligible, and in a given direction (Makowski et al., 2019a,b). Values should be larger than 0.5 to indicate practical significance; a cut-off of 0.75 was used here to be conservative. Based on the recommended cut-offs and interpretations for PD, PS, and ROPE, we determined that with PS >0.50, PD >0.75, and ROPE <0.25 warranted reporting here.

Fecal coliform levels were positively associated with Enterococcus and total coliform levels (Table 1, Supplementary Table 3). While E. coli levels were associated with levels of other FIB, land use and spatial factors, such as rurality, that were positively associated with fecal coliform and Enterococcus levels tended to be negatively associated with E. coli levels, and vice versa (Tables 1–3). Of the 31 samples tested for host-specific microbial source tracking markers (MST), 22 (71%) and 11 (35%) samples were positive for human and ruminant markers. Detection of human markers were positively associated with detection of ruminant markers (PD = 0.91; PS = 0.80); however, the magnitude of this effect is uncertain (ROPE = 0.16; Supplementary Table 3). Additionally, neither marker was associated with levels of any of the other microbial targets considered here (PS <0.50, PD <0.95, and ROPE >0.025). This may indicate that (i) human and ruminant fecal contamination were coming from the same sources, (ii) that sources of FIB were different from sources of human and ruminant contamination, and/or (iii) that other host-specific sources of contamination not tested here were driving fecal contamination in the present study. For example, past studies have found that FIBs, such as E. coli can naturalize and survive in non-host environments (e.g., algal mats, water, soil; (Hendricks, 1967; Byappanahalli et al., 2003; Ksoll et al., 2007; Nautiyal et al., 2010; Goto and Yan, 2011; NandaKafle et al., 2018)). As a result, FIB levels may not reflect recent fecal contamination events, resulting in a lack of association between host-specific markers and FIB levels. It is also interesting to note that the findings reported here are inconsistent with a previous study conducted on Onondaga Creek (NY, US), which found a significant positive association between human marker detection and Enterococcus levels, and a significant inverse association between detection of human and ruminant markers (Green et al., 2019). This discrepancy may be due to the fact that this former study collected samples from Onondaga Creek in 2015, while the present study represents a larger number of waterways and a longer sampling period (2008–2017). Water quality is known to vary over space and time, and to vary by the spatial scale of analysis (e.g., Goyal et al., 1977; Pandey et al., 2012; Partyka et al., 2018a; Fluke et al., 2019; Weller and Jordan, 2020; Weller et al., 2020b). It may also indicate that processes driving fecal contamination differ for each of the watersheds sampled here, complicating detection of watershed-specific patterns of contamination (which were the focus of the previous study on Onondaga Creek, NY, US). However, since only three samples from Onondaga Creek and 31 samples total were tested for MST markers, examining watershed-specific patterns of association were outside the scope of the study reported here, and should be investigated in future studies on Onondaga Lake, and other watersheds.

The distribution of FIB was non-uniform with Enterococcus and fecal coliform levels being higher and E. coli levels being lower, on average, within the city of Syracuse, NY, compared to outside it (i.e., urban vs. rural samples; Figure 1; Table 2); neither the probability of detecting human or ruminant MST markers were associated with rurality (Supplementary Table 3). In fact, while E. coli levels were, on average, 0.35 log10 CFU/100-mL (89% Credibility Interval [CI] = −0.61, 0.00; PF = 0.94; ROPE = 0.10) lower, and Enterococcus and fecal coliform levels were, on average, 0.32 log CFU/100-mL (CI = 0.13, 0.52; PD = 0.99; ROPE < 0.01) and 0.51 log CFU/100-mL (CI = 0.44, 0.59; PD = 1.00; ROPE < 0.01) higher in samples collected from urban, compared to rural sites. Since latitude increased and elevation decreased along the rural-urban gradient, it is not surprising that fecal coliforms and Enterococcus were positively associated with latitude and negatively associated with elevation (Syracuse, NY is farther north and at a lower elevation, than surrounding rural areas), while the reverse patterns were observed for E. coli, and human and ruminant MST markers (Table 2). Overall, these patterns are also consistent with the land use associations observed here. Specifically, the only association between land use and E. coli levels was a positive association with pasture cover (Table 3), which mirrors the positive association between probability of detecting human and ruminant MST markers and pasture cover (Supplementary Table 4). Conversely, Enterococcus and fecal coliform levels were negatively associated with agricultural land uses and forest-wetland cover, but positively associated with developed cover (Table 3). Interestingly, while the presence and magnitude of the association between E. coli levels and pasture cover appeared robust to buffer size, the associations between Enterococcus and fecal coliform levels and pasture cover, and between ruminant detection and pasture cover were only evident for 1,098 and 122 m buffers, respectively (Table 3, Supplementary Table 4). The relationship between total coliform levels and cropland, human marker detection, and cropland, and ruminant marker detection and developed and forest-wetland cover also seemed dependent on buffer size. Thus, while there were associations between specific land uses and each fecal indicator, these associations appeared dependent on the scale of analysis. This is not unexpected, as a previous study that examined nutrient pollution in Maryland waterways also found that spatial scale of analysis affected model results (Weller and Jordan, 2020). Similarly, a study in Arkansas, US that characterized associations between E. coli levels and riparian land use found that the size and area of the riparian buffer affected the change-point (i.e., the point within the land use parameter that demarcates a significant change in E.coli levels; (Scott et al., 2017)). Since the buffers used here were selected to replicate the recommended distance between large livestock operations (>1,000 head) and agricultural water sources in the Leafy Green Marketing Agreement (i.e., 366 m), this finding has direct implications for produce safety, and highlights the importance of considering spatial scale (e.g., using 122, 366, or 1,098 m buffers), when developing strategies and/or guidance for identifying potential fecal contamination risks.

Regardless of the scale for analysis, the study reported here identified strong associations between adjacent land use and fecal contamination, which is consistent with past studies conducted in the Northeastern US, and other areas (Shiels and Guebert, 2010; Wilkes et al., 2011; Duris et al., 2013; Scott et al., 2017; Wu, 2019; Weller et al., 2020b). For instance, we found that E. coli levels increased by 0.30 (89% CI = 0.17, 0.39; PD = 1.00; ROPE <0.01) log10 CFU/100-mL for each percent increase in pasture cover within 1,098 m of the sampling site (Table 1); this positive linear association is consistent with the findings of the aforementioned Arkansas (US) study (Scott et al., 2017). Conversely, Enterococcus and fecal coliform levels both decreased as pasture cover increased, but increased as developed cover increased (Table 1). Given the pattern of land use associations observed here, our findings suggest that the sources of Enterococcus and fecal coliform contamination differ from sources of E. coli or MST markers, with the former coming from developed areas, and the latter from agricultural areas. While we would expect human markers to be associated with urban areas and developed land uses, the association between human markers and agricultural land uses and rural areas may indicate the presence of human fecal contamination sources in rural hinterlands. For example, areas within Syracuse are on municipal water and sewer, while the rural areas to the south largely utilize septic; as such one source of fecal contamination in rural areas may be failing septic. While such a hypothesis is supported by past studies conducted in Queensland, Australia (Ahmed et al., 2005), Georgia, US (Sowah et al., 2014, 2017), Michigan, US (Verhougstraete et al., 2015), and New York, US (Rao et al., 2015), data on the age, location, and density of rural fecal contamination sources, such as septic systems, were not available and this could not be investigated in the present study. This highlights a specific research gap, the impact and source of rural fecal contamination, that future studies in the Onondaga Creek Watershed should address.

Despite differences in land use associations for the various fecal indicators, levels of all indicators were higher in Onondaga Creek, compared to the other waterways sampled (Table 2, Supplementary Table 3). For instance, E. coli levels in Kimber Brook were almost 2 log10 (89% CI = −3.29, −0.31; PD = 0.97; ROPE<0.01) lower, on average, than levels in Onondaga Creek. These findings support previous studies' conclusions that Onondaga Creek was responsible for most of the bacterial and sediment loads entering Onondaga Lake (Steven, 1996; Prestigiacomo et al., 2007). Given the consistency across fecal indicators, including the host-specific markers, this suggests that, independent of land use near each sampling site, fecal contamination patterns varied between the waterways sampled here with Onondaga Creek being more impaired than the other waterways. Thus, while this finding is not surprising, it is important for guiding efforts to mitigate fecal contamination of surface waters that flow into Onondaga Lake. Specifically, since waterways sampled here drain into Onondaga Lake, this finding suggests that mitigation efforts should target Onondaga Creek, as opposed to Ley Creek, Harbor Brook, or Hopper Brook.

All FIB tested for in the present study were positively associated with nutrient (e.g., nitrate, phosphorous), chemical (e.g., salinity, conductivity), or sediment (e.g., turbidity, total suspended solid) contamination. For example, the concentration of E. coli, Enterococcus, and fecal coliforms increased by 0.16 (CI = 0.06, 0.26; PD = 0.99; ROPE = 0.11), 0.68 (CI = 0.08, 1.24; PD = 0.96; ROPE = 0.01), and 0.33 (CI = 0.26, 0.41; PD = 1.00; ROPE < 0.01) log10 CFU/100-mLs, respectively, for each log10 increase in nitrate levels (Table 1, Supplementary Table 3). Similar patterns were observed for conductivity, salinity, total organic carbon, total phosphorous, total dissolved solids, total suspended solids, and turbidity (Table 1). In fact, Enterococcus and fecal coliform levels increased by ~1 log for each log10 increase in total suspended solid levels (89% CI = −0.48, 2.39; PD = 0.89; ROPE = 0.03) and for each log10 increase in total organic carbon levels (89% CI = 0.71, 1.48; PD = 1.00; ROPE < 0.01), respectively (Table 1). The discovery of positive associations between physicochemical and microbial water quality parameters is consistent with the previous studies (Soupir et al., 2010; Viau et al., 2011; Wilkes et al., 2011; Rao et al., 2015; Verhougstraete et al., 2015). For example, a survey of 64 Michigan, US rivers found that nutrient levels were able to account for 48% of variance in E. coli levels (Verhougstraete et al., 2015). Additionally, a study that identified factors associated with detection of parasites of fecal origin in an Ontario, Canada watershed, found that Giardia detection was associated with both conductivity and total phosphorous levels (Wilkes et al., 2011). Viau et al. surveyed the microbial and physicochemical quality of Hawaiian streams, and found evidence of positive associations between multiple microbial and physicochemical parameters, including E. coli and total phosphorous levels, and between enterococci levels and turbidity (Viau et al., 2011). In fact, a positive association has been repeatedly identified between microbial water quality and turbidity levels in studies conducted in different water types (e.g., reservoirs, canals, ponds), regions (e.g., the Northeastern and Southeastern US, Ecuador), and where different study designs were used (e.g., microbial targets, sampling strategies; (Viau et al., 2011; Francy et al., 2013; Rao et al., 2015; Topalcengiz et al., 2017; Partyka et al., 2018b; Weller et al., 2020a)). Given the reproducibility of the turbidity-microbial water quality relationship, multiple studies have suggested using turbidity as a supplemental indicator for identifying when a waterway may be contaminated by FIB. Overall, the positive association between microbial water quality (e.g., FIB levels), and nutrient (e.g., nitrate) and sediment (e.g., turbidity) levels in this and other studies, may indicate that similar processes are driving microbial, sediment and nutrient contamination. Such a conclusion is supported by literature on sources, fate, and transport of microbial, nutrient, and sediment contaminants. For instance, it is well-established in the scientific literature that rain events facilitate the release and transport of microbes, nutrients, and sediment from sources to surface waterways (Nagels et al., 2002; Muirhead et al., 2004; Zhou et al., 2017; Fluke et al., 2019). Indeed, rainfall antecedent of sampling events, and wind speed (a proxy for storm events) were associated with elevated FIB levels in this and previous studies (Francy et al., 2013). Additionally, the reported study also highlights that even though sampling in phases two and three of the present study only occurred during dry weather, bacterial pollution was still persistent in sampled waterways. Moreover, bacteria are often bound to particles, with one study reporting that 49% of all Enterococci detected were attached to soil or manure particles (Soupir et al., 2010), and streambed sediments have been repeatedly identified as an in-channel store for FIBs, whose resuspension can result in elevated FIB levels in the water column (Goyal et al., 1977; Muirhead et al., 2004; Rehmann and Soupir, 2009; Kim et al., 2010; Zhou et al., 2017). Thus, differential persistence of FIBs in streambed sediments, may account for the lack of correlation between FIB groups in the present study. From an applied perspective, these findings suggest that strategies to mitigate fecal contamination may concurrently reduce nutrient and sediment contamination; this is of specific interest for land managers in the Syracuse (NY, US) area, given the aforementioned interest in improving the quality of water entering Onondaga Lake; as well as, in regions with similarly impaired waterways.