The data used for this study were originally collected and processed by Marquette University researchers (Beattie et al., 2018) based on (1) the number and spatial distribution of samples and (2) the high detection rate of ARGs in riverbed sediment and surface water. The samples were obtained from Kewaunee, Ahnapee, and East Twin Rivers in Kewaunee County, Wisconsin during five sampling events representing four seasons and 20 sampling sites (sampling sites are depicted as circles in Figure 1). Sites were selected based on the impacts of expected variability from the dense livestock agriculture in Kewaunee County and public access to the sites (Beattie et al., 2018). DNA was extracted from samples, and quantitative PCR was used to quantify AMR-associated genes [i.e., erm(B), qnrA, tet(W), sul1, and intI1] and the 16S rRNA gene. Antibiotic resistance gene selection is detailed in a previously published study (Beattie et al., 2018); briefly, genes resistant to antibiotics commonly used in agriculture (tetracycline, macrolides, and sulfonamides; tet(W), erm(B), and sul1), genes found on mobile genetic elements (intI1 and sul1), and genes conferring resistance to clinically important antibiotics (fluoroquinolones; qnrA) were chosen to explore the diversity of environmental resistance in the study area. Values measuring below the detection limit (< 8% of the data points; see Supplementary Table S1) were set to half the detection limit. Detailed sampling methods, DNA extraction, qPCR protocols, and the full ARG dataset (relative and absolute abundances) can be found in a previously published study (Beattie et al., 2018). River network, precipitation, temperature, and source location data and processing can be found in previously published studies (Wiesner-Friedman et al., 2021a,b).

This study uses a physically meaningful LUR model (i.e., source terms are not allowed to have negative regression coefficients) implemented in FIT without modifiers (i.e., attenuators or amplifiers) to focus on characterizing source contributions to each ARG response, y_i (Wiesner-Friedman et al., 2021b):

The ARG response, y_i, represents log10 of the relative abundance z_i, where z_i is the ratio of ARG copies to 16S rRNA gene copies. The first three variables, namely, P1_i, P2_i, and Freezing_i, represent recent precipitation, antecedent precipitation, and a seasonal, Bernoulli-distributed variable representing whether the monthly average temperature was freezing for each i^thsample. The first term, β₀, and the last term, ε_i, represent the regression intercept and random error, respectively. β₁, β₂, and β₃ are associated regression coefficients, and ∑u=1Uβusi(u)(α(u)) is a sum of contributions from spatially distributed sources. Each source term, si(u)(α(u)), is defined as a function of ground hauling, overland, and downstream decay hyperparameters; α^(u), the flow-connected distances (i.e., overland and downstream distances) from the spatial locations of sources of type u (e.g., WWTPs) to the sampling locations of each i^thsample; and sample site flow (proxied by Strahler stream order) (Wiesner-Friedman et al., 2021b). The sources can also be weighted by information associated with their scale (e.g., size of the land cover area, gallons of manure, or equally weighted) (Wiesner-Friedman et al., 2021a). Because hyperparameters describe the distance decay of dimensionless quantities that denote the scale of sources, source terms describe dimensionless source contributions to sampling sites. All source terms si(u) are SPMs equal to the z-scored Sum of Exponentially Decaying Contributions (SEDC) so that a one standard deviation increase in the u^th SEDC represents a β_u increase in the response. An RAR expresses the ratio of relative abundances for a one standard deviation increase in source u. In other words, a one-unit increase in SPM si(u) results in a (100(RAR^(u)−1)) percentage increase in relative abundance z (Wiesner-Friedman et al., 2021b). See Supplementary material S5 for details.

Multiple spatial databases are available to represent potential sources in Kewaunee County. While many sources of elevated ARGs exist globally, we identified thirteen categories of potential sources that may be important to Kewaunee County (u = 1,2,…,13).

The first eight potential sources of the u^th categories were spatially related to sampling site locations with SPMs, which are the z-scored overland and river-distance with flow (ORF) SPM: (1) AFOs, (2) manure application fields, (3) septic systems, (4) industrial land application sites, (5) municipal or septage (i.e., domestic-originating) land application sites, (6) WWTPs, (7) low-intensity developed land cover representing rural accumulations on imperviousness, and (8) high-intensity developed land cover representing more urban/residential accumulations on imperviousness. The next three potential sources were soil sources represented by the dominant soil type within a 1-km radius of the sampling location: (9) Type A represents soils with the highest infiltration rate when saturated, likely consisting of sand, sandy loam, loamy sand, or gravel soil types; (10) Type C represents soils with a low infiltration rate when saturated, likely consisting of clay loam, silty clay loam, sandy clay, or silty clay; and (11) Type D represents soils with a very low infiltration or high run-off potential, likely consisting of clay loam, silty clay loam, sandy clay, or silty clay. The last two source categories were those related to the sampling site locations using the ground hauling, overland, and river distance with flow (GORF) SPM. Here, the SPM leverages two source locations to capture how land application occurs on land disproportionately closer to the origin of the waste [e.g., manure is hauled from AFOs (u = 1) to manure application fields (u = 2), and AFOs will minimize hauling distances for cost purposes (Hadrich et al., 2010)]. The sources defined with this SPM are 12) AFOs via the ground transport or hauling of manure to fields and 13) domestic-originating land application sites (u = 5) with waste that is applied more greatly in proximity to denser residential areas represented by septic system locations (u = 3). See previous study for SPM equations and descriptions (Wiesner-Friedman et al., 2021b).

The microbial FIT framework is a three-step approach (Wiesner-Friedman et al., 2021b) using the physically meaningful LUR (Equation 1). The three stages of the FIT framework, Find, Inform, and Test, were applied independently to the 10 ARG responses (i.e., sediment and water measurements of five ARGs). This modeling was implemented in MATLAB 2020b (https://scicrunch.org/resolver/RRID:SCR_001622).

In Find reliable databases of spatially distributed sources, many candidate databases can be explored for their ability to reproduce source-term relationships for unseen data. Here, for each of the potential ARG source categories (u = 1,2,…,13), we explored the candidate options for representing the sources based on (a) available databases [e.g., the Wiscland-2 land cover database, the Wisconsin Pollution Discharge Elimination System (WPDES) database, and county databases available for manure storage and septic systems], (b) coding options, (c) different classes of data, and (d) weighting options (see Supplementary material S2 for details).

The find stage obtains a reliability score for each candidate database option. The reliability score is calculated by obtaining hyperparameters for the SPM for a training set of response data. Using training hyperparameters for a test set (see Supplementary material S4), we calculated a reliability score with test-set regression coefficients. The score rewards candidate databases with the most consistently positive test-set regression coefficients [i.e., sign stability score (SSS)], the largest sum of test-set coefficients (M), and the lowest variability of the test-set coefficients (Wiesner-Friedman et al., 2021b). See Supplementary material S2 for the reliability score equation.

In the inform stage, source terms are informed with transport-characterizing hyperparameters in ORF and GORF SPMs. In this stage, hyperparameter values, α^(u) (i.e., describing average ground hauling distance γ_G, decay overland α_O, and downstream α_R), were obtained by independently maximizing the RAR (10^β(u)) for every source type. This maximization was subject to a penalty on very low or high values of overland decay and ground hauling hyperparameters, α_O and γ_G, that yielded poor regression qualities (Wiesner-Friedman et al., 2021b). The spatial predictor models leveraging the hyperparameters, α^(u), are described in detail in the microbial FIT framework (Wiesner-Friedman et al., 2021b).

The goal of the last stage is to test the predictive ability of informed source terms to identify key sources of elevated ARGs in Kewaunee County Rivers and Streams. Before model selection, if the correlation of informed source terms was large (ρ≥0.7), those yielding the highest univariate R-squared were chosen over correlated options (Dormann et al., 2013). Then, seasonal, precipitation, and source terms were stepwise selected with AIC. The prediction of individual source impacts on unsampled locations was initially assessed through the find stage based on the robustness of the source term's consistent contributions and the ARG's log10 relative abundance across training and test sets. Ultimately, the prediction of total source impacts was assessed with the adjusted R² resulting from the test stage.

To understand which ARGs are likely to be shed by bovines based on antibiotic usage, we reviewed information from a thorough systematic review outlining and quantifying antibiotic usage on Wisconsin dairy farms (Pol and Ruegg, 2007). Since the publishing of that study, practices may have changed, and a previous ARG study found local knowledge to be beneficial (Rogers et al., 2018). We conducted a small survey with dairy cattle veterinarians and asked Wisconsin dairy cattle veterinarians about the changes they have observed in antibiotic use over their careers. Google Maps was used to search for veterinarians, veterinary services, and “large animal” services. In total, 12 veterinary offices were identified as serving the dairy cattle industry, representing the practices of over 30 veterinarians. Surveys were sent to all 12 offices, and four returned our surveys. See Supplementary material S6.1 for survey questions and veterinarian responses.

According to the four veterinarians we interviewed, ceftiofur, a beta-lactam (Dowling, 2004), was the most frequently prescribed antibiotic for disease treatment (Dowling, 2004). Other antibiotics that veterinarians prescribe are enrofloxacin, florfenicol, tulathromycin, and oxytetracycline. These interviews (see Supplementary material S4 for interview details and discussion) reflected antibiotic use for disease treatment that was consistent with the findings of the 2007 study on antibiotic usage at AFOs in Wisconsin (Pol and Ruegg, 2007). These antibiotics belong to the broader classes of fluoroquinolones, sulfonamides, macrolides, and tetracyclines, which correspond well with resistance encoded by or co-occurring with the ARGs used for this LUR modeling (Pal et al., 2015; Beattie et al., 2018). Our interviews did indicate that some antibiotics are still used preventatively (e.g., tetracycline flushes) because operations have difficulty monitoring large herds. Preventative antibiotic use is therefore a missing component in understanding the frequency and dose associated with these classes of antibiotics. However, our interviews and the 2007 study provide a biological link between the panel of ARGs from our study and dairy AFOs. Additionally, although intI1 has been identified as a marker for anthropogenic pollution more broadly, clinical class I integrons can carry ARG cassettes, conferring multidrug resistance with the ability to spread rapidly through horizontal gene transfer (Gillings et al., 2015).

We could not directly obtain information on clinical antibiotic usage in rural Wisconsin. However, one wastewater study (Karthikeyan and Meyer, 2006) indicates a total of six classes of antibiotic compounds found in the influent of municipal wastewater from across Wisconsin. Based on the frequency of the detection of different antibiotic classes, the study indicates that the ARGs for this study may also well represent the anthropogenic impacts of clinical antibiotic use.

After implementing the FIT framework across five ARGs in sediment and surface water, we found that all five of the ARG responses (see Table 2) are positively associated with bovine sources (i.e., AFOs, AFOs via the ground hauling of manure to application fields, or manure application field locations) in at least surface water or sediment. For ARG responses in sediment, FIT selected the GORF AFO source term, which represents the contributions from AFOs via the ground hauling of manure to application fields. A one standard deviation increase in GORF AFO contributions was associated with RARs of 1.58 and 2.01, representing 58% and 101% increases (p < 0.10) in the relative abundance of erm(B) and tet(W).

Across ARG responses in surface water, FIT selected the GORF AFO source term for erm(B). A one standard deviation increase in these AFO contributions was associated with a 45% (p < 0.05) increase in the relative abundance of erm(B) in surface water. For tet(W) and qnrA responses, FIT selected the ORF AFO source term, representing the contributions directly from AFO locations. A one standard deviation increase in AFO contributions was associated with 77% (p < 0.05) and 49% (p < 0.05) increases in the relative abundance of tet(W) and qnrA, respectively. For sul1 and intI1, FIT selected the ORF manure application fields representing contributions directly from field locations (i.e., irrespective of AFOs). A one standard deviation increase in manure application field contributions was associated with 41% (p < 0.10) and 36% (inclusion lowers AIC) increases in the relative abundance of sul1 and intI1, respectively.

In sediment and surface water, the magnitude of the association between tet(W) and AFOs was greatest compared to other ARG responses. A greater association may indicate that more tet(W) genes are located at AFOs or may suggest that oxytetracycline is used frequently to prevent diseases at dairy AFOs compared to other antibiotics. Some studies suggest that tet(W) and tetracycline resistance may be more specific to dairy feces than other ARGs and antibiotic resistance phenotypes (Srinivasan et al., 2008; Kyselková et al., 2015). The strength of the associations from our study supports that tet(W) is correlated with dairy manure.

This is the first study to characterize ARG contributions with GORF or ORF SPMs for bovine sources. Our findings imply that overland and downstream transport and dilution from flow are key processes in disseminating AMR from AFOs and manure fields. This study also suggests that manure hauling to application sites is a factor in elevated ARGs in the environment. Our novel spatial predictors and modeling approach are in agreement with AFO's association with ARGs found by Pruden et al. (2012) using a different SPM and applied to a different geographical region. Pruden et al. (2002) found that the relative abundance of sul1 in sediment was positively correlated (R² = 0.35, p < 0.001) with the average upstream capacities of AFOs. However, these authors found no significant relationship with the relative abundance of tet(W) in sediment (Pruden et al., 2012), which could reflect differences in livestock and antibiotic use. Our study quantified seven novel associations between bovine sources and ARG responses. We have quantified the associations between bovine sources and levels of erm(B), tet(W), qnrA, sul1, and intI1 in surface water (i.e., five novel associations) and the association between bovine sources and levels of erm(B) and tet(W) in sediment (i.e., two novel associations), which can inform ecological studies of AMR and microbial risk assessment.

The find stage of FIT enables the exploration of many databases of collocated, spatially distributed sources. The result is the set of source locations and associated information that consistently represent sources to the response from a cross-validation approach. Here, we report the find stage results for bovine sources because they were the most consistently selected source from the test stage of FIT. Across ARG responses in surface water, FIT modeled contributions from AFOs with the county manure storage option and contributions from manure application fields with the crop rotation land cover option.

For ARG responses in sediment, FIT modeled contributions from AFOs with WPDES CAFO locations, but differences existed in CAFO weightings. The difference in weighting CAFOs by animal units for erm(B) and tet(W) and equal weighting for sul1 may reflect differences in the transport processes resulting in elevated relative abundance, independent of the number of animal units (e.g., some other selective pressure emanating from sources).

One key difference in FIT's database selection was that, to represent AFOs, the WPDES CAFOs option was selected for sediment responses and the county's manure storage option for surface water responses. One explanation for this difference is that sediment sampling was impossible at three sites during one high-flow event (Beattie et al., 2018). These manure storages are known to overflow during high-flow events (Burch et al., 2021), and the additional surface water samples may have better captured transient contamination from manure storages.

After implementing FIT across 5 ARGs in sediment and surface water, we found that three of the ARG responses (Table 2) are positively associated with land-applied waste sources in sediment (sul1 and intI1) or surface water [tet(W) and intI1]. ARG responses in sediment, municipal waste, or septage land application characterized land-applied residential waste sources. For sul1, FIT selected the ORF land-applied residential waste source term representing land-applied waste from residential origins (i.e., land application of municipal waste or septage consisting of solid or semi-solid residue generated during the treatment of domestic sewage via primary, secondary, or advanced wastewater treatment and the wastewater contents of septic or holding tanks, dosing chambers, grease interceptors, seepage beds/pits/trenches, privies, or portable restrooms (Wis. Admin, 2021). A one standard deviation increase in contributions from land-applied waste from residential use was associated with a 63% (p < 0.05) increase in the relative abundance of sul1. For intI1, FIT selected the GORF land-applied waste source term, representing land-applied residential waste weighted by the density of nearby septic systems. A one standard deviation increase in land-applied waste from residential use was associated with a 43% (p < 0.05) increase in the relative abundance of intI1.

For ARG responses in surface water, industrial waste's land application characterized these five ARGs' secondary sources. For tet(W), FIT selected the ORF land-applied industrial waste source term representing by-product solids from the animal product or food processing industry (i.e., remains of butchered animals, paunch manure, cheese production waste, and vegetable waste materials). A one standard deviation increase in land-applied waste from industrial use was associated with a 36% (p < 0.10) increase in the relative abundance of tet(W). Then, for intI1, FIT selected the ORF land-applied industrial waste source term [i.e., “both the by-product solids from the animal product or food processing and liquid waste such as silage, leachate, whey, whey permeate, whey filtrate, contact cooling water, cooling or boiler water containing water treatment additives, and wash water generated in industrial, commercial, and agricultural operations” (Wis. Admin, 2021)]. A one standard deviation increase in land-applied waste from industrial use was associated with a 41% (p < 0.10) increase in the relative abundance of intI1.

Our findings are consistent with current knowledge that ARGs are enriched in biosolids from treatment processes (i.e., primary, secondary, or advanced wastewater treatment) (Chen and Zhang, 2013; Burch et al., 2014; Pepper et al., 2018). This is the first study to report an association between modeled contributions from spatially distributed land-applied waste and ARGs recovered from riverbed sediment and surface water. This is also the first study to show that septage and municipal or industrial waste disposal on land pollute and correspond with a quantifiable environmental impact on antibiotic resistance levels in sediment and surface water. The WPDES database lists the facility names associated with the land application sites. After searching on company websites for the products associated with each facility producing industrial wastewater or sludge destined for land application, we found that 88.8% of the industrial land-applied waste sites are associated with dairy and meat products. None (i.e., 0%) of the facilities were associated with pharmaceuticals. This suggests that the disposal of industrial wastes from dairy and meat processing extends the polluting ability of industrial livestock agriculture and that industrially produced, land-applied pharmaceutical waste is not a source of elevated ARGs in this region.

The bottom of Table 2 shows the results from the Inform and Test stages corresponding to the soil as a source. For surface water responses, only soil type D (clay loam, silty clay loam, sandy clay, or silty clay) coverage was associated with a 70% increase in the log10 relative abundance of erm(B). In sediment responses, soil type A (i.e., sand, sandy loam, loamy sand, or gravel) coverage was associated with 152% and 79% increases in sul1 and intI1 log10 relative abundances, respectively.

Soil type D contains clay, and we expect ARGs to correlate with clay based on previous research (Mao et al., 2014; Wang et al., 2016). A previous study also suggests that clay microbial communities are more resilient to change from anthropogenic sources compared to other soil types, like sand (Neumann et al., 2013). Therefore, this contribution to erm(B) relative abundance of clay soil sources may suggest the long-term impacts of agricultural sources on clay microbial communities. However, we found different results from the surface water results. Overall, soil appears to significantly contribute to the relative abundance of these ARGs, and the differences in soil sources for surface water and sediment in this study have many potential explanations (e.g., adsorption, desorption, and absorption between overland soil, riverbed sediment, and surface water). More observational data and controlled mesocosm-scale experiments are needed to validate these findings and characterize these complex dynamics.

WWTPs were expected to show associations with sul1 and intI1, as sul1 is often conserved in integron-integrase mobile genetic elements (Pruden et al., 2012; Amos et al., 2015), but the FIT model selected neither municipal nor industrial WWTPs associating with any of the five relative abundance responses across surface water or sediment. One explanation is that in this rural area of ~20,000 people (US Census Bureau, 2022), only five of the WWTPs were flow-connected to the 20 sampling sites. Our study area consists of a larger bovine–human ratio (Borchardt et al., 2021; Burch et al., 2021; Wiesner-Friedman et al., 2021a) compared to the study areas of previous ARG research in the South Platte River Basin in Colorado, United States (Pruden et al., 2012), and the Thames Watershed in Oxfordshire, United Kingdom (Amos et al., 2015).

In our study, both sul1 and intI1 were associated with land-applied municipal waste and septage. The land-applied municipal waste and septage represent the aggregation of treated septage and wastewater, suggesting that some ARGs may originate from WWTPs, but most likely, the persistent application of biosolids on land represents a more significant source than WWTP effluent or septic systems in this region. Quantitation methods with low detection limits for intI1, more flow-connected sampling sites to sources, or a different ARG panel may be needed to detect the impacts of WWTP effluent or septic systems.

Bovine source terms associated with elevated ARGs in sediment had exponential influence ranges of α_O<13 km on the river network, indicating that decayed contributions would still be detected in the river when manure fields were up to 13 km away from the river network. For all bovine source terms associated with ARG responses in surface water, decayed contributions would still be detected when sources were up to 10 km from the river network. For all land-applied waste source terms associated with ARGs, the influence range was up to 8 km for sediment and 10 km for surface water responses.

These overland, exponential influence ranges, α_O, were determined from the inform stage of FIT. In a previous study (Wiesner-Friedman et al., 2021a), we remarked that this hyperparameter captures more than average overland transport. This hyperparameter characterizes the extent to which a microbial response can capture a signal from sources so that longer overland influence ranges may indicate either A) longer transport and B) an increased probability of detection, or both. Previously, we found longer overland influence ranges around sources for host-associated Bacteroides in sediment (i.e., α_O> 1 km) vs. surface water (i.e., α_O < 1 km) (Wiesner-Friedman et al., 2021a).

In this study, on average, we found long overland influence ranges [i.e., a mean value of αO^=7.07 km (95%CI: 4.22 km, 9.92 km)] for both sediment and surface water. ARGs can exist naturally in soils (D'Costa et al., 2011) and be carried by either aerobic or anaerobic bacteria (Xu et al., 2021); They can be transferred to other bacteria by several genetic mechanisms (Davies and Davies, 2010). These factors may increase the transport of and ability to detect ARGs in surface waters compared to anaerobic Bacteroides gene markers. Furthermore, Bacteroides persist briefly outside of the gut of their hosts (Ballesté and Blanch, 2010), so shorter overland influence ranges would be expected. One implication of longer influence ranges around bovine sources for ARGs compared to host-associated Bacteroides genes is that host-associated markers may underrepresent the risks associated with fecal contamination from bovine sources in surface water. Additionally, this research and previous studies have found an increased probability of detecting microbial genes in sediments compared to surface water (see Supplementary material S2) (Kasich et al., 2012; Wiesner-Friedman et al., 2021a), which suggests that risks associated with particular sources may be underrepresented from sampling transient surface water compared to time-integrated polluted sediment.

The FIT model of ARGs in the environment is the first to characterize overland influence ranges around these sources. A setback distance of 34–67 m from surface water has previously been recommended for manure and slurry land application under experimental conditions (Hall et al., 2020). The exponential influence ranges from our study indicate that sources very close to the river will greatly contribute but that distant sources up to 13 km away from the river network can also impact ARG levels. A factor that may lead to long influence ranges in this region is the karst geology, where fractures, sinkholes, caves, disappearing streams, and springs may provide direct pathways for contaminants, including antibiotics and ARGs, to reach ground and surface waters (Stange and Tiehm, 2020; Xiang et al., 2020). Our findings are consistent with a study in a karst region in Germany, where elevated ARGs and human-specific fecal markers were detected in a spring 9 km away from the suspected source (Stange and Tiehm, 2020).

Since this is the first study to report LUR results for erm(B) in sediment and erm(B), tet(W), qnrA, sul1, and intI1 in surface water, we identified the databases that most reliably report the spatial location of bovine and land-applied waste sources associated with elevated relative abundance of these ARGs and their corresponding regression coefficients (Table 2), and we have shown the spatial impact of these pollution sources (Figure 1).

Figure 1A shows the geographical distribution of an increase expected in the relative abundances of erm(B), tet(W), and sul1 in sediment associated with risk-conservative contributions from AFOs via the application of manure. This figure, which results from a GORF SPM, suggests that elevated ARGs in riverbed sediment are localized around manure application fields. However, the extent to which those sources (viz., differently sized black squares) qualify as polluters relates to their proximity to AFOs and the scale of the operation (viz. differently sized blue diamonds). This localized pollution may be influenced by ARG soil attachment (Barrios et al., 2021). In previous LUR, AFOs were associated with sul1 in sediment, but the spatial localization was not reported or depicted (Pruden et al., 2012). Our LUR/FIT modeling showing localized sediment ARG pollution is consistent with the field studies of the enrichment of ARGs in AFO manure application and dissemination into the environment (Fahrenfeld et al., 2014; Wallace et al., 2018).

Figure 1B shows the geographical distribution of an increase expected in the relative abundances of tet(W) and intI1 in surface water associated with risk-conservative contributions from land application sites. This figure, which results from an ORF SPM, suggests that land application sites have a dispersed impact on elevated ARGs in surface water.

These maps indicate locations where additional monitoring may be needed to understand the impacts of different sources on environmental and public (i.e., water users') health. The upper-bound hyperparameter values and database information also serve other regions with similar geography, land use, agricultural practices, and population to preliminarily define monitoring locations for AMR studies.