Predicting Nutrient Incontinence in the Anthropocene at Watershed Scales

Quantifying nutrient attenuation at watershed scales requires long-term water chemistry data, water discharge, and detailed nutrient input chronicles. Consequently, nutrient attenuation estimates are largely limited to long-term research areas or modeling studies, constraining understanding of the ecological characteristics controlling nutrient attenuation and complicating efforts to protect or restore water quality in developed and developing regions. Here, we combined long-term data and a broad suite of biogeochemical parameters from 49 watersheds in northwestern France to test how well instantaneous measurements can predict nitrogen (N) and phosphorus (P) attenuation at watershed scales. We evaluated 13 biogeochemical and 12 hydrological proxies of hydrological flowpaths, residence time, and biogeochemical transformation. Across the 49 watersheds, nutrient attenuation ranged from 88 to −2% for N and 99–96% for P. The strongest biogeochemical proxies of N attenuation were NO3- isotopes, rare earth elements (REEs), radon, and turbidity, together explaining 75% of observed variation. For P attenuation, REEs, NO3- isotopes, molecular weight of dissolved organic matter, and radon were the strongest proxies, but only explained 27% of observed variation. However, a single hydrological parameter—annual runoff—explained 91% of N attenuation and the relative abundance of schist bedrock explained 56% of P attenuation. We discuss how runoff both controls and reflects watershed hydrology, biogeochemistry, and nutrient attenuation. For example, runoff was correlated with long-term decreases in nutrient concentration, demonstrating how leakier watersheds recover more quickly from nutrient saturation. Given the immense fertilization capacity of modern society, we propose that eutrophication can only be solved by reducing nutrient inputs, though hydrochemical proxies can provide valuable information on where to carry out essential food production activities.

Quantifying nutrient attenuation at watershed scales requires long-term water chemistry data, water discharge, and detailed nutrient input chronicles. Consequently, nutrient attenuation estimates are largely limited to long-term research areas or modeling studies, constraining understanding of the ecological characteristics controlling nutrient attenuation and complicating efforts to protect or restore water quality in developed and developing regions. Here, we combined long-term data and a broad suite of biogeochemical parameters from 49 watersheds in northwestern France to test how well instantaneous measurements can predict nitrogen (N) and phosphorus (P) attenuation at watershed scales. We evaluated 13 biogeochemical and 12 hydrological proxies of hydrological flowpaths, residence time, and biogeochemical transformation. Across the 49 watersheds, nutrient attenuation ranged from 88 to −2% for N and 99-96% for P. The strongest biogeochemical proxies of N attenuation were NO − 3 isotopes, rare earth elements (REEs), radon, and turbidity, together explaining 75% of observed variation. For P attenuation, REEs, NO − 3 isotopes, molecular weight of dissolved organic matter, and radon were the strongest proxies, but only explained 27% of observed variation. However, a single hydrological parameter-annual runoff-explained 91% of N attenuation and the relative abundance of schist bedrock explained 56% of P attenuation. We discuss how runoff both controls and reflects watershed hydrology, biogeochemistry, and nutrient attenuation. For example, runoff was correlated with long-term decreases in nutrient concentration, demonstrating how leakier watersheds recover more quickly from nutrient saturation. Given the immense fertilization capacity of modern society, we propose that eutrophication can only be solved by reducing nutrient inputs, though hydrochemical proxies can provide valuable information on where to carry out essential food production activities.
Keywords: nitrogen, phosphorus, attenuation, biogeochemical proxy, nutrient legacy INTRODUCTION Since the Industrial Revolution, humans have more than doubled reactive nitrogen (N) inputs (Gruber and Galloway, 2008) and quadrupled phosphorus (P) inputs (Elser and Bennett, 2011) into the Earth's ecosystems (Seitzinger et al., 2010;Foley et al., 2011;Abbott et al., 2018a). Consequently, approximately 80% of freshwater and coastal ecosystems now experience eutrophication induced by anthropogenic inputs of N and P (Howarth et al., 2000;Galloway et al., 2003;Poisvert et al., 2017;Le Moal et al., 2019). Though recognized as a planetary priority (Foley et al., 2011;Steffen et al., 2015;Le Moal et al., 2019), efforts to reduce eutrophication have had mixed results, partly because of two challenges that emerge at watershed scales (i.e., 1-10,000 km 2 ). First, it is difficult to quantify the overall residence time of nutrients in complex watersheds, which ranges from minutes to millennia as nutrients may be recycled or stored in plant biomass, soil, and groundwater Sebilo et al., 2013;Marçais et al., 2018;Carey et al., 2019;Kolbe et al., 2019). Second, the capacity of ecosystems to remove or permanently retain nutrients via vertical processes such as denitrification or diagenesis is highly variable, and the socioecological drivers (e.g., watershed characteristics and agricultural practices) of nutrient removal are poorly understood Abbott et al., 2016;Dupas et al., 2018;Goyette et al., 2018;Jarvie et al., 2019).
There are three general fates for nutrients moving through the soils, riparian zones, surface waters, and aquifers of a watershed: (1) Retention (i.e., a long or short-term delay) by biological and physical processes, including nutrient uptake, sorption, or hydrological residence time (Covino et al., 2010;Sebilo et al., 2013;Van Meter et al., 2016;Dupas et al., 2017;Ehrhardt et al., 2019), (2) Vertical removal to the atmosphere or lithosphere, including denitrification, aeolian transport, or mineral precipitation with various metals (Groffman et al., 2006;Seitzinger et al., 2006;Pinay et al., 2018;Randall et al., 2019), and (3) Longitudinal export from the watershed via surface or subsurface flow (Burt and Pinay, 2005;Seitzinger et al., 2010;Abbott et al., 2018a). The reactivity and mobility of organic and inorganic nutrients depend on and influence biogeochemical conditions Bernhardt et al., 2017), meaning that the fate of carbon, N, and P can vary substantially through time (e.g., storm events or seasons) and in space (e.g., different watersheds or biomes) (Dupas et al., 2016;Moatar et al., 2017;Musolff et al., 2017;Minaudo et al., 2019). Consequently, watersheds with similar nutrient inputs often have completely different nutrient export regimes, particularly for headwater watersheds that make up most of the terrestrial-aquatic interface (Bishop et al., 2008;Abbott et al., 2018b;Helton et al., 2018;Wollheim et al., 2018). This variability in nutrient attenuation capacity is likely associated with differences in surface and subsurface characteristics, including differences in water storage capacity and residence times, abundance and activity of biotic nutrient sinks (e.g., plant or microbial assimilation, dissimilatory microbial metabolism), and abiotic factors (e.g., high sorption capacity in soils, mineral precipitation, presence of chemical reducers in bedrock) (Hansen et al., 2002;Aquilina et al., 2012Aquilina et al., , 2018Thomas and Abbott, 2018;Kolbe et al., 2019). All these factors are influenced by changes in hydrological connectivity, uneven distribution of reactants and organisms due to the coevolution of surface and subsurface characteristics, and the stochastic nature of human and natural disturbance (Hansen et al., 2000;Thomas et al., 2015;Covino, 2017;Moatar et al., 2017).
Despite advances in understanding nutrient dynamics, it remains difficult to quantify nutrient retention and removal on timescales matching hydrological and nutrient residence times (Vitousek, 2004;Pinay et al., 2015;Ehrhardt et al., 2019). Additionally, many nutrient attenuation studies focus on single-nutrient dynamics, despite evidence that eutrophication is caused by interactions among multiple nutrients and factors (Carpenter et al., 1998;Elser et al., 2007;Paerl et al., 2016;Hobbie et al., 2017;Le Moal et al., 2019). Consequently, there is no straightforward way to predict a watershed's sensitivity to high nutrient inputs, seriously limiting our ability to prevent eutrophication or improve water quality where it has already been degraded. Quantifying nutrient attenuation at the watershed scale involves lengthy data acquisition and costly infrastructure (Burt, 1994;Howden et al., 2010;Burt et al., 2011), which are not always available in developing nations where agriculture and urbanization are intensifying fastest (Seitzinger et al., 2010;FAO, 2017;Dupas et al., 2019b).
In this context of variability, two conceptual approaches for predicting nutrient attenuation are the hot spots/control points concept, which predicts where and when biogeochemical reactions are more likely to occur (McClain et al., 2003;Bernhardt et al., 2017) and the Damköhler number, which uses the ratio of residence time to reaction time to predict how much of a reactant can be transformed or retained (Ocampo et al., 2006;Zarnetske et al., 2011;Oldham et al., 2013). Pinay et al. (2015) and Abbott et al. (2016) proposed to combine these concepts in the HotDam framework by combining multiple proxies of biogeochemical transformation, hydrological flowpaths, and hydrological residence time. Combining or crossing multiple proxies such as solute and isotopic concentrations, hydrograph properties, and organic matter can illustrate terrestrial and aquatic conditions across nutrient flowpaths Abbott et al., 2016;Shogren et al., 2019), potentially paving the way for a more systematic understanding of what controls nutrient attenuation in watersheds.
Here, we apply the HotDam framework to identify controls on attenuation of N and P in 49 small to medium watersheds in Brittany, France using a crossed-proxy approach. Our overarching goals were to understand the drivers of nutrient attenuation in agricultural watersheds and test how well nutrient attenuation could be predicted with easily measurable proxies. These proxies included nutrient stoichiometry, organic matter biodegradability, dissolved gases, rare earth elements (REEs), nitrate isotopes, and hydrograph parameters (Figure 1). We hypothesized that nutrient attenuation would be controlled by the hydrological properties of the watersheds (e.g., water residence time and dominant flowpaths), surface and subsurface characteristics (e.g., land use, topography, geology), and biogeochemical conditions (e.g., stoichiometry and spatiotemporal distribution FIGURE 1 | Conceptual diagram of the crossed-proxy approach we used to investigate controls on watershed-scale nutrient attenuation. We hypothesized that watershed-scale nutrient attenuation was a function of water flowpath (red), residence time (purple), and biogeochemical reactivity (pink).
of electron donors and acceptors; see Table 1). We tested these hypotheses by calculating mass balances of N and P for each watershed using 10 years of local, regional, or national agency data, and then sampling the watersheds three times across flow conditions and seasons to analyze the broad suite of proxies.

Site Description and Experimental Design
The Brittany region of northwestern France has a rich repository of environmental data generated by academic and governmental research. For example, 27 of the 49 watersheds in this study had a nearby surface water monitoring station, which provided the concentration and discharge data necessary to calculate annual nutrient fluxes and long-term trends (Fovet et al., 2015;Thomas et al., 2019), and estimates of nutrient inputs were available for all watersheds (Poisvert et al., 2017). More generally, these intensively-managed, agricultural landscapes experience high but decreasing nutrient inputs, providing insight into how watershed-level nutrient fluxes respond to changes in nutrient loading Sutton et al., 2013;Poisvert et al., 2017;Dupas et al., 2018). The region is a part of the Armorican massif which is composed of metamorphic and igneous rock, primarily granite, schist, and micaschist (Aquilina et al., 2012;Goderniaux et al., 2013;Kolbe et al., 2016). The climate is temperate oceanic, with a mean annual temperature of 11.2 • C and mean annual rainfall ranging from 1,400 mm in the west to 600 mm in the east, relatively well-distributed throughout the year (Gascuel-Odoux et al., 2010;Thomas et al., 2015Thomas et al., , 2019. The area has an average stream density of about 1 km km −2 , relatively shallow groundwater, and hydromorphic riparian soils that cover about 20% of the land surface (Mourier et al., 2008;Dupas et al., 2013;Marçais et al., 2018). Land use is dominated by row crops, indoor pig and poultry husbandry, and pastureland for cows (a mean of 80% agricultural cover across the study watersheds; Table S1), making Brittany one of the highest density regions in France and Europe for animal breeding (Gascuel-Odoux et al., 2010;Poisvert et al., 2017;Kim et al., 2019). N and P concentrations in many Brittany watersheds are decreasing, attributable primarily to reduction of point sources such as wastewater and feedlot effluent Abbott et al., 2018b;Dupas et al., 2018).
We analyzed stream water chemistry in 49 agricultural watersheds ranging from 2.38 to 2,080 km 2 (Figure 2). Most of the watersheds were small to medium sized (mean area = 232 km 2 , median area = 18.2 km 2 ), which was a consequence of Brittany's geography as a peninsula and our design to capture variability in headwater watersheds that make up the majority of the land-water interface (Burt and Pinay, 2005;Bishop et al., 2008;Heathwaite, 2010). Sampling points were typically near the coast, but above the zone of tidal influence (Figure 2). To capture variable flow conditions and seasonal differences, we collected samples during field campaigns in November 2015 (lowest flow), March 2016 (highest flow), and June 2018 (moderate flow; Figure S1). During each field campaign, we  (Bertin and Bourg, 1994) SUVA 254 x x − + SUVA 254 is an estimate for aromaticity of organic molecules (Weishaar et al., 2003); higher SUVA 254 may correlate with lower bioavailability S R x x + − S R is inversely related to molecular weight of CDOM (Helms et al., 2008); higher S R correlates with lower molecular weight and may indicate more bioavailable compounds (Ewing et al., 2015) BDOC x x + − BDOC measures biodegradability of DOM (McDowell et al., 2006); higher BDOC could stimulate denitrification but also nutrient mineralization Stoichiometric ratios of C and N can indicate whether C or N attenuation is more likely (Sterner and Elser, 2002), and also indicate hydrologic flowpath because deep flowpaths have lower C:N ratios pH x ? ?
As pH increases, organic colloids become more electronegative and release adsorbed phosphate particles (Gu et al., 2019), and C and N are more available (Glass and Silverstein, 1998)  Proxies were chosen because they are indicative of hydrologic flowpath, residence time, and/or biogeochemical transformations. Predicted relationships between proxies and nutrient attenuation as the proxy value increases (e.g., higher concentration of DSI is hypothesized to increase N and P attenuation).
sampled all 49 sites within 1 week to capture watershed signals in comparable hydrological conditions. Seventeen of the watersheds were independent drainage basins and 32 were nested within the Couesnon and Rance watersheds (23 and 9 nested subwatersheds, respectively), allowing us to assess nutrient attenuation controls across a greater range of watershed sizes (Figure 2) and take advantage of previous research on those sites (Abbott et al., 2018a;Thomas et al., 2019). Though the sub-watersheds of the Couesnon and Rance watersheds are geographically close, they span the Brittanywide range of observed land use and watershed characteristics ( Table S1).

Conceptual Approach
Nutrient retention and removal (hereafter "attenuation") involve multiple hydrological and biogeochemical processes that are difficult to characterize because many of them are not directly observable due to long timescales or inaccessibility (e.g., groundwater processes) Kolbe et al., 2019), 4-dimensional variation in subsurface characteristics (Sebilo et al., 2013;Musolff et al., 2015), and nutrient legacies (Van Meter and Basu, 2017;Ehrhardt et al., 2019). Consequently, to identify the ecological drivers of nutrient attenuation at watershed scales, we selected tracers or proxies that could be associated with hydrological flowpath, residence time, and biogeochemical reactions Abbott et al., 2016). Informed by the ecological control points concept, which assesses reaction rates in a spatiotemporal context (McClain et al., 2003;Bernhardt et al., 2017); and the Damköhler approach which assesses overall attenuation capacity (Ocampo et al., 2006;Oldham et al., 2013), we attempted to quantify biogeochemical and hydrological controls of N and P attenuation at watershed scales with a crossed-proxy approach . We were particularly interested in why relationships between land use and stream nutrient concentrations and fluxes often break down at small scales (Burt and Pinay, 2005;Heathwaite, 2010), and how well watershed characteristics and easily measured proxies could predict nutrient attenuation and shine light on the relative importance of surface and subsurface attenuation processes (Ben Maamar et al., 2015;Dupas et al., 2019a;Kolbe et al., 2019) and hydrological time lags in soils, sediments, and aquifers (Thomas et al., 2012;Sebilo et al., 2013;Van Meter et al., 2016). To address these questions, we used a diverse set of physicochemical parameters described in detail in the Supplementary Information (SI: Proxy Toolbox) and briefly outlined below (Figure 1). set by the bedrock, but redox conditions can cause selective changes (Dia et al., 2000;Gruau et al., 2004). Specifically, cerium (Ce) readily oxidizes to Ce +4 and precipitates in the presence of oxygen (Moffett, 1990;De Carlo et al., 1997;Braun et al., 1998), creating a negative Ce anomaly in water that has experienced consistently oxidizing conditions (Gruau et al., 2004;Pinay et al., 2015). Because redox conditions and organic matter availability strongly influence N and P attenuation processes (e.g., denitrification; P adsorption by Fe-oxyhydroxides) (Stumm and Sulzberger, 1992;Pinay et al., 2015;Gu et al., 2019), we tested the relationship between the Ce anomaly of stream waters with nutrient attenuation at watershed scales. To assess water residence time, we used DSi concentration, which has been found to strongly correlate with subsurface residence time in many hydrogeological contexts (Ayraud et al., 2008;Marçais et al., 2018). Our DSi estimates of residence time, calculated using the empirical relationship derived in Marçais et al. (2018), agreed with estimates derived from chlorofluorocarbons (CFCs) and other dissolved gases from other studies in this region (Molénat et al., 2013;Ben Maamar et al., 2015;Kolbe et al., 2016). Radon-222 ( 222 Rn) is another tool to constrain groundwatersurface water interactions (Bertin and Bourg, 1994;Cable et al., 1996;Stieglitz et al., 2010). A product of natural radioactive decay in igneous bedrock, 222 Rn has a half-life of 3.82 days, making it an ideal tracer of deep flowpaths (Oyarzún et al., 2014). Because deep and long flowpaths increase the likelihood of encountering redox conditions suitable for N removal pathways such as denitrification, we measured 222 Rn concentration in all stream water samples.
We used multiple characteristics of DOM to assess the degree of biogeochemical processing, nutrient source, and multielemental interactions. DOM has been described as a master variable that influences multiple nutrient cycles (e.g., it is a major source of inorganic N and P in nutrient-poor ecosystems) and general physicochemical conditions (McDowell, 2003;Zarnetske et al., 2018). DOM consists of dissolved organic carbon (DOC), N, P, and other nutrients in molecular forms ranging from complex molecules to simple compounds. The molecular composition of the DOM influences its biodegradability and photoreactivity, affecting its persistence in the ecosystem and influence on nutrient cycles (Wymore et al., 2018;Harjung et al., 2019). In addition to biodegradability incubations (details below), we calculated two optical proxies of DOM composition: specific ultra-violet absorbance at 254 nm (SUVA 254 ) and the spectral ratio (S R ) of slopes within the 275-290 and 350-400 nm range (Weishaar et al., 2003;Helms et al., 2008;Vonk et al., 2015). DOM concentration and characteristics can also indicate hydrological flowpaths because DOM is less abundant and more microbially altered in groundwater (Shen et al., 2015;Mu et al., 2017;Coble et al., 2019).
Nutrient stoichiometry, which is based on conservation of mass and constant proportions in many organisms, allows prediction of retention or release of different compounds based on availability and relative demand (Sterner and Elser, 2002;Allen and Gillooly, 2009;Helton et al., 2015). We used nutrient ratios as metrics of flowpath and biogeochemical transformation. For example, a negative relationship between DOC and NO − could simply be caused by a negative correlation between sources, where watersheds that favor deeper hydrological flowpaths have a carbon-poor and nitrogen-rich signal (Abbott et al., 2018b). Therefore, we predicted higher nutrient attenuation in catchments with greater DOC:NO − 3 ratios.

Hydrograph Metrics
The characteristics of river flow can indicate fundamental hydrological properties at watershed scales (Fang and Shen, 2017;Moatar et al., 2017). Higher peak flows during storm events can indicate greater near-surface runoff, while less responsive hydrographs and higher base flows between events can indicate longer residence time and greater proportion of subsurface flow (Feijoó et al., 2018;Kirchner, 2019). In this context, we calculated several, non-redundant hydrological metrics (see Table S2) based on daily stream flow: (i) the mean, (ii) coefficient of variation, (iii) skewness, (iv) kurtosis, (v) the autoregressive lag-one correlation coefficient (AR1), (vi) the amplitude, (vii) the phase of the seasonal signal (Archfield et al., 2014), and (viii) the W2. The W2 is an index of hydrologic reactivity that is the percentage of annual discharge that occurs during the highest 2% of flows (Walsh and Lawler, 1981;Moatar et al., 2013Moatar et al., , 2017.

Field and Laboratory Analysis
To quantify the proxies described above, we collected water samples and measurements from 49 watersheds throughout Brittany, France (Figure 2). We selected the 49 sites based on accessibility, availability of historical data (nutrient input and export chronicles and land-use data), and to cover a range of watershed sizes. Field campaigns were in early November of 2015, late March of 2016, and late June of 2018. At each site, we collected one 5-L sample of stream water for immediate sensor readings and eight smaller samples for laboratory analyses. From the first sample, we used a handheld multiparameter probe (YSI, incorporated; Yellow Springs, USA) to measure dissolved oxygen, redox, pH, temperature, and turbidity. We determined SUVA 254 and S R from the same sample with a field-deployable spectrophotometer (s::can; Vienna, Austria). For the lab analyses, we immediately filtered subsamples using a 50 mL syringe and two 250 mL filter towers. We used a 0.2 µm cellulose acetate syringe filter to prepare samples for the analysis of cations, REEs, and NO − 3 isotopes. For the first filter tower, we used a 0.45 µm cellulose acetate filter to prepare samples for Molybdate-Reactive Phosphorus (MRP), anions, and DOC analysis. For the second tower, we used a 0.7 µm glass fiber filter, which removes most particulates but allows many bacteria to pass (Vonk et al., 2015), to prepare samples for the biodegradable DOC (BDOC) bioassay experiment. All filters, towers, and syringes were pre-rinsed with de-ionized water and flushed with sample prior to collecting final samples. We also collected an unfiltered, bubble-free 200 mL sample for 222 Rn analysis. These samples were analyzed for 222 Rn within 12 h of sampling using a radon detector (Durridge RAD7 analyzer, Billerica, USA), with most samples analyzed immediately in the field. Delayed samples were adjusted for time lags to correct 222 Rn decay. The 222 Rn values for the spring sampling were lost due to operator error.
MRP concentrations were determined colorimetrically via reaction with ammonium molybdate (Murphy and Riley, 1962), with a precision of ±4 µg l −1 (Gu et al., 2018). Nitrate isotope samples were frozen immediately and shipped to the UC Davis Stable Isotope Facility for analysis of δ 15 N and δ 18 O of NO − 3 by bacterial denitrification assay (McIlvin and Casciotti, 2011). Isotope ratios of δ 15 N and δ 18 O were measured using a ThermoFinnigan GasBench + PreCon trace gas concentration system connected to a ThermoScientific Delta V Plus isotope-ratio mass spectrometer (Bremen, Germany) with a precision of ± 0.4% and 0.5% for δ 15 N and δ 18 O, respectively. Cations and REE samples were analyzed by inductively coupled plasma mass spectrometry (ICP-MS; Agilent 7700×, Santa Clara, USA). Calibration curves and accuracy controls were performed following best practices (Yeghicheyan et al., 2013), using river water reference material for trace elements with a wide compositional range (SLRS-5, National Research Council of Canada). De-ionized water purified with a Milli-Q (Millipore, Darmstadt, Germany) system was used for blanks. Total relative uncertainties were ±5%.
We conduced BDOC bioassay experiments using four replicates of 100 mL aliquots of filtered stream water for each site and incubated in a dark incubation chamber for 28 days at 20 • C (Vonk et al., 2015). We sampled each replicate at the beginning (t 0 ) and end of the incubation (t 28 ). The t 0 samples were acidified with 6M HCl to a pH of 2 and stored in the refrigerator at 4 • C until the t 28 sampling. After the t 28 samples were acidified, the t 28 and t 0 samples were analyzed together for DOC (Shimadzu TOC-5050A, Kyoto, Japan, precision ±5%) within 1 week. We calculated the percent BDOC for each sample using the following equation: where DOC t0 and DOC t28 are the concentrations of DOC at t 0 and t 28 . We reported the mean of the 4 replicates as the overall BDOC value per sampling for each site.

Watershed Characteristics and Nutrient Trends
To test how climate and landscape characteristics affect nutrient attenuation, we delineated all watersheds using ArcMap (ESRI) and extracted landscape characteristics including vegetation cover, land use, bedrock type, river network density, and flow regulation. Mean annual temperature and precipitation were downloaded from the WorldClim database (Fick and Hijmans, 2017) using the "raster" package in R (Version 3.5.2; R Core Team, 2018). To obtain long-term nutrient and hydrological discharge data, we selected monitoring stations using the "Near" tool in ArcMap. We obtained long-term data for 27 of the 49 watersheds and calculated annual runoff, nutrient flux, and attenuation for those 27 watersheds using 10 years of flow and water chemistry data (NO − 3 and TP) from July 2008-July 2018. We calculated the eight flow metrics described above from the flow time series using the EflowStats R package.
We calculated mean annual specific discharge by averaging the flow per year and dividing by watershed area. We calculated N and P fluxes exported from the watershed using the dischargeweighted concentration method Raymond et al., 2013): where Ci and Qi represent concentration and runoff at the time of sampling, Q is mean annual runoff, and k is a conversion factor to obtain fluxes in kg N or P ha −1 yr −1 . We calculated apparent attenuation of NO − 3 and TP using the following massbalance equation: where Flux is calculated as described above (Equation 2) and Surplus is the watershed-specific nutrient surplus based on fertilizer inputs and crop yields (Poisvert et al., 2017;Dupas et al., 2018). Because not all watersheds had nearby water quality and discharge stations, we calculated another metric of nutrient attenuation: the residuals from a regression of NO − 3 or TP concentration against percent agricultural cover. While this metric is much coarser than the mass-balance estimates of attenuation (Equation 3), it provided an independent metric of attenuation that was quantifiable even for watersheds without discharge and long-term data. Additionally, it allowed calculation of seasonal dynamics in the relationship between land use and nutrient concentration, whereas the mass-balance-based nutrient attenuation metric provided a single value derived from the entire decadal time series. We reasoned that because most nutrient inputs in this region come from agricultural activity, watersheds with nutrient concentrations above the trend line would be less attenuative and points below the line would be more attenuative, comparatively to the average data set (Figure 3). We used NO − 3 concentration for these calculations because it accounted for 92% of measured total N, and TP to account for all measured P ( Figure S2). One weakness of both the mass-balance and residual metrics of nutrient attenuation is that apparent attenuation can occur if significant time lags exits between nutrient inputs and export (Basu et al., 2011;Ehrhardt et al., 2019). Over the past two decades, agricultural practices have improved and many nutrient point sources have been eliminated in western France (Poisvert et al., 2017;Abbott et al., 2018b;Dupas et al., 2018), and we reasoned that the slope of the decline in nutrient concentration would be negatively correlated with hydrological nutrient legacy (i.e., watersheds with more groundwater and soil nutrient storage would have slower rates of nutrient decrease after inputs ceased or were reduced). For the 27 watersheds with long-term data, we calculated Theil-Sen slope estimates using the Siegel method, which provides an estimate of change through time that is highly robust to extreme values due to real variability (e.g., during high or low flows) or analytical errors (Abbott et al., 2018b). We correlated these slope estimates with the nutrient attenuation metrics to assess how much of the apparent attenuation could be due to hydrological nutrient legacy rather than more permanent removal and retention processes.

Quantifying Spatial Stability and Identifying Drivers of Nutrient Attenuation
We determined the persistence of spatial patterns through time for proxies and water chemistry across watersheds using the spatial stability concept (Abbott et al., 2018a;Dupas et al., 2019b). The Spearman rank correlation between each pair of sampling dates for each parameter indicates how much of the spatial structure is preserved through time. Spatial stability of water chemistry can result from high spatial variability and synchronous temporal variability (Abbott et al., 2018a;Dupas et al., 2019b).
To disentangle the controls on nutrient attenuation, we regressed each proxy against nutrient attenuation (mass balance and residuals) and fluxes. We used pairwise correlation and multiple linear regression (MLR) to identify which proxies were the most closely associated with nutrient attenuation. Because many of the proxy data were not linearly related with attenuation and fluxes, we used Spearman correlations to quantify the strength of relationships. For the MLR models we grouped predictors into two categories: biogeochemical proxies and hydrological and watershed characteristics (hereafter referred to as the proxy and hydro models, respectively). To test our hypotheses about controls on nutrient attenuation, we categorized each predictor based on whether it was most associated with hydrological and watershed characteristics or biogeochemical reactions (Figure 1, Table 1). We standardized all predictors (mean = 0 and standard deviation = 1) to allow comparison of parameter coefficients as a measure of relative contribution to model prediction of the response variable (nutrient mass balance estimates), and we checked for FIGURE 3 | Visualization of how we calculated residuals in the relationship between agricultural land use and nutrient concentrations for each catchment. Because agriculture is the predominant nutrient source in the Brittany region, we used departure from the relationship between agriculture and nutrient concentration as a metric of nutrient attenuation capacity. Points above the regression line represent less attenuation capacity because there is more nutrient in the system than would be expected with this rough estimate of nutrient inputs.
Frontiers in Environmental Science | www.frontiersin.org multicollinearity among predictors using the variance inflation factor (Abbott and Jones, 2015). We ran the models for each category (i.e., proxy and hydro models) and selected the most parsimonious model for each by stepwise regression. Though model selection techniques continue to be controversial in ecology (reviewed in Malone et al., 2018), our primary goal was to assess relative influence of the major predictors. As such, we simply discuss the overall trends (e.g., what parameters appeared repeatedly in multiple models) of individual and multiple regression results (Tables 2, 3) and we abstain from interpreting the inclusion or exclusion of less influential predictors.

Nutrient Context and Spatial Stability of Water Chemistry
The three samplings captured distinct hydrological conditions, with low but variable discharge among sites in the fall sampling (November of 2015), highest discharge in the spring (March of 2016), and low and consistent discharge among sites in the summer (June of 2018; Figure S1). Nutrient concentrations and stoichiometry varied substantially across the watersheds and through time (e.g., 0-15 mg L −1 of N-NO − 3 and 0-0.6 mg L −1 of MRP). We observed a negative relationship between DOC and NO − 3 and a positive relationship between DOC and MRP, which varied somewhat seasonally (Figure 2b, Figure S1). Inorganic forms of N and P dominated total concentrations for these nutrients across watersheds, with NO − 3 making up 93% of total N and MRP making up 74% of TP ( Figure S2).
For the 27 watersheds with decadal nutrient concentration data on a monthly time step, the long-term trends were fundamentally different for C, N, and P (Figure 4). Theil-Sen slope estimates for DOC concentration indicated slight increases through time for 75% of the watersheds, decreases for NO − 3 concentration for all but three watersheds, and no change for PO 3− 4 and TP concentrations (Figure 4).
Spatial stability (persistence of spatial rankings through time) varied substantially among proxies and solutes ( Figure 5). Conductivity, δ 15 N, TP, DSi, and many major ions showed high spatial stability (i.e., more than half the spatial pattern among the three sampling dates was preserved; ρ > 0.7). All the DOM properties, 222 Rn, O 2 , temperature, and NO − 2 showed very low stability (i.e., ρ < 0.4) indicating substantial seasonal and potentially interannual variability, with the rest of the parameters showing moderate stability (Figure 5).

Differences in Attenuation and Fluxes in Agricultural Watersheds
Median nutrient attenuation as calculated by mass balance was 58.1% for N and 98.6% for P across the watersheds. N and P mass balance results were weakly correlated (R 2 = 0.18, p < 0.001), indicating that watersheds with high attenuation for N tended to also have high attenuation for P (Figure 6). Results for N mass balance showed substantially higher variability than P mass balance, ranging from −2 to 88.4% for N but only 96.1 to 99.5% for P. As another estimate of nutrient attenuation, we calculated the residuals of the regression of percent agricultural cover and measured NO − 3 and TP concentrations (Figure 3). N-NO − 3 was significantly correlated with agricultural cover (R 2 = 0.21, p < 0.001) and residuals from that relationship ranged from −7.23 to 6.93 mg L −1 with a median of −0.37 mg L −1 . TP was not significantly correlated with agricultural cover (R 2 = 0.001, p = 0.70) and residuals from that non-significant relationship ranged from −0.54 to 0.047 mg L −1 with a median of 0.016 mg L −1 . N residuals showed a weak, positive correlation with N mass balance estimates (R 2 = 0.15, p < 0.001; Figure 6), but we observed no relationship between P residuals and P mass balance (Figure 6).
N fluxes varied by an order of magnitude (6.1-64 kg N ha −1 yr −1 ) and were negatively correlated with mass balance estimates (R 2 = 0.73, p < 0.001), with fluxes increasing as mass balance estimates decreased (Figure 6). P fluxes varied by a factor of  We considered models with 1-6 predictors and bolded the most parsimonious model in each category based on AIC.
Frontiers in Environmental Science | www.frontiersin.org  We considered models with 1-6 predictors and bolded the most parsimonious model in each category based on AIC. All predictors were standardized (mean = 0 and standard deviation = 1) to allow comparison of parameter coefficients as a measure of relative contribution to model prediction of the response variable (greater absolute value means more influential). Full models contained all variables within a category (e.g., measured proxies or watershed characteristics calculated or extracted from GIS) and the most parsimonious model (emphasized in bold) was chosen using the Akaike information criterion (AIC).

FIGURE 4 |
The distribution of Theil-Sen nutrient slopes for C, N, and P for the 27 watersheds with long-term nutrient data. Horizontal gray lines indicate zero slope (no temporal trend).

Individual Predictors of Nutrient Attenuation
Based on the pairwise Spearman correlations with hydrological parameters, annual runoff was strongly negatively correlated with the N and P mass balance estimates (ρ = −0.89 and −0.41, respectively), and the hydrological reactivity index (W2) was strongly positively correlated with N and P mass balance (ρ = 0.60 and 0.32, respectively). For biogeochemical proxies, N mass balance estimates was correlated with δ 15 N and δ 18 O, 222 Rn, Ce/Ce * , and DOC:NO − 3 stoichiometry (ρ > |0.40|) and had a weaker relationship with SUVA 254 (ρ = −0.27; Figure 7). N residuals showed similar results for δ 15 N and δ 18 O (ρ > 0.40) but had a much weaker relationship with 222 Rn, a much stronger relationship with DOC:NO − 3 (ρ = −0.20 and 0.71, respectively), and no significant relationship with Ce/Ce * or SUVA 254 . N fluxes showed similar but opposite relationships as mass balance, and a weak relationship with DSi was also observed (ρ = 0.23). P FIGURE 5 | Spatial stability of biogeochemical tracers and chemical concentrations across all 49 watersheds. Parameters are grouped thematically by color. Solutes with a Spearman's ρ above the horizontal line maintain more than half the spatial pattern for that parameter in pairwise comparisons among the three sampling dates. mass balance estimates was positively correlated with δ 15 N and δ 18 O (ρ = 0.41 and 0.29, respectively; Figure 8) and negatively correlated with 222 Rn (ρ = −0.35). P residuals were strongly negatively correlated with δ 15 N and 222 Rn (ρ = −0.52 and −0.44, respectively) and had weaker relationships with δ 18 O, Ce/Ce * , and S R (ρ < |0.25|). The P fluxes had similar but opposite relationships with the same proxies as mass balance, in addition to a negative relationship with DOC:NO − 3 (ρ = −0.40). To test how nutrient legacy could be related to apparent nutrient attenuation, we correlated the decadal Theil-Sen concentration slopes (Figure 4, Figure S3) with estimates of nutrient attenuation. N mass balance estimates and residuals were positively correlated with decadal trends for NO − 3 , indicating that the watersheds experiencing the slowest decreases in NO − 3 through time tended to have higher apparent attenuation (R 2 = 0.17 and 0.10, respectively, p < 0.05; Figure S3). N flux was negatively correlated with decadal NO − 3 trends, indicating that watersheds with faster decreases in NO − 3 had higher N fluxes (R 2 = 0.25, p < 0.05). PO 3-4 slopes were positively correlated with P mass balance estimates (R 2 = 0.15, p < 0.05), indicating that the watersheds with the slowest decreases or greatest increases in PO 3-4 had higher apparent P attenuation, but PO 3-4 slopes were not significantly correlated with P residuals. TP fluxes were negatively correlated with NO − 3 trends (R 2 = 0.27, p < 0.05).

Multiple Linear Regression Models of Nutrient Attenuation
For the hydro MLR models (only hydrological and watershed characteristics parameters), the models with only one predictor explained 91% of the variance explained by annual runoff for N mass balance estimates and 56% of the variance explained by the abundance of schist bedrock for P mass balance (see Tables 2,  3). The most parsimonious (based on AIC) hydro model for N mass balance estimates also included other hydrological and land use variables (e.g., mean annual flow and amplitude, and % artificial land cover, wetlands, and stream length), which slightly enhanced the model performance (adjusted R 2 = 0.94, AIC = 12.19). For P mass balance estimates, the initial hydro model (i.e., only relative abundance of schist) was improved by adding river density (adjusted R 2 = 0.62, AIC = 4.84). The final hydro model results are shown in Figure 9.
The proxy MLR models that used biogeochemical proxies to predict attenuation did not explain as much variance as the hydro models, but they indicated other processes besides runoff and bedrock type that are important drivers of attenuation. For example, the single most important proxy for N was δ 15 N, which accounted for 37% of the variation in N mass balance estimates. However, the most parsimonious model included 5 other proxies (ΣREE, 222 Rn, Ce/Ce * , δ 18 O, and FIGURE 6 | Relationships among our metrics of nutrient attenuation. Here, "Mass Balance" refers to the mass balance estimates of attenuation, "Flux" refers to the fluxes of those nutrients from the watersheds based on concentration and discharge data, and "Residuals" refers to the remaining variation in nutrient concentration after accounting for differences in agricultural land cover (Figure 3). Trend lines for significant correlations (p < 0.05) are shown and R 2 -values are reported. turbidity), together accounting for 75% of the variation in N mass balance estimates. For P mass balance estimates, the single most important proxy was also δ 15 N, accounting for 9% of the variation. However, the most parsimonious model included 4 other proxies (ΣREE, S R , 222 Rn, and Ce/Ce * ), together explaining 27% of the variation in P mass balance estimates. Furthermore, the shared proxies between N and P attenuation models (e.g., δ 15 N, ΣREE, 222 Rn, and Ce/Ce * ) had surprisingly similar relationships in both direction and magnitude, apart from Ce/Ce * which had a positive relationship with N mass balance estimates and a negative relationship with P mass balance estimates. Proxy model results are shown in Figure 9.

DISCUSSION
We hypothesized that hydrological properties, surface and subsurface characteristics, and biogeochemical conditions would interact to determine nutrient attenuation at watershed scales. Based on our analysis of 49 watersheds, 27 of which had longterm nutrient and discharge estimates, hydrological properties set the initial attenuation capacity, with secondary effects from biogeochemical conditions and land-use parameters. These results corroborate findings from other regions where runoff strongly controls nutrient attenuation and flux (Covino et al., 2010;Zarnetske et al., 2018;Ehrhardt et al., 2019). However, we point out that runoff is a high-level parameter that interacts FIGURE 7 | Correlation analysis of N attenuation metrics and biogeochemical proxies. Trend lines are shown when p < 0.05 for Spearman rank correlations between the parameters (relationships are not necessarily linear). This non-parametric test is robust to non-linearity in the relationships, which are only depicted linearly to indicate direction of the relationship. with other hydrological metrics, watershed characteristics, and biogeochemical reactions at multiple spatiotemporal scales. Below, we discuss how runoff influences and is influenced by many ecological dynamics including long-term nutrient legacies, redox conditions, and land use, and how this could inform our local and global efforts to solve eutrophication. FIGURE 9 | N and P attenuation (based on mass balance estimates reported as percentages) predicted by hydrological and watershed characteristics, and biogeochemical proxies. The "hydro" model included hydrological and watershed characteristics and the "proxy" model included biogeochemical proxy data from the three field campaigns (see Table 2).

Tradeoffs Between Nutrient Attenuation and Nutrient Recovery
At larger spatiotemporal scales, untangling proximate and ultimate drivers of nutrient flux and attenuation is exceedingly difficult because socioecological systems co-evolve based on shared and dynamic conditions (Thomas et al., 2015;Bogaart et al., 2016;Malone et al., 2018). Consequently, at medium to large scales, many risk factors for nutrient flux are co-linear, including climate, land use, soil type, ecosystem stature, and flow regime (Knoben et al., 2018;Lin et al., 2019;Smits et al., 2019). This has been observed within the boundaries of our study region, where agriculture tends to be more intense in areas underlain with micaschist, which is better suited for mechanical agriculture, but also more prone to nutrient export because of its soil and topographic properties (Thomas et al., 2015).
Of all the hydrochemical parameters we measured, annual runoff had the strongest single relationship with our long-term nutrient attenuation metrics (mass-balance, residuals, and nutrient flux). While this correlation is not surprising for flux, it reinforces a growing understanding that the routing and amount of lateral water flow fundamentally regulates hydrological connectivity and nutrient transport Zarnetske et al., 2018). Perhaps more interestingly, runoff was correlated with long-term trends in nutrient concentration (i.e., rate of nutrient recovery) in our dataset. This means that the high-runoff watersheds that have low nutrient attenuation tend to recover faster when nutrient inputs are decreased. This suggests that hydrological attenuation (Ehrhardt et al., 2019) is an important contributor to nutrient attenuation in these watersheds. It also makes sense in the context of soil and groundwater nutrient legacies-i.e., a system with more flow can be flushed faster-depending on the ratio between storage and flow (Basu et al., 2011;Van Meter and Basu, 2017;Abbott et al., 2018b). However, this inverse relationship between attenuation and recovery in agricultural watersheds raises the question of whether there is an inevitable tradeoff between nutrient resistance (i.e., nutrient removal capacity) and nutrient resilience (i.e., fast recovery rates; Goyette et al., 2018;Dupas et al., 2019a). Thinking comparatively between the two main nutrients in our study, the higher P attenuation values, which we know are associated with accumulation in soil and sediment rather than permanent removal (Hansen et al., 2002;Sharpley et al., 2013), and the much lower ratio of stock to flux (i.e., watershed P >> annual export of P) suggests that P follows this prediction of high resistance and extremely low resilience (Goyette et al., 2018;Haas et al., 2019).

Subsurface Processes as Key Regulators of Surface Concentrations and Fluxes
Though one valid interpretation of our findings is that hydrological dynamics, especially flow and river network density, dominate nutrient attenuation capacity, several lines of evidence suggest that a more complex range of factors regulates nutrient behavior in these watersheds. Water residence time, a fundamental ecohydrological property (Zarnetske et al., 2011;Kolbe et al., 2016;Thomas et al., 2016) was not associated with our attenuation metrics nor with recovery rates based on the long-term trends. This is likely because only a portion of the total residence time has the requisite conditions for biogeochemical processes to retain or remove nutrients. This concept, termed exposure time (Murphy et al., 1997;Oldham et al., 2013;Pinay et al., 2015), is central to the HotDam framework, because storage or transport in nonreactive zones only affects hydrological attenuation (i.e., time lags; Ehrhardt et al., 2019), not biogeochemical attenuation or removal . For example, groundwater denitrification increases with depth in ∼80% of watersheds where it has been quantified because of more abundant electron donors in deeper, unweathered bedrock . This dynamic also relates to the previously discussed relationship with runoff, because weathering rates are higher when there is greater water flow through the watershed (Marçais et al., 2018). Greater weathering could depress the vertical horizon where electron donors can sustain denitrification , meaning that higher runoff watersheds not only have shorter residence times, but that their biogeochemical reaction capacity may be lower.
Another indicator that annual runoff is not the sole control of nutrient dynamics is that several biogeochemical proxies associated with deep flowpaths and groundwater were strongly related to both N and P attenuation metrics. Nitrate isotopes, 222 Rn, REEs, and DOC:NO − 3 were influential proxies of attenuation and fluxes, pointing to the importance of deep flowpaths and anoxia as drivers of nutrient removal (i.e., denitrification) and attenuation via biological and abiotic processes (Gu et al., 2018(Gu et al., , 2019Kolbe et al., 2019). While P is typically not transported through deep groundwater flowpaths because of physicochemical properties (Hansen et al., 2002), P varies systematically with discharge in many catchments with total P often increasing and MRP often decreasing as flow increases . Additionally, excess NO − 3 in groundwater can trigger sulfate release and iron oxidation in the presence of pyrite, potentially mobilizing P via links with iron and sulfur in some environments (Smolders et al., 2010;Tang et al., 2016;van Dijk et al., 2019). In watersheds that are discharging their nutrient legacies (i.e., current nutrient loading < historical nutrient loading), NO − 3 concentrations are often higher in groundwater than in near-surface water (Abbott et al., 2018b;Dupas et al., 2018), which could create an indirect link between deep water flow and P attenuation dynamics (Dupas et al., 2015).
A surprising interaction between hydrology and biogeochemistry was the positive relationship between the W2 index (percentage of cumulative discharge that occurs during the highest 2% of daily discharge values) and nutrient attenuation. Contrary to major predictions of ecosystem ecology that nutrient attenuation decreases with hydrological pulses or floods (Fisher et al., 1998;Raymond et al., 2016;Wollheim et al., 2017), the watersheds in this region with higher W2 were relatively more attenuative of nutrients. However, it is important to note that Brittany has temperate hydrology with few floods (Thomas et al., 2019), meaning that the watersheds with relatively higher W2 are still not particularly flashy or hydrologically reactive compared with other regions . In this instance, the positive correlation between nutrient attenuation and the W2 index could be associated with the weathering mechanism described previously (i.e., watersheds with a larger proportion of surface flow vs. groundwater flow could have less weathered and more reactive aquifers; Kolbe et al., 2019) or it could be associated with sediment legacies in and near the river network. In this region and many others, large stocks of nutrient-laden sediments have accumulated in streams, riparian zones, and small reservoirs (Song and Burgin, 2017;Feijoó et al., 2018). These sediments can be important or even primary sources of nutrients, particularly during low flow periods (Dupas et al., 2015;Gu et al., 2018). Watersheds with more powerful or more frequent floods could have flushed out these sediments, effectively increasing net nutrient attenuation by reducing internal loading.
Together, these multi-proxy findings suggest that subsurface characteristics, including hydrological flowpaths and the location of biogeochemical activity, are fundamental to regulating nutrient attenuation and export at watershed scales.

Comparison With Nutrient Attenuation in Different Biomes and Land Use Regimes
The high level of variability in N attenuation that we observed has also been observed in different climatic and anthropogenic contexts. Watersheds in the Northeastern United States have NO − 3 mass balance estimates that range from 9 to 74% (Campbell et al., 2004). While point-source P loading is relatively wellconstrained (Kronvang et al., 2007;Grizzetti et al., 2012), watershed-scale estimates of P mass balance remain less common (Withers and Jarvie, 2008;Sharpley et al., 2013;Goyette et al., 2018). In natural ecosystems, P attenuation is often very high, though this varies strongly with ecosystem age and disturbance regime as well as stoichiometric conditions (Vitousek and Reiners, 1975;Verry and Timmons, 1982;Vitousek, 2004;Elser et al., 2007). In contrast, urban watersheds are often extremely leaky to P, with P mass balance estimates ranging from −7 to 74% with a mean of 22% (Hobbie et al., 2017). These urban P losses are attributable to the prevalence of impervious surfaces that decrease residence time, increase flashiness of water flow, and increase fluctuations between oxic and anoxic conditions (Hale et al., 2015(Hale et al., , 2016Blaszczak et al., 2019).

Preventing Rather Than Curing Nutrient Incontinence
The Anthropocene is characterized by concurrent and connected socioecological crises caused by human interference with many of the Earth's biological and abiotic cycles (Vitousek et al., 1997;Steffen et al., 2015;Abbott et al., 2019;Lade et al., 2019). Many proposed solutions to these crises seek to manage ecosystem response rather than modify human activity, essentially treating human demand as immutable (Jaggard et al., 2010;Garnier et al., 2014;Abbott et al., 2019). This tendency to cure rather than prevent is particularly prevalent in the fight against eutrophication, where researchers, policymakers, and land managers are going to great lengths to retain or remove anthropogenic nutrients in virtually every component of the watershed (Sharpley et al., 2008;Pu et al., 2014;Bol et al., 2018;Wollheim et al., 2018). Given how critical food production is to human wellbeing (Foley et al., 2011;Sutton et al., 2013;Rasul, 2016), these efforts are certainly justified and some nutrient mitigation practices such as cover crops, no-till cultivation, hedgerows, and riparian concentration can improve nutrient attenuation and enhance other ecosystem benefits (Hansen et al., 2000;Roley et al., 2016;Pinay et al., 2018;Thomas and Abbott, 2018). However, we believe that managing ecosystems to enhance nutrient removal will never solve eutrophication alone. Instead, a growing body of research suggests that managing human nutrient demand-i.e., reducing nutrient inputs into the Earth's watersheds-while also protecting nutrient attenuation capacity across the land to water gradient is the only comprehensive solution to eutrophication Bol et al., 2018;Leaf, 2018;Dupas et al., 2019a;Le Moal et al., 2019).
One reason for our position that nutrient prevention should be prioritized over attenuation is simply the immense fertilization capacity of our global society (Gruber and Galloway, 2008;Seitzinger et al., 2010). The capacity of ecosystems to assimilate or remove N and P can easily be overwhelmed by human nutrient inputs, and recovery from nutrient saturation can take decades to millennia Goyette et al., 2018;Haas et al., 2019;Randall et al., 2019). This raises the question of how can nutrient inputs be decreased while also achieving the sustainable development goal of eliminating malnutrition? One solution would be to only use N and P fertilizers to feed humans (Sutton et al., 2013). While population is projected to increase through 2050, most of the projected increases in global nutrient use stems from shifts in diet and per-capita consumption rather than population growth per se (Seitzinger et al., 2010;Liu et al., 2012;Sutton et al., 2013;Godfray et al., 2018). More than half of the global human nutrient and water footprints are due to livestock (Pelletier and Tyedmers, 2010;Herrero et al., 2013;Mekonnen and Hoekstra, 2015;Abbott et al., 2019), even though animal products represent a small portion of global nutrition and actually have a net negative effect on human health in many developed countries (Lassaletta et al., 2014;Godfray et al., 2018;Willett et al., 2019). This means that decreasing meat and dairy consumption could substantially reduce nutrient inputs while improving human health (Sutton et al., 2013;de Vrese et al., 2018;Godfray et al., 2018;Willett et al., 2019). Additionally, the creation of biofuels is incredibly nutrient intensive, though only marginally advantageous energetically and economically (Dominguez-Faus et al., 2009;Yang et al., 2011). Reducing livestock husbandry and biofuel cultivation-i.e., limiting fertilizer use to crops destined for direct human consumption-could decrease the global nutrient footprint by more than half, while also alleviating pressure on wildlife habitat and other Earth systems (Steffen et al., 2015;Springmann et al., 2018;Abbott et al., 2019).
While large-scale changes in diet and energy production have social and technical challenges of their own (Godfray et al., 2018;Willett et al., 2019), they could substantially reduce eutrophication, especially when combined with expanded control of nutrient point sources (e.g., wastewater treatment plants), which has been extremely effective at reducing nutrient export (Musolff et al., 2015;Abbott et al., 2018b;Le Moal et al., 2019). Assessing nutrient attenuation capacity at watershed scales with hydrochemical proxies could also contribute to these efforts by informing redistribution of essential human activities to more resistant or resilient parts of the landscape (Thomas et al., 2015;Abbott et al., 2018a;Dupas et al., 2019b;Refsgaard et al., 2019).

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
BA, GP, FM, and AL conceived the sampling design. RF, BA, SG, TK, and RD carried out field and laboratory work. RF, BA, and RD developed the core ideas of the paper with substantial contributions from all authors.

ACKNOWLEDGMENTS
We thank François Rouault at Agrocampus Ouest for his field assistance and expertise, and we thank Gabriel Cano at BYU for his guidance on model selection. We also thank Nathalie Le Bris and Frédérique Pallois for assistance with incubations and chemical analyses at the Plateforme Analytique ECOBIO, and Martine Bouhnik-Le Coz and Eliot Chatton for assistance with rare earth analyses at Géochimie Environnementale (GrEEN), in Rennes, France.