Influence of foliar traits, watershed physiography, and nutrient subsidies on stream water quality in the upper midwestern United States

The relationship between nutrient cycling and water quality in mixed-use ecosystems is driven by interactions among biotic and abiotic processes. However, the underlying processes cannot always be directly observed or modeled at broad spatial scales. Numerous empirical studies have employed land use patterns, variations in watershed physiography or disturbance regimes to characterize nutrient export from mixed-use watersheds, but simultaneously disentangling the effects of such factors remains challenging and few models directly incorporate vegetation biochemistry. Here we use structural equation models (SEMs) to assess the relative influence of foliar chemical traits (derived from imaging spectroscopy), watershed physiography, and human land use on the water quality (summer baseflow nitrate-N and soluble reactive phosphorus concentration) in watersheds across the Upper Midwestern United States. We use an SEM to link water quality (stream nitrate-nitrogen and dissolved phosphorus) to foliar retention (AVIRIS-Classic derived foliar traits related to recalcitrance), watershed retention (wetland proportion, MODIS Tasseled Cap Wetness), runoff (agricultural and urban land use), and watershed leakiness (AVIRIS-Classic foliar nitrogen, nitrogen deposition). The SEMs confirmed that variables associated with foliar retention derived from imaging spectroscopy are negatively related to watershed leakiness (standardized path coefficient = −0.892) and positively to watershed retention (standardized path coefficient = 0.705), with features related to watershed retention and runoff exerting the strongest controls on water quality (standardized path coefficients of −0.270 and 0.331 respectively). Comparing forested and agricultural watersheds, we found significantly increased importance of foliar retention to watershed leakiness in forests compared to agriculture (standardized coefficients of −1.004 and −0.764 respectively), with measures of watershed retention more important to runoff and water quality in agricultural watersheds. The results illustrate the capacity of imaging spectroscopy to provide measures of foliar traits that influence nutrient cycling in watersheds. Ultimately, the results may help focus development and restoration policies towards building more resilient landscapes that take into consideration associations among functional traits of vegetation, physiography and climate.

The relationship between nutrient cycling and water quality in mixed-use ecosystems is driven by interactions among biotic and abiotic processes. However, the underlying processes cannot always be directly observed or modeled at broad spatial scales. Numerous empirical studies have employed land use patterns, variations in watershed physiography or disturbance regimes to characterize nutrient export from mixed-use watersheds, but simultaneously disentangling the effects of such factors remains challenging and few models directly incorporate vegetation biochemistry. Here we use structural equation models (SEMs) to assess the relative influence of foliar chemical traits (derived from imaging spectroscopy), watershed physiography, and human land use on the water quality (summer baseflow nitrate-N and soluble reactive phosphorus concentration) in watersheds across the Upper Midwestern United States. We use an SEM to link water quality (stream nitrate-nitrogen and dissolved phosphorus) to foliar retention (AVIRIS-Classic derived foliar traits related to recalcitrance), watershed retention (wetland proportion, MODIS Tasseled Cap Wetness), runoff (agricultural and urban land use), and watershed leakiness (AVIRIS-Classic foliar nitrogen, nitrogen deposition). The SEMs confirmed that variables associated with foliar retention derived from imaging spectroscopy are negatively related to watershed leakiness (standardized path coefficient = −0.892) and positively to watershed retention (standardized path coefficient = 0.705), with features related to watershed retention and runoff exerting the strongest controls on water quality (standardized path coefficients of −0.270 and 0.331 respectively). Comparing forested and agricultural watersheds, we found significantly increased importance of foliar retention to watershed leakiness in forests compared to agriculture (standardized coefficients of −1.004 and −0.764 respectively), with measures of watershed retention more important to runoff and water quality in agricultural watersheds. The results illustrate the capacity of imaging spectroscopy to provide measures of foliar traits that influence nutrient cycling in watersheds. Ultimately, the results may help focus development and restoration policies towards building more resilient landscapes that take into consideration associations among functional traits of vegetation, physiography and climate.

Introduction
Stream water quality in mixed-use watersheds is a function of the interacting biotic and abiotic factors that control ecosystemlevel nutrient cycling (Griffith 2002;Meador and Goldstein 2003;Buck et al., 2004;Ertaş et al., 2022). The forested and agricultural components of landscapes display inherently different nutrient cycling regimes, especially as consequence of anthropogenically mediated changes in ecosystem processes (e.g., nitrogen deposition, fertilizer use), soil properties (tillage), moisture regimes (irrigation) and broad-scale environmental changes related to climate (Schindler and Bayley 1993;Vitousek et al., 1997;Mosier 1998;Compton and Boone 2000). Nutrient management in mixed-use landscapes is important to the reduction of high nitrogen (N) and phosphorus (P) concentrations in streams that contribute to eutrophication and acidification of receiving waters and the loss of biological diversity and estuarine productivity, and present a public health concern (Vitousek et al., 1997;Friedland et al., 2021). Agricultural fertilization is the primary nonpoint source of nitrate-nitrogen (NO 3 -N) to receiving waters, but forest functional properties and disturbances to forests can also affect nutrient export, especially of nitrate-N K. N. Eshleman et al., 2009;Shortle et al., 2021). Like nitrate, streamwater phosphorus originates from a variety of sources, but in general is derived from nonpoint sources associated with agricultural and urban land use, and point sources associated with urbanization. Phosphorus is best predicted by urban and industrial activity in most cases and in some cases agriculture (Pieterse et al., 2003;Zampella et al., 2007;Coskun et al., 2008;Rosov et al., 2020).
The functional properties of vegetation affect nutrient dynamics at the watershed scale, even when decoupled from environmental factors such as N deposition and soil C: N ratios (Huang et al., 2011). The response of plants to nutrient availability depends on their physiology and climatic constraints on seasonal N availability (Arain et al., 2006;Seaton et al., 2019). Plant N availability and subsequent variations in local N cycling rates are also functions of competition within local species assemblages and the attendant differences in physiological trait adaptations (Reich et al., 1999;Reich et al., 2003;Huang et al., 2011). Therefore, it is notable that even within forested ecosystems, nutrient cycling rates can differ between deciduous and conifer dominated ecosystems (Aber and Driscoll 1997;Melvin et al., 2015), and may vary considerably within landscapes dominated by a single functional type . This is partly because of differences in foliar nitrogen and lignin/cellulose concentrations that characterize leaf lifespan (Wright et al., 2005;Shipley et al., 2006;Santiago 2007), influence soil C: N ratios  and, in combination with disturbance events (Eshleman et al., 2000;McNeil et al., 2007), may be strong controllers of nutrient cycling (Aber et al., 1991;Fortunel et al., 2009b;de Bello et al., 2010). Further, N output (or "leakage") from forested ecosystems has also been linked to disturbance and land use legacies (Aber and Driscoll 1997;Chen et al., 2004;McNeil et al., 2008). For example (McNeil et al., 2007), demonstrated that N leakage from watersheds may be mediated by insect-related defoliation and could maintain N limitation in temperate forest ecosystems. Similar results have been obtained in the Hubbard Brook Experimental forest (Li et al., 2004) and the Chesapeake Bay watersheds (Eshleman et al., 2000;Townsend et al., 2004;Eshleman et al., 2009;Eshleman and Sabo 2016).
In most mixed land cover watersheds, agriculture is the most significant contributor of nutrients to receiving waters (Howarth et al., 1996;Carpenter et al., 1998;Howarth 1998;Hutchins et al., 2010;Roberts and Prince 2010). The major underlying factors are fertilizer application and land management practices (Mattikalli and Richards 1996;Cicek et al., 2010;Hutchins et al., 2010;Roberts and Prince 2010) and resultant direct nutrient runoff during rainfall events (Reay et al., 1992;McCarty et al., 2008;Morari et al., 2012). For pasture-dominated landscapes, N and P export has also been linked to subsidies from animal waste (Worrall and Burt 1999;Pieterse et al., 2003;Buck et al., 2004). Agricultural nutrient management thus forms an important part of water resources protection strategies (Shepard 2005) and it has been shown that even small improvements in agricultural practices can result in significant improvements in water quality indicators (Diebel et al., 2008;Diebel et al., 2009;Zimmerman et al., 2019).
Numerous studies have shown landscape-scale variables to be good predictors of aquatic habitat quality (Lyons et al., 1996;Fitzpatrick et al., 2001;Rooney and Bayley 2011) and have provided evidence of the importance of riparian zones for sustaining diverse fish communities in streams (Fennessy and Cronk 1997;Lowrance et al., 1997;Fitzpatrick et al., 2001;Meador and Goldstein 2003). A growing body of literature has empirically linked water quality in receiving waters with landcover associations (L. B. Johnson et al., 1997;Basnyat et al., 2000;Griffith et al., 2002;Buck et al., 2004;Stanley and Maxted 2008;Singh et al., 2013), but disentangling the relative influence of vegetation traits (as opposed to cover type), climate, disturbance and watershed physiography remains to be explored due to the complex and possibly nonlinear interactions between these factors.
Multispectral and hypertemporal satellite sensors have long been employed to generate land use and land cover maps, spatially explicit estimates of forest disturbance (Kennedy et al., 2010;Townsend et al., 2012;Soulard et al., 2017), and to derive parameters describing vegetation phenology (Reed et al., 1994;Stöckli et al., 2008;Schleip et al., 2009) that we hypothesize are related to water quality. Recent research has shown that spectroscopic methods can be used to characterize key plant functional traits such as foliar nitrogen (Townsend et al., 2003;Majeke et al., 2008;Martin et al., 2008;Singh et al., 2015;Wang et al., 2019), specific leaf area (SLA), foliar lignin and cellulose (Wright et al., 2002;Majeke et al., 2008;Kokaly et al., 2009;Dybzinski et al., 2013) and potentially δ 15 N concentrations (Kleinebecker et al., 2009;Elmore and Craine 2011;Singh et al., 2015) making it possible for these functional traits to be incorporated into analyses of the drivers of water quality. However, the interactions among these potential drivers have not been comprehensively evaluated, especially given the complexity of existing process models and the absence of data needed to quantify these traits across most regions. So, while a large body of research has linked the role of physiological trait associations (Santiago et al., 2004;Shipley et al., 2006;Fortunel et al., 2009a;de Bello et al., 2010), climatic factors (Fausey et al., 1995;Krysanova et al., 1998;Hall et al., 1999;Hu and Ou 2013), land use patterns (Basnyat et al., 2000;Griffith et al., 2002;Chen et al., 2007) and disturbance (Eshleman et al., 2000;Townsend et al., 2004;Eshleman et al., 2009) to nutrient cycling rates ranging from the stand to the watershed scale, the relative role of canopy foliar biochemical and structural traits remain to be explored due to lack of concurrent spatially-explicit data. A major objective of this research was to test how measurements of foliar traits derived from NASA's Airborne Visible/Infrared Imaging Spectrometer (AVIRIS-Classic) relate to stream water nutrient export.
Here, we utilize a structural equation modeling approach to assess the relative influences of foliar biochemistry, watershed physiography and human land use patterns on water quality in watersheds across the Upper Midwestern United States ( Figure 1). We explore the influence of four broad groups of variables on stream water quality: 1) nutrient retention due to foliar biochemistry, 2) nutrient retention due to watershed physiography, 3) an index of human activity (landscape composition dominated with urban and agricultural land use), and 4) indicators of watershed-scale nutrient 'leakiness'. We hypothesize that: A) higher foliar recalcitrance has a direct positive effect on nutrient retention in watersheds and indirectly with overall water quality (i.e., lesser nutrient export), and is negatively correlated with indicators of watershed leakiness and indicators of human activity; B) indicators of watershed leakiness and high human activity have negative direct, influences on water quality (i.e., high nutrient export), as moderated indirectly by retention in watersheds and foliar recalcitrance. The structural model proposed is shown in Figure 2 with arrows pointing in the direction of hypothesized influence. Variables constituting the latent factors in Figure 2 are described in detail in the methods section.

Data Watersheds
Watersheds for this study were identified within 53 swaths that were imaged by AVIRIS-Classic between 2008 and 2011 in the upper Midwestern United States (Figure 1). Watersheds were selected to reflect the diversity in forest functional associations, landscape physiognomy and landcover composition. We selected first-through third-order watersheds from the NHD-plus database (NHDPlus 2010) that fell more than 90% within the AVIRIS swaths. NHD-plus stream layers intersected with selected watersheds were overlaid with TIGER (TIGER/Line 2011) road networks to identify~350 potential sampling locations. A total of 216 watersheds were selected in the field for sampling after screening for safe access and presence of water in the channel. Streams were sampled in 2010 and 2011, always within 1 year of AVIRIS-Classic acquisition, with 28 watersheds sampled in both years ( Figure 1). Our sampling was limited to 2 years to correspond as closely as possible with AVIRIS overflights. Watershed sizes ranged from 9.8 ha to 7,685 ha with a median of 725 ha and mean of 1,123 ha (S.D. 1,329 ha, Supplementary Appendix Figure S1). Approximately 25% of watersheds are predominantly forested (>70% forest), around 13% of watersheds comprise at least 25% agricultural and pastoral use. There is a broad mix of wetlands across the state but predominantly in the forested north, The maximum proportion of urban landcover is 4%.

Water quality sampling
A total of 216 water samples from wadeable streams were obtained in the late summer (August-September) of 2010 and 2011, 5-50 m upstream from culverts depending on local accessibility. Streamwater samples were filtered on-site using 0.45 micron glass-fiber filters (Whatman Plc Piscataway, NJ), stored on ice in 60 ml Nalgene bottles (Nalge Nunc International Corporation, Rochester NY), and frozen until the chemical analyses were performed. Samples were analyzed for nitrate-N (NO 3 -N) and soluble reactive phosphorus (SRP) using protocols described in the Manual of analytical methods (WSLH 1993). Nitrate-N concentrations averaged 0.353 mg/L (0.0004-4.97 mg/ L, median 0.042) and SRP concentrations averaged 0.0123 mg/L (Below detection limit −0.179 mg/L, median 0.0045). Distributions of both nitrate-N and SRP concentrations were strongly lognormal and were log-transformed prior to all analyses (Supplementary Appendix Figure S2).

Watershed characteristics
Broad indicators of watershed physiography (stream length, stream density) were obtained from NHD-plus stream network data (NHDPlus 2010) clipped to watershed boundaries obtained by terrain analysis of the National Elevation Dataset (Gesch et al., 2002) at field sampled locations. Average soil infiltration capacities for delineated watersheds were derived from the STATSGO database (Schwarz and Alexander 1995) by identifying major hydrologic soil groups, linking the hydrologic soil groups with infiltration capacities (NRCS 2007) and weighting by area within the respective watersheds. Contribution of groundwater to baseflow was obtained from Wolock (2003) and averaged within watersheds.

Climate and vegetation phenology
Daily precipitation and maximum and minimum temperatures were obtained from the DayMet database (Thornton et al., 2012), and averaged annually within each watershed to characterize an average climatological record for each watershed. Mean parameters of vegetation phenology for each watershed for the year sampled were obtained by fitting a double logistic model (Eq. 1) to NDVI time series obtained from the MODIS MOD09A1 product (500m, 8-day NDVI composites). Phenological parameters were obtained for each year and averaged over all pixels contained in each watershed, following Butt et al. (2011): Where: a = minimum NDVI, b = max-min NDVI, r i = maximum rate of NDVI increase at start of season, SOS = start of season date (first inflection point), r d = maximum rate of NDVI decline at end of season, EOS = end of season date (second inflection point). Note that only SOS was used as a parameter as most sampling occurred during the middle of the growing season.

Landcover and foliar traits
Data on proportional landcover for each watershed were obtained from the National Land Cover Database (Fry et al., 2011). Maps of predicted foliar traits characterizing canopy biochemical (%N, %C, Lignin, Cellulose) and structural parameters (leaf mass per area, LMA) were obtained from results of concurrent AVIRIS campaigns (Singh et al., 2015) and averaged for each watershed. AVIRIS campaigns were conducted during the midsummer peak greenness period spanning 2008 through 2011 and measured −800 GB of imagery for this landscape. Imagery was acquired from NASA's ER-2 platform at an altitude of ca. 14,000-20,000 m with pixel sizes ranging from 12 to 18 m.

Direct inputs and disturbance
In the absence of a continuous data record from Landsat (45% of all Landsat five and Landsat ETM + imagery was cloud contaminated in summer months), we employed disturbance indices (Healey et al., 2005) derived from the MOD09A1 product to characterize disturbance between the year of sampling and the previous year. Similar to other spatial products, averaged values of disturbance indices were assigned to each watershed for each Frontiers in Environmental Science frontiersin.org year. Data on atmospheric nitrogen inputs to watersheds were obtained from the North American Nitrogen Deposition program (NADP 2014) for years corresponding to the sampling year and averaged within watershed boundaries.

PLS path models
Structural equation models (SEMs) are a set of statistical methods that aim to estimate a network of causal relationships (Vinzi et al., 2010b). Structural relationships are constructed as recursive linkages between (unmeasured) latent complex concepts, each measured through observable indicators. Overall, the intent is to study the complexity of a system using a causality concept among latent constructs (latent variables, LVs) while describing each LV by measured observations called manifest variables (MVs). The partial least squares (PLS) approach to SEMs, also known as path modeling (PLS-PM), represents an intersection of path analysis (Tukey 1964;Alwin and Hauser 1975;Vinzi et al., 2010b) and confirmatory factor analysis (Thurstone 1931). Proposed as a component-based alternative to covariance based structural equation modeling (CB-SEM) estimation procedures by Wold (Wold 1966), the PLS-PM technique iteratively solves for blocks of the measurement model in the first step (the relation of LVs to MVs), and proceeds with the estimation of the structural model (the interrelationships between LVs) in the second step. These steps are iterated until the aggregated residual error is minimized (Dijkstra 2010). The PLS-PM approach attempts to explain the residual variance of the latent and manifest variables rather than modeling the sample covariance matrix. As such, the PLS-PM approach relaxes strict distributional and sample size requirements of data when compared to covariance-based SEMs (Lohmoller 1989). In contrast to CB-SEM analyses, the PLS-PM approach allows for formative indicator constructs to be identified, i.e., situations when the manifest variables are supposed to cause changes in latent variables (as opposed to reflective indicators, in which the relationship is reversed). Confidence intervals for parameter estimates are finally obtained empirically by bootstrapping techniques. Details of PLS-PM are described in Vinzi et al. (2010a).
We defined five latent constructs (Figure 2) that represent a conceptual description of watershed function. These were: 1) nutrient retention due to foliar recalcitrance (RETENF), characterized by: carbon to nitrogen ratio, lignin to nitrogen ratio, fiber to cellulose ratio and leaf mass per unit area; 2) nutrient retention in wetlands (RETENW), characterized by: MODIS tasseled cap wetness index, percentage wetland landcover, percentage water area, stream length, groundwater contribution to baseflow, and soil infiltration capacity and mean phenological start of season date; 3) watershed nutrient 'leakiness' (LEAKGE) , characterized by: foliar N concentration, stream density, MODIS tasseled cap brightness index, total atmospheric N deposition, MODIS disturbance index, index of aridity (ratio of precipitation to potential evapotranspiration), proportion of area under coniferous forest and the composite terrain index derived from a digital elevation model (Gesch et al., 2002); 4) nutrient runoff from human dominated landcover (RUNOFF), characterized by: proportion of urban built up area, proportion of area under agriculture and pasture, and finally 5) water quality (WTQUAL), characterized by measured streamwater concentration of nitrate-N (NO 3 -N) and soluble reactive phosphorus (SRP). Table 1 presents the hypothesis associated with each variable, and Table 2 lists basic statistical information on all variables, and presents the hypothesis associated with each. In brief: for foliar retention (RETENF) we assume that the combination of high C: N ratios, high Nitrogen to Lignin ratios, high fiber to cellulose ratios, and high leaf LMA will in combination reduce nutrient exports from watersheds. For retention in watersheds (RETENW), a high tasseled cap wetness index, high proportion of wetland landcover, large proportion of water bodies, longer stream lengths, higher soil infiltration capacity, and an early start of season will reduce export. Note that we do not assign a hypothesis for baseflow contribution. For the potential for fertilizer-related subsidies from human dominated landcover (RUNOFF), we hypothesize that higher proportions of agricultural, urban (via lawns and golf courses), and pasture (via export from livestock operations) will cause higher exports. For watershed 'leakiness' (LEAKGE), we assume high foliar N concentrations (with higher N in cropland), higher stream density (i.e., shorter in-stream processing), higher tasseled cap brightness index (as an indicator of exposed soil), higher disturbance (derived from MODIS), higher aridity, and higher N subsidies from atmospheric deposition will cause higher exports, while larger proportions of coniferous landcover, and high values of the composite terrain index indicating wetlands or other lowlands areas will reduce nutrient export. For the inner model (Figure 3), we hypothesize that foliar retention (RETEN) will reduce nutrient export, especially in combination with forested and wetland landcover (RETENW). We hypothesize that higher proportion of agricultural and pastoral landcover (RUNOFF) will increase nutrient export but should be mediated by foliar retention and watershed retention. We hypothesize that watershed properties influencing nutrient leakage (LEAKGE) should increase export, especially in combination with agricultural land cover, but should also be mediated by foliar and wetland retention.
We specified RETENF, RETENW, and LEAKGE as formative constructs (i.e., latent variables defined as being 'caused' by manifest variables), and RUNOFF and WTQUAL as reflective constructs (i.e., latent variables that 'cause' increases in manifest variables). We bootstrapped the model 500 times to obtain uncertainty estimates on all path coefficients. All analyses were conducted using the plspm package (Sanchez 2013) in R (R 2020).
Initial models were evaluated for specification appropriateness by inspecting measures of unidimensionality of the latent blocks. Unidimensionality metrics ensure satisfaction of the implicit assumption that manifest variables are better related to their own latent variable than others, and Kaiser's rule that the first Eigenvalue of the correlation matrix should be higher than one while others are smaller. We used Dillon-Goldstein's Rho (Chin 1998) and principal components analysis of each block to check for unidimensionalty following Vinzi et al. 2010a andSanchez (Sanchez 2013). We also report the classical Cronbach's alpha for each model construct. The model was termed suboptimal if it failed unidimensionality tests, namely if 1) the Dillon-Goldstein's Rho was less than 0.7, and 2) if the first and second Eigenvectors were both higher than 1.0 (Kaiser's rule). Model fits were further considered suboptimal if 3) the average communality index of a construct was less than 0.5 and 4) if the goodness-of-fit was less than 0.5. In brief, the communality index measures how much of the variability in a manifest variable is explained by the variability in its latent variable score (Vinzi et al., 2010b). The average communality of all MVs in a latent variable should be at least greater than 0.5. The goodness-of-fit statistic is simply the geometric mean of the average communality index and the average R 2 of each latent variable.
To optimize the SEM, we 1) inspected the model outer correlation matrix to identify the manifest variables having higher correlations with latent constructs other than the ones they were initially assigned to, and 2) evaluated bootstrap confidence intervals of weights of under-performing manifest variables. We dropped those manifest variables that had low predictive power (non-significant path weights) and loaded highly on latent variables outside of their construct (high cross-loadings that did not make ecological sense, for example, if stream length were to load highly on foliar traits.) We report both the fully specified model and the final reduced model from this iterative process.
We conducted bootstrap tests with replacement to test whether nutrient cycling mechanisms differed between watersheds that were predominantly forested (>70% forest)

FIGURE 3
Fitted (reduced variable) structural model. Path colors denote whether the antecedent latent variable positively (blue) or negatively (red) effects the descendent latent variable. For example, RETENF increases RETENW, but decreases LEAKGE. Numbers adjacent to links denote standardized coefficient sizes. See Table 6 for significance levels.

Results
Unidimensionality tests indicated suboptimal fits in the fully specified model (Supplementary Table S2). Inspection of the outer correlation matrix revealed that multiple MVs loaded higher on other variables than the proposed construct indicating some amount of model miss-specification (Supplementary Table S1) and bootstrapped estimates of path weights (Supplementary Table S2) revealed that most, if not all, such MVs did not have significant path weights (p-value < 0.05). We iteratively dropped each non-significant MV and refit the model until all unidimensionality measures were satisfied (Table 3). After iteratively dropping non-significant manifest variables, we were left with two manifest variables in each construct (Table 4). The reduced model indicated an adequate fit with the Dillon-Goldstein's Rho >0.7 for both reflective LVs (RUNOFF, WTQUAL, Table 3B). The reduced model also satisfied Kaiser's rule (second Eigenvectors of all LVs less than 1.0; Table 3) and inspection of the outer correlation matrix confirmed the unidimensionality of the construct with every MV loading highly on its own LV (Table 4). Average communalities for all LVs were also greater than 0.5 (RETENF: 0.939, RETENW: 0.539, RUNOFF: 0.922, LEAKGE: 0.866, WTQUAL: 0.656). The reduced model had a higher goodness-of-fit (0.697) compared to the fully-specified model (0.536) ( Table 3).
The reduced model contained two manifest variables per each latent variable (Table 5). Nutrient retention due to foliar recalcitrance (RETENF) was best explained by foliar C: N ratios followed by foliar Lignin:N ratios. Retention in watersheds (RETENW) was best explained by the tasseled cap wetness index followed by proportional landcover classified as lakes.
Runoff from human-dominated land use (RUNOFF) was equally well explained by proportional area under agriculture and pasture. Watershed 'leakiness' (LEAKGE) was best explained by foliar nitrogen concentration and atmospheric N deposition.
When paths were aggregated over all possible linkages (i.e., Total effects, Table 7), the model indicates that nutrient retention due to foliar retention had a significant effect on reducing watershed-scale nutrient export (RETENF → ↓WTQUAL; p < 0.0001), and reducing nutrient export from human-dominated landcover (RETENF → ↓RUNOFF; p < 0.0001). The model did not exhibit significant linkages between human-dominated land use patterns and watershed nutrient leakiness, or between watershed leakiness and water quality. Table 8 shows inner model fits that describe how well manifest variables described latent variables. Note the fit RETENF is not Comparing path coefficients between mostly forested (>70%) and other watersheds (Table 9; Figure 4), we found that nutrient retention due to foliar recalcitrance was the predominant factor mitigating nutrient leakage in forested watersheds (RETENF → ↓LEAKGE; p = 0.002). However, in other watersheds (dominated by agriculture in our study area), wetland retention played in more prominent role in ameliorating water quality directly (RETENW → ↓WTQUAL; p = 0.034) as well as through reduction of runoff (RETENW → ↓RUNOFF; p < 0.0001) and leakage (RETENW → ↓LEAKGE; p = 0.022).
Latent variable scores exhibited strong latitudinal gradients; foliar and watershed retention increased with latitude (RETENF*Latitude r = 0.83, RETENW*Latitude r = 0.74; both p < 0.0001) and runoff and watershed leakage declined (RUNOFF*Latitude r = −0.71, LEAKGE*Latitude = −0.91 respectively; both p < 0.0001), consequently translating to a large latitudinal gradient in water quality indices (WTQUAL*Latitude r = −0.71, p < 0.0001) that parallels agricultural land use patterns in the Upper Midwest (% [Agriculture + Pasture]*Latitude r = −0.72, p < 0.0001). Interestingly, this pattern is independent of the gradient in forest cover in surveyed watersheds (%Forest*Latitude r = 0.065, p = 0.387). These results show that foliar retention as derived from AVIRIS-Classic has a stronger relationship to ecosystem processes than traditionally used landscape variables such as percent land cover. Maps of latent variable scores derived from applying path weights to manifest variables show large amounts of spatial variation in projected litter recalcitrance between predominantly human dominated agriculture-woodlot landscapes (predicted low foliar recalcitrance Figure 5A), mixed-hardwood and wetland dominated forests (predicted moderate foliar recalcitrance Figure 5B) and mixed-coniferous wetland dominated landscapes (predicted high foliar recalcitrance Figure 5C) across the Midwest.

Discussion
We explored the relative influence of canopy foliar traits, watershed physiography, and human land use patterns on measures of water quality in first-through third-order streams in the Upper Midwest, United states. We employed structural equation modeling to relate water quality measures (nitrate-N and soluble reactive phosphorus in stream water) with four broad latent variables: 1) nutrient retention due to foliar biochemistry, 2) nutrient retention due to watershed physiography, 3) an index of human activity (landscape composition dominated with urban and agricultural land use), and 4) indicators of watershed-scale nutrient 'leakiness'. First, we found that higher foliar recalcitrance, as derived from AVIRIS-estimated C: N and lignin, had a direct positive effect on nutrient retention in TABLE 7 Total effects of the final (reduced variable) model, estimates denote whether an antecedent latent variable ("From" column) positively or negatively affects the descendent variable ('To' column) while accounting for all paths connected through other latent variables (i.e., direct + indirect effects). Path weights are presented as estimates from 500 bootstrap estimates, accompanied by interquartile confidence intervals (LCI, UCI), associated T statistics and P-values. Effects that are not significant (p < 0.05) are grayed out for clarity. watersheds and indirectly with overall water quality (i.e., lesser nutrient export), and was negatively correlated with indicators of watershed leakiness and indicators of human activity. Second, indicators of watershed leakiness, as derived from NADP deposition and AVIRIS-estimated %N, and high human activity (derived from NLCD land cover) had direct negative influences on water quality (i.e., high nutrient export) which could be moderated indirectly by retention in watersheds (Tasseled Cap wetness and water body cover) and foliar recalcitrance. The comparison of structural models for mostly forested watersheds with agricultural watersheds showed that nutrient retention due to foliar recalcitrance was the major factor determining watershed nutrient leakiness in forested ecosystems (Table 9; Figure 4). This finding agrees with earlier studies that have associated recalcitrant foliar traits with reducing litter decomposability and affecting eventual nitrogen cycling rates in several ecosystems (Aber et al., 1991;Verchot et al., 2001;Fortunel et al., 2009a;Kazakou et al., 2009), especially ones with high rates of N deposition and disturbance.
Similarly, we found that indicators of wetland landcover mitigated nutrient export to streams. Importantly, the comparison of agricultural and forested watersheds (Table 9; Figure 4) revealed that watershed retention formed the most important factor mitigating runoff, leakage and eventual water quality in agricultural watersheds. Natural wetlands are known to be important sinks of nutrients (Howard Williams 1985;Johnston 1991;Uusi-Kamppa et al., 2000;Saunders and Kalff 2001), and have led to the proliferation of constructed wetlands for ameliorating water quality issues (Mitsch et al., 1995;Uusi-Kamppa et al., 2000;Mitsch et al., 2005;Vymazal 2007). While the MODIS wetness index was a strong predictor of nutrient retention in watersheds, proportional area mapped as wetlands (from NLCD 2006) was not. Indicators of human dominated land use also behaved as expected with higher proportional areas under agriculture and pastures indicating higher nutrient export to watersheds. Intensive agricultural (Allan et al., 1997;Johnes and Heathwaite 1997;Howarth 1998;Boesch et al., 2001;Hutchins et al., 2010) and animal feeding operations (Carpenter et al., 1998;Pieterse et al., 2003;Buck et al., 2004) have long been identified as detrimental to the water quality of receiving waters, and have attracted considerable attention for targeted management and establishment of total maximum daily loads (Reckhow et al., 2005;Migliaccio and Srivastava 2007;Srivastava et al., 2007;USEPA 2010). Indicators of nutrient leakage from watersheds, however, were not significant predictors of water quality (Table 6, 7). The leakage latent variable consisted of measures of foliar N concentration (possibly an indicator of fertilizer application and/or differential nutrient

FIGURE 4
Comparison of inner path coefficients of models built for predominantly forested watersheds (>70% forested, green) vs. more agricultural ones (red). *p < 0.05, **p < 0.01, *p < 0.001 .   TABLE 9 Comparison of path coefficients between predominantly forested (>70% forested) and more agricultural watersheds. Estimates of Global path coefficients ('Estimate' column in Table 6) are presented along with differences, the bootstrap T-statistic of the difference (T), the degrees of freedom and the associated p-value. Comparisons that were not significant grayed out for clarity.

Path
Global uptake by vegetation) and N deposition (direct inputs), both ostensibly strong predictors of nitrogen export (Aber and Driscoll 1997;Fenn and Poth 1999;Ollinger et al., 2002;Paerl et al., 2002). These effects may have been better captured by the latent vector for human land use (RUNOFF), which explicitly coded for patterns of agricultural cover types. Post-hoc tests revealed that this might also be an effect of the strong correlation of spatial N deposition patterns in the Upper Midwest with agricultural land use intensity (r = 0.66, p < 0.0001) and proportional land use under pasture (r = 0.55, p < 0.0001). As such, the results do not suggest that foliar N and N deposition are unrelated to water quality, only that patterns in these variables are weaker predictors than measures of land use. Our findings are in general agreement with other studies conducted in the Upper Midwest in general and Wisconsin in particular (Fitzpatrick et al., 2001;Robertson et al., 2006;Rogers et al., 2009;Singh et al., 2013), but our work builds on these studies by incorporating measures related to the role of vegetation functional traits in watershed nutrient export. Our effort was made possible by recent developments in imaging spectroscopy technology (Ustin et al., 2004;Ollinger and Smith 2005;Kokaly et al., 2009;Singh et al., 2015;Wang et al., 2019) that enable landscape-scale estimation of foliar biochemical and morphological traits known to be strong drivers of ecosystem function and nutrient cycling (Reich et al., 1992;Shipley and Lechowicz 2000;Wright et al., 2004;Shipley et al., 2006). High foliar lignin to nitrogen ratios (Robinson and Jolidon 2005;Hobbie et al., 2007;Johnson et al., 2007) and leaf dry matter content (Garnier et al., 2004;Kazakou et al., 2006;Quested et al., 2007;Fortunel et al., 2009a;Fortunel et al., 2009b) slow litter decomposition rates, but linking measures of litter recalcitrance with landscape-scale processes has remained difficult. We demonstrate that information obtained from imaging spectrometry can be used to simultaneously link ecosystem characteristics such as foliar biochemistry with anthropogenic stressors such as land use patterns to assess landscape-scale responses of ecosystems to perturbations. Indeed, information obtained from imaging spectrometers is increasingly being used for assessing ecosystem attributes such as foliar biochemistry (Asner et al., 2007;Kokaly et al., 2009;Singh et al., 2015;Wang et al., 2019), nutrient cycling (Martin and Aber 1997;Ustin et al., 2004;Martin et al., 2008;Kokaly et al., 2009), invasive species (Lawrence et al., 2006;Underwood et al., 2006;He et al., 2011;Somers and Asner FIGURE 5 Maps of latent variable RETENF (nutrient retention due to foliar recalcitrance) derived by applying path weights (Estimate column, Table 5) to AVIRIS-derived maps of C: N ratio and lignin: N ratio. Brighter colors indicate projected low litter recalcitrance due to high litter quality. Panels on the right are land cover maps derived from NLCD 2006. Letters refer to location of sites in Figure 1. Polygons indicate watersheds used in the study.

Frontiers in Environmental Science
frontiersin.org 2013) and mapping canopy fuels (Ha et al., 2006;Jia et al., 2006), and will be increasingly available to scientists and managers with data from missions such as PRISMA and HISUI, and planned missions such as EnMAP (Stuffler et al., 2007), CHIME and SBG. CHIME and SBG have potentially great potential, as these missions are intended to secure global coverage with biweekly or better repeat visits, enabling characterization of dynamics linkages between landscape or vegetation processes and water quality. As well, airborne programs such as the U.S. National Ecological Observatory Network's (NEON) Airborne Observation Platform (Kampe et al., 2010) can be linked to ongoing aquatic sampling by NEON to further test relationships between watershed variables from imaging spectroscopy and aquatic response variables. Land use and ecosystems link in multiple ways that influence water quality, although the connections may be indirect (Maloney and Weller 2011). Whereas it is relatively easy to monitor human-driven ecosystem stressors such as land use patterns as a basis to infer fertilizer application rates, the assessment of the role of landscape-scale, apparently unobservable determinants (e.g., decomposition rates) can be difficult. Our study provides a methodology that explicitly links measures of human land use, watershed physiography, ecosystem-wide nutrient subsidies and foliar biochemistry. This could facilitate a spatially explicit approach to targeting management interventions across large regions. We leverage the power of structural equation models to allow the formation of latent constructs that, although 'unobservable' in the strict sense, are based on measurements of ecosystem and anthropogenic indicators that have found wide empirical support in previous research and are used in guiding watershed management efforts worldwide.

Conclusion
Overall, we found that canopy foliar chemical and structural traits affected nutrient export across a broad region. These findings are made possible by the use of imaging spectroscopy data. Once such imagery becomes widely available, they will provide the basis for much more comprehensive assessments of the drivers of water quality at landscape and broader scales. Our results provide a potential framework to guide landscape management decisions such as designing nutrient conserving landscapes (Diebel et al., 2009). For example, 'high leakage' watersheds could be amended with strategically located wetlands or riparian buffers (Mitsch et al., 1995;Diebel et al., 2009;de Souza et al., 2013) when topographic modifications (e.g., grading) are not possible. Management agencies could target watersheds with varying amounts of disturbance regimes, soil properties and/or physiography (Smith et al., 1997;Brakebill and Preston 2003;Brakebill et al., 2010) by focusing on those factors that are controllable by either natural or management interventions in specific landscapes. For example, the comparisons in Figure 4 suggest that foliar traits play a larger role in mediating nutrient runoff and mitigating effects of N deposition in forests than in largely agricultural watersheds. In contrast, watershed retention (i.e., wetlands) may play a larger role in mitigating runoff and ameliorating effects of N deposition in agricultural compared to forested watersheds. Managing and monitoring ecosystems to ensure balanced delivery of multiple ecosystem services is a key challenge for applied ecology (de Bello et al., 2010). Our study provides a step in this direction: by identifying ecological associations in addition to landscape-scale physiographic and climatological variables, we enable insights into a range of drivers of water quality that may help focus development and restoration policies towards building more resilient landscapes. At the same time, utilizing recent advances in satellite and airborne imaging technologies may make this process more standardized and available to a broader audience of managers and stakeholders.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions
AS and PT designed the study. AS collected field data, processed and analyzed imagery, conducted analyses, and wrote the initial manuscript. PT contributed to interpretation of the results, as well as writing, editing and proofreading the manuscript.

Funding
This work was conducted under NASA Terrestrial Ecology grant NNX08AN31G and Applied Sciences grant NNX09AO15G to PT, with additional support from NSF ASCEND Biology Integration Institute award DBI-2021898 and the NASA Surface Biology and Geology (SBG) mission through funding from CalTech Jet Propulsion Laboratory award 1673139.

Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher's note All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.