Relating spatial patterns of stream metabolism to distributions of juveniles salmonids at the river network scale

. Understanding the factors that drive spatial patterns in stream ecosystem processes and the distribution of aquatic biota is important to effective management of these systems and the conservation of biota at the network scale. In this study, we conducted ﬁ eld surveys throughout an extensive river network in NE Oregon that supports diminishing populations of wild salmonids. We collected data on physical habitat, nutrient concentrations, bio ﬁ lm standing stocks, stream metabolism (gross primary production [GPP] and ecosystem respiration [ER]), and ESA-listed juvenile salmonid density from approximately 50 sites across two sub-basins. Our goals were to (1) to evaluate network patterns in these metrics, and (2) determine network-scale linkages among these metrics, thus providing inference of processes driving observed patterns. Ambient nitrate-N and phosphate-P concentrations were low across both sub-basins ( < 40 l g/L). Nitrate-N decreased with watershed area in both sub-basins, but phosphate-P only decreased in one sub-basin. These spatial patterns suggest co-limitation in one sub-basin but N limitation in the other; experimental results using nutrient diffusing substrates across both sub-basins supported these predictions. Solar exposure, temperature, GPP, ER, and GPP:ER increased with watershed area, but bio ﬁ lm Chl a and ash-free dry mass (AFDM) did not. Spatial statistical network (SSN) models explained between 70% and 75% of the total variation in bio ﬁ lm Chl a , AFDM, and GPP, but only 21% of the variation in ER. Temperature and nutrient concentrations were the most supported predictors of Chl a and AFDM standing stocks, but these variables explained little of the total variation compared to spatial autocorrelation. In contrast, solar exposure and temperature were the most supported variables explaining GPP, and these variables explained far more variation than autocorrelation. Solar exposure, temperature, and nutrient concentrations explained almost none of the variation in ER. Juvenile salmonids — a key management focus in these sub-basins — were most abundant in cool stream sections where rates of GPP were low, suggesting temperature constraints on these species restrict their distribution to oligotrophic areas where energy production at the base of the food web may be limited.


INTRODUCTION
Analysis of spatial patterns can provide insights into ecological processes that may not be apparent through traditional sampling and statistical techniques (Jeltsch et al. 1999, McIntire andFajardo 2009). This is particularly evident in stream and river ecosystems where evaluation of inherent network structure helps elucidate processes that can be missed in assessments at smaller spatial scales (Fausch et al. 2002). Spatially explicit data have recently been used to infer processes regulating network-scale stream temperatures (Isaak et al. 2017b), stream water chemistry (McGuire et al. 2014), hydrologic processes (Segura et al. 2019), the distribution of biota (Filipa et al. 2017, Isaak et al. 2017a, Saunders et al. 2018, and stream metabolism (Rodriguez-Castillo et al. 2018). The processes of primary production and ecosystem respiration (collectively referred to as stream metabolism) are integral in determining community structure of aquatic biota and ecosystem function through controls on nutrient dynamics, biogeochemical cycling, and energy flow to consumers (Bernhardt et al. 2017). Gross primary production (GPP) and ecosystem respiration (ER) are processes of particular interest in stream and river environments because they integrate physical patterns in watershed structure with ecosystem processes throughout river networks, while at the same time acting as a driver of biotic organization across the landscape (Bernhardt et al. 2017). Given these linkages, determining spatial patterns of stream metabolism can improve our understanding how stream communities are assembled and can aid in the management of these ecosystems. In this study, we use spatially explicit data throughout a stream network to evaluate relationships between nutrient dynamics, biofilm standing stocks, stream metabolism, and the distribution of two native fish species of high interest to managers (Chinook salmon, Oncorhynchus tshawytscha; Steelhead, O. mykiss) in a large Columbia River Basin tributary.
Localized experimental and observational studies have found GPP to be regulated by light availability (Bilby andBisson 1992, Roberts et al. 2007), nutrient loading (Peterson et al. 1993), and temperature (Demars et al. 2011). Stream reachscale GPP is also often closely related to standing stocks of benthic periphyton (Morin et al. 1999, Dodds 2006, Bernot et al. 2010; however, this is not always the case as the strength of the relationship between GPP and periphyton standing stocks may be affected by grazing pressure, selfshading, and time since disturbance (Uehlinger 2006, Roberts et al. 2007). Regional studies exploring relationships between GPP and abiotic stream characteristics have found that GPP is most closely associated with light flux (Mulholland et al. 2001, Bernot et al. 2010, Finlay 2011. Conversely, attempts to correlate nutrient concentrations with rates of GPP in observational field studies rarely find significant positive relationships in riverine systems (see Bernhardt et al. 2017), despite nutrient availability being clearly identified as limiting factors in reach-scale experimental nutrient additions and localized bioassays (Peterson et al. 1993, Tank and Dodds 2003, Johnson et al. 2009). This disconnect may in part be attributed to the interaction between metabolism and nutrient concentrations; through their influence on nutrient uptake and retention, these ecosystem processes both affect, and are affected by, stream nutrient concentrations , Tank et al. 2017a. In river reaches with high rates of GPP and ER, high nutrient demand can lead to a reduction in the availability and concentration of limiting nutrients (Tank et al. 2017a, b). Given the longitudinal connectivity of streams, localized areas of high uptake in one area can result in reduced supply downstream, which can in-turn create periodicity in the longitudinal patterns of nutrient concentrations , Dong et al. 2017. Exploring network patterns of stream GPP and ER (collectively stream metabolism) together with patterns of stream nutrient concentrations could reconcile the disconnect between local GPP and local stream nutrient availability.
In addition to the importance of stream GPP and stream ER on patterns in stream chemical attributes throughout a watershed, these metabolic processes integrate energy flow of linked aquatic-terrestrial food webs, and spatial variability in metabolism rates has the potential to drive the distribution and productivity of stream biota. At broad spatial scales, the ratio of GPP to ER (hereafter GPP:ER) in streams has been correlated with aquatic macroinvertebrate production in streams (Marcarelli et al. 2011). More recently, Saunders et al. (2018) found that rates of GPP throughout a stream network were positively correlated with juvenile Chinook salmon (Oncorhynchus tshawytscha) and steelhead (O. mykiss) density in a NE Oregon basin. At the stream reach scale (e.g., 100-1000 m), experimental increases in GPP and ER have been linked with increased secondary production of invertebrates and fish (Bilby and Bisson 1992, Peterson et al. 1993, Cross et al. 2006. Given these linkages, understanding spatial variation in stream metabolism is necessary for understanding how ecosystems are structured and to effectively manage land-use practices at large spatial scales (Naiman et al. 2012). However, our understanding of stream metabolism across stream networks is limited, especially regarding processes driving large-scale spatial patterns.
Recent research in the analysis of stream networks has suggested that the inclusion of spatial autocorrelation functions can account for nonindependence of samples and improve understanding of how explanatory variables are related to response variables (Isaak et al. 2014). These spatial statistical network (SSN) modeling techniques have been developed to apply autocorrelation functions that are specific to the unique spatial structures of streams (e.g., connected by flow paths). SSN models have been applied to a number of stream physical metrics (Isaak et al. 2017b, Scown et al. 2017, and they are increasingly being applied to biological data (Frieden et al. 2014, Filipa et al. 2017, Isaak et al. 2017a), but their use in application to integrated ecosystem processes such as GPP or ER has been limited (but see Rodriguez-Castillo et al. 2018). In addition to assessing potential factors influencing GPP and ER, the use of SSN models allows for prediction of GPP and ER at unsampled locations, thereby increasing spatial resolution across stream networks. Quantifying GPP and ER is logistically taxing, and accurate prediction of these metrics at unsampled locations may be particularly useful for managers and researchers interested in guiding management actions.
We evaluated and analyzed spatially explicit data to explore patterns of stream temperature, solar radiation, nutrient availability, biofilm standing stocks, stream metabolism, and the distribution of juvenile salmonids (O. tshawytscha; O. mykiss) at approximately 50 sites throughout two sub-basins of the Grande Ronde River in NE Oregon, USA. Our goals were to (1) evaluate network patterns in these metrics and (2) determine network-scale linkages among these metrics, thus providing inference of processes driving observed patterns. To accomplish these goals, we first evaluated relationships between watershed area and all variables to determine broad patterns associated with increasing stream size. This approach has been utilized in a number of studies evaluating stream nutrient dynamics, GPP, and ER (Finlay 2011, Hoellein et al. 2013. We then evaluated relationships among explanatory variables (i.e., nutrients, light, and temperature) and response variables (i.e., periphyton biomass, GPP, and ER) using SSN models and further determined how the spatial patterns of explanatory variables and how the underlying network structure of the response variables (spatial autocorrelation) could account for these processes across the watershed (Rodriguez-Castillo et al. 2018). We expected light, nutrients, and temperature to all be important factors accounting for periphyton standing stocks and GPP estimates, and we expected temperature and nutrients to be important factors predicting ER (Acuna et al. 2008, Demars et al. 2011). In addition, we expected spatial autocorrelation to be an important additional predictor variable for periphyton standings stocks, GPP, and ER as this metric can encompass a wide range of potential unmeasured factors that may not be clearly associated with the abiotic predictor variables that we selected.
This high-resolution network-scale analysis provided an opportunity to evaluate predictions of how nutrients, primary producers, and stream metabolism are spatially structured within a stream network. During periods of low flow, nutrient demand is expected to be highest and nutrient supply lowest, increasing the potential for biotic controls on nutrient concentrations (Wollheim et al. 2018). We therefore expected limiting nutrient concentrations to decrease with watershed area due to biological uptake and depletion of nutrients as has been observed in other study systems during low-flow conditions . In contrast, light availability and temperature are expected to increase with watershed area as streams widen, and the River Continuum Concept (Vannote et al. 1980) evokes these changes along a continuum of increasing stream size as drivers of increasing biofilm standing stocks, GPP and GPP:ER. However, others have postulated that local conditions within a watershed (e.g., geomorphology, lithology, and climate) and anthropogenic modifications to stream and riparian ecosystems disrupt processes along this continuum (Minshall et al. 1983, 1985, Ward and Stanford 1983, Finlay 2011. In this study with a full network ❖ www.esajournals.org perspective of 50 sites, we can evaluate these frameworks. In addition, by quantifying light, nutrients, and temperature across this network, we gain insight into mechanisms driving patterns in periphyton standing stocks and metabolism at the network scale. We also explored relationships between stream productivity and key biota in this river system (i.e., ESA-listed salmonid fish) in a comparison of spatial patterns in GPP with the spatial patterns of juvenile salmon and steelhead throughout the network. In this analysis, we sought to determine whether the areas with high densities of juvenile salmonids overlapped with areas of low vs. high background productivity. In some systems, areas of high productivity coincide with high fish abundance (Saunders et al. 2018), but if that is not the case, this could help to identify opportunities for enhancing the effectiveness of stream management.

Study sites
The study was conducted in two large subbasins of the upper Grande Ronde River in NE Oregon: the upper Grande Ronde River mainstem system (hereafter UGR), and Catherine Creek (hereafter CC). The Grande Ronde River flows north from its headwaters in the Blue Mountains into the Snake River, which then flows into the Columbia River. The climate in this region is characterized by cold, moist winters and warm, dry summers with mean daily air temperatures near the city of La Grande average À0.42°C in January and 21°C in July. Average annual precipitation across the basin ranges from 36 cm in the valleys to 152 cm in the mountains, with most of the precipitation in the mountains falling as snow (McCullough et al. 2016). Annual streamflow runoff is therefore mostly reliant on winter snowpack, with peak flows typically occurring in the spring and minimum flows occurring in the late summer prior to the onset of winter precipitation (Kelly and White 2016).
A total of 54 sites were sampled in the summer of 2016 across the two sub-basins (Fig. 1). All sites were associated with the Columbia Habitat Monitoring Program (CHaMP 2016), a program characterizing tributary spawning and rearing habitat of Columbia River salmonids. CHaMP sites were selected based on a spatially balanced design (generalized random tessellation stratified [GRTS]; Stevens and Olsen 2004) throughout current, historical, and potential habitat for Chinook Salmon and Steelhead. Sites sampled in 2016 ranged from first-order tributaries to mainstem sections; therefore, sites in this study represent the variety of habitats within these sub-basins. CHaMP site reach lengths were approximately 20 times the bankfull stream width, ranging from 120 to 600 m.

Physical variables
Estimates for physical variables demonstrated to directly or indirectly influence periphyton standing stocks, GPP, and ER based on previous research (i.e., solar access and stream temperature; see Larned 2010, Bernhardt et al. 2017 were obtained for each site using the procedures outlined in CHaMP (CHaMP 2016). Solar access -the percentage of sunlight reaching a stream surface after accounting for sun angle, topographic shading, and vegetative shading-was estimated using a SunEye (Solmetric, Sebastopol, California, USA). To obtain comparable temperature values for all sites, we used 2002-2011 August mean temperature estimates from the NorWeST stream temperature model (Isaak et al. 2017b). To evaluate NorWeST model predictions, we also empirically quantified temperature at 30 of the 54 sites in summer 2016 at 1-h intervals and found a high degree of correlation (R 2 = 0.985; P < 0.001). Watershed area for each site was obtained using the STARs package (v2.0.4; Peterson and Ver Hoef 2014) in ArcGIS (v10.3.1).

Nutrient concentration
Nutrient concentrations were sampled throughout each basin during mid-August (8/6-8/9). Nutrients were collected from all 2016 CHaMP sites (n = 54) as well as additional sites (n = 27) such as tributary junctions and unsampled tributaries to increase spatial resolution of nutrient concentrations. At each site, three replicate samples were taken from flowing water using a 60-mL syringe. Water was filtered through 25mm Whatman GF/F filters into 15-mL plastic vials. Samples were stored on ice and frozen within 10 h. Samples were kept frozen until May 2017, when they were analyzed using a Dionex 1500 Ion Chromatograph (detection limit = 2 lg/L) ❖ www.esajournals.org for nitrate-N and phosphate-P. Check standards were run along with samples to ensure machine accuracy. Nutrient samples were also collected from a subset of sites in mid-June to determine whether spatial patterns in nutrients differed as flows decreased through summer. Patterns were generally consistent between sampling events, although concentrations were slightly greater in mid-June. Given this consistency and the greater number of sites sampled in August, we focus analyses on August nutrient concentrations but provide June concentrations in Appendix S1: Fig. S1.

Nutrient limitation
Nutrient limitation was assessed using nutrient diffusing substrate (NDS) bioassays. Five sites were selected across each sub-basin (n = 10 total) to capture a range of temperature and landscape positions (Fig. 1). However, the NDS arrays at one site in UGR were tampered with, resulting in a total of four sites in UGR. At each site, a metal L-bar containing 12 poly-con cups comprised of different nutrient treatments was placed in a rifle at the downstream end of each site (Tank et al. 2017b). Cups were filled with 2% agar and one of four potential treatments; control (unamended), nitrogen addition (1 mol/L N; NH 4 Cl), phosphorus addition (1 mol/L P; KH 2 PO4), and nitrogen and phosphorus combined (1 mol/L N, 1 mol/L P). Glass fritted disks were placed on top of agar and lids with 25 mm diameter holes were closed so that the disk was firmly attached. Nutrient diffusing substrate bioassays were deployed for 21 d. Upon retrieval, chlorophyll a (Chl a) on each glass fritted disk was measured using an in situ fluorometer (bbe-moldaenke BenthoTorch, Schwentinental, Germany). Prior to measurement, cups were positioned in a plastic container filled with stream water and placed in the shade for a 30 min minimum acclimation period to account for light effects on fluorometric Chl a estimation (Kaylor et al. 2018). For each 12-cup array, the three cups of each treatment were averaged to obtain mean treatment values.

Periphyton standing stocks
Standing stocks of periphyton Chl a and ashfree dry mass (AFDM) were quantified at 50 sites following methods outlined in Kaylor et al. (2018). At each site, 11 evenly spaced transects were established and a single rock was collected from each transect, except for transect 11, where ❖ www.esajournals.org two rocks were gathered for a total of 12 rocks. Rocks were collected from the 25th, 50th, or 75th percentile of stream wetted width, and this location was altered systematically at each successive transect. The surveyor walked to the approximate stream location and then without looking down, grabbed the first rock they touched that was between 10 and 25 cm in width. A 7 cm diameter PVC pipe section was used to define a circular area to scape periphyton from the top of each rock. Periphyton from four rocks was rinsed into a single container and Chl a, and AFDM sub-samples were drawn from the pooled sample using a modified 60-mL syringe. This resulted in three replicates of Chl a and AFDM per site. The subsample was filtered in the field through 47-mm Whatman GF/F filters using a handheld vacuum pump. Filters were folded in aluminum foil and either flash frozen using dry ice or put on ice and frozen within 6 h. Chl a was quantified using acetone extraction and fluorometric analysis (Arar and Collins 1997). Ash-free dry mass samples were dried at 60°C for 24 h and then weighed to the nearest mg. Samples were then combusted at 500°C for 2 h and reweighed. The difference between dried mass and ashed mass was divided by the proportional area sampled to obtain AFDM (g/m 2 ).

Gross primary production and ecosystem respiration
Stream metabolism was estimated using single-station, open-channel methods (Grace and Imberger 2006) which utilize diurnal changes in stream dissolved oxygen (DO) to estimate rates at which oxygen is contributed to (i.e., GPP), and consumed from (i.e., ER) streams. MiniDOT optical dissolved oxygen and temperature sensors (Precision Measurement Engineering, Vista, California, USA) were deployed at summer baseflow conditions in a total of 49 sites-32 in UGR and 17 in CC-during 22 July-26 August 2016. There was limited cloud cover during the summer and no substantial precipitation events occurred during the sampling period. All sensors were placed at the downstream end of each CHaMP site in stream sections with flowing but non-turbulent water (Siders et al. 2017).
Because of limitations on the number of available loggers (n = 11), sensors were rotated throughout the 49 sites by deploying them for a minimum of three cloudless 24-h periods, and then moving them to new locations. This approach is consistent with Rodriguez-Castillo et al. (2018), in which GPP was similarly quantified at the network scale by collecting DO data at each site for three days. Consequently, metabolism was estimated at sites over slightly different time intervals. To describe temporal trends over the sampling interval, three longer-term stations were established to monitor metabolism throughout the duration of sampling other sites (n = 2 in UGR and n = 1 in CC; see Appendix S1: Fig. S2). At each of these stations, photosynthetically active radiation (hereafter PAR; Odyssey sensors; Dataflow Systems, Christchurch, New Zealand) and barometric pressure (Extech RHT50 sensors; FLIR Commercial Systems, Nashua, New Hampshire, USA) were measured at 5-min intervals.
Metabolism was estimated using a Bayesian single-station estimation program (BASE v3.0; Grace et al. 2015) in the program R (R Core Team 2015). The BASE program simultaneously estimates GPP, ER, and the reaeration coefficient (K) through an iterative process using a Markov Chain Monte Carlo (MCMC) method. Additional parameters for light saturation (p) and temperature dependence (h) can optionally be estimated by the model. Five-parameter models, in which estimates of GPP, ER, K, p, and h were derived, always outperformed (i.e., lower AIC values) 3-parameter models where values for p and h were fixed. Therefore, 5-parameter models were used for all metabolism estimates (Grace et al. 2015). We used default model priors for GPP, ER, p, and h. The prior for K for each site was derived from nighttime regression utilizing diurnal DO concentrations and DO saturation (Hall and Hotchkiss 2017). We applied a standard deviation of the K prior distributions that allowed for deviation from the nighttime regression estimate but reduced extreme daily K estimate outliers and improved consistency among days at each site. The default model priors for BASE restrict daily K values to be less than 40 d À1 . However, some sites in this study were high gradient with K values that may exceed this value. We therefore set the maximum K to 60 d À1 , so we did not force K to be lower than expected based on nighttime regression.
Model input requirements include light flux (PAR lmolÁm À2 Ás À1 ), dissolved oxygen concentration (mg/L), temperature (°C), barometric pressure (ATM), and salinity (ppt). Dissolved oxygen and temperature were obtained empirically at each site. Due to limited sensors, PAR could not be recorded at all stations; PAR was continuously recorded at permanent stations and extra sensors (n = 3) were deployed at additional sites. When PAR was not empirically recorded at a site, data from the nearest site (<5 km) were used. Barometric pressure was recorded at each of the three permanent stations, and data were modeled for other sites based on differences in elevation. Salinity was assumed to be zero for all sites.
All models were run with 20,000 MCMC iterations (10,000 burn-in iterations). For each day, model performance was assessed using criteria outlined in Grace et al. (2015). Only models with all R-hats <1.1, posterior predictive checks between 0.1 and 0.9, and r 2 of predicted DO values vs. measured DO values >0.7 were included in subsequent analysis. We excluded sites where model-estimated K exceeded 45 d À1 (n = 4) as high K values make it difficult to obtain reliable GPP and ER estimates Hotchkiss 2017, Appling et al. 2018). Model estimates of GPP and ER for multiple days at each site were averaged to obtain mean values per site. These values were multiplied by mean stream depth to obtain aerial rates of GPP and ER (g O 2 Ám À2 Áday À1 ).

Juvenile salmonid surveys
Abundance of salmonids was estimated at CHaMP sites during the period of summer low flow as described in . Depending on stream size, one or two snorkelers moved in an upstream direction the length of each reach while enumerating fish and communicating with each other their observations to avoid double counting. Snorkel counts at each site were expanded to abundance estimates using a correction factor (Jonasson et al. 2015) developed from paired mark-recapture and snorkel survey data to account for fish that were not observed by snorkelers; the correction factor was specific to habitat type (pools, riffles, runs). Aerial density of salmonids (fish/m 2 ) was calculated by dividing the corrected abundance estimates by the total surface area within the reach as determined from reach metrics collected during from CHaMP surveys.

Statistical analysis
To examine how explanatory and response variables relate to stream size, we first plotted watershed area (km 2 ) against each variable. This follows an established approach to infer ecological processes (e.g., nutrient dynamics) from spatial patterns in response variables plotted against watershed area, and to identify thresholds marking drastic shifts in measured variables .
We assessed nutrient limitation on nutrient diffusing substrates at each site using one-way analysis of variance (ANOVA) to compare Chl a accrual across the four treatments (C, N, P, N+P). We then applied Tukey's post hoc test to determine significant differences among treatments. These differences were used to assess the type of nutrient limitation at each site, and more broadly to determine the dominant nutrient(s) limiting periphyton accrual in each sub-basin (i.e., N limitation, P limitation, primarily N limited and secondarily P limited, primarily P limited and secondarily N limited, or co-limited by both N and P).
We used SSN models (Isaak et al. 2014) to evaluate potential factors influencing periphyton standing stocks, GPP and ER at the network scale. Prior to SSN model formulation, spatial data were formatted and processed using the STARs package (v2.0.4; Peterson and Ver Hoef 2014) in ArcGIS (v10.3.1). A preconditioned stream network layer (e.g., continuous stream network with all stream segments flow-oriented toward a single drainage point) was downloaded from the National Stream Internet (https://www. fs.fed.us/rm/boise/AWAE/projects/NationalStream Internet/NSI_network.html). Spatial processing to produce an SSN object was conducted following procedures outlined in Ver Hoef et al. (2014).
SSN models are based on the multiple linear regression model framework, with fixed-effect predictors (e.g., habitat or landscape covariates) and spatial autocovariance functions as random effects. Variation is partitioned into fixed effects, spatial covariance, and residual variation known as the nugget. There are three potential covariance structures for SSN models. Tail-up covariance (TU) accounts for covariance with points upstream of the designated location within the stream network, tail-down covariance (TD) accounts for spatial autocorrelation with points ❖ www.esajournals.org downstream of the designated location within the stream network, and Euclidean covariance represents autocorrelation based on linear distances not associated with a stream network. Although including multiple covariance structures (e.g., tail-up and tail-down in same model) can improve model predictions (Garreta et al. 2010), each additional covariance structure is associated with additional costs in terms of parameter estimation (n ≥ 2 parameters per covariance structure). Given our low-sample size (~50), we restricted models to tail-up covariance structures, as this autocorrelation structure was expected to best represent the processes evaluated (e.g., GPP should be influenced by upstream processes). Covariance structures can be modeled with exponential, spherical, Mariah, or linear-with-sill forms (Peterson and Ver Hoef 2010). Preliminary comparisons of models with different covariance forms (e.g., exponential, spherical) revealed little differences in root mean square prediction error (RMSPE) among models relative to differences observed between spatial and non-spatial models. We therefore used the spherical form for all covariance structures.
We used four fixed-effect variables-solar access, nitrate-N concentration, phosphate-P concentration, and stream temperature-to predict biofilm biomass, GPP and ER, as mechanisms for their control over response variables are wellestablished. All correlations of the four explanatory variables were evaluated for collinearity and exhibited Pearson's correlation coefficients <0.6. For each response variable (Chl a, AFDM, GPP, and ER), we formulated a set of candidate models based on all combinations of the four explanatory variables (15 total model structures for each response variable). To achieve assumptions of normality, nitrate-N, phosphate-P, and all response variables were natural-log transformed. Visual inspection of plotted data did not provide evidence for non-linearity among covariates or inflated variance with increasing mean of dependent variables and the analysis proceeded incorporating only linear relationships. Similarly, there was no evidence of either increasing or decreasing variance as values of the x-variable increased. Next, all candidate models were fit and model assumptions of normality and constant variance were checked using model residuals.
Candidate models were ranked with Akaike's information criterion (AIC) adjusted for small sample size (AIC c ; Hurvich and Tsai 1989). To evaluate the potential importance of fixed-effect variables in explaining variation in response variables, we calculated the relative importance of each variable based on the sum of Akaike weights for all models containing each variable (Burnham and Anderson 2004). We present the three top-ranked models (lowest AIC c values).
In summer 2016, we obtained reliable data from 45 to 50 CHaMP sites for each response variable; however, the usage of SSN models allows for prediction at unsampled locations as long as the fixed-effect explanatory variables in each model are also available at each prediction location. We therefore predicted each response variable at 52 unsampled CHaMP sites, where explanatory covariates were available. For each response variable, the model with the lowest AIC c value from the set of candidate models was used to predict values at unsampled locations.

Physical attributes
The UGR drains a larger area than CC, and therefore, sites within UGR encompassed a greater range of watershed areas (22-285 km 2 vs. 38-1405 km 2 ). Despite differences in watershed area, discharge at the farthest downstream site in CC was 0.85 m 3 /s, whereas discharge at the farthest downstream site in UGR was 1.34 m 3 /s. Solar exposure and temperature increased with watershed area in both sub-basins (Fig. 2a, b), but solar exposure and temperature were typically lower for a given watershed area in CC compared to UGR. Solar exposure ranged from 34% to 83% (mean = 57%; SD = 12.8) and from 31% to 98% (mean = 68%; SD = 16.8) in CC and UGR, respectively. Mean August temperature ranged from 10.5°to 16.6°C (mean = 13.5°C; SD = 1.9) and from 11.8°to 19.8°C (mean = 15.8°C; SD = 2.4) in CC and UGR, respectively.

Nutrient spatial patterns
Spatial patterns in nutrient concentrations differed between UGR and CC (Fig. 2c). In UGR, phosphate-P concentrations were generally higher (>10 lg/L) in sites with watershed areas of less than 100 km 2 compared to sites with ❖ www.esajournals.org watersheds greater than 100 km 2 where concentrations were less than 10 lg/L. Nitrate-N concentrations were low throughout UGR but were elevated in some sites with watershed area less than 100 km 2 . Although the nitrate-N:phosphate-P ratios were always <16, the very low nutrient concentrations and decreasing trend of both nutrients would suggest co-limitation of biofilm chlorophyll a (lg/cm 2 ), (f) biofilm ash-free dry mass (g/m 2 ), (g) gross primary production (g O 2 Ám À2 Áday À1 ), (h) ecosystem respiration (g O 2 Ám À2 Áday À1 ), (i) the ratio of gross primary production to ecosystem respiration, and (j) salmonid density (# m À2 ). Red points indicate sites from Catherine Creek, and blue points indicate sites from the upper Grande Ronde River. N and P in UGR, especially in sites with watershed areas greater than 100 km 2 . In CC, nitrate-N decreased with watershed area, whereas phosphate-P remained elevated across sites independent of watershed area, suggesting greater relative demand for N and therefore potential N limitation. Elevated phosphate-P concentrations in CC may be indicative of younger, more erodible basalts compared to older, more weathered underlying geology in UGR.

Nutrient diffusing substrates
Nutrient diffusing substrates were used to test predictions of nutrient limitation inferred from spatial patterns of nitrate-N and phosphate-P throughout the stream network. In UGR, we observed shifts from N limitation in the upper reaches to co-limitation and then primary N limitation with secondary P limitation in sites with the largest watershed area (Fig. 3a). At site 1 in UGR, which had the smallest drainage area (63 km 2 ), Chl a accrual on N-amended and N+Pamended substrates were significantly different (Tukey's post hoc comparisons, P < 0.05) from C-amended and P-amended substrates but were not different from each other (P > 0.05), indicating N limitation. At site 2 (356 km 2 ), Chl a on N-amended substrates was greater than C-amended and P-amended substrates, but this was not significant at a = 0.05. N+P-amended substrates were significantly different than C-, N-, and P-amended substrates indicating co-limitation. At sites 3 (512 km 2 ) and 4 (1005 km 2 ), Chl a on N-amended substrates was significantly greater than C-amended substrates (P < 0.05) and N+P-amended substrates were significantly greater than N-amended substrates indicating primary N limitation and secondary co-limitation.
In Catherine Creek, NDS responses to treatments consistently demonstrated N limitation (Fig. 3b). At every site (ranging in watershed area from 24 to 279 km 2 ), Chl a accrual on N-amended substrates was significantly greater than control substrates (P < 0.05), and Chl a accrual on N+Pamended substrates were significantly greater than controls at 4 of the 5 sites, but were not significantly different from N-amended substrates at any of the sites. P-amended substrates were not significantly (P > 0.05) different from control substrates at any site.

Periphyton standing stocks
Chl a and AFDM were positively correlated (r 2 = 0.73) and exhibited similar spatial patterns (Figs. 2e, f and 4a, b). Chl a ranged from 0.4 to 13.1 lg/cm 2 and AFDM ranged from 5.1 to 48.1 g/m 2 ; however, Chl a and AFDM were considerably lower in CC with maximum values of 2.5 lg/cm 2 and 17.6 g/m 2 , respectively. Across both sub-basins, there were no consistent trends between watershed area and Chl a or AFDM (Fig. 2e, f). However, we were able to predict 70-73% of the variation in Chl a and AFDM across this stream network using spatial autocorrelation, temperature, and nutrient concentrations (Table 1). Of this variation, the proportion attributed to fixed effects (e.g., temperature and nutrients) was only 0.16-0.20 and 0.24-0.30, while the proportion attributed to autocovariance was 0.60-0.63 and 0.44-0.45 for Chl a and AFDM, respectively. Relative variable importance averaged across all models indicated that temperature (0.55), nitrate-N (0.26), and phosphate-P (0.17) were the most important fixed-effect variables explaining Chl a (Table 2). Temperature was also the highest ranked fixed-effect explaining AFDM (0.59), followed by phosphate-P (0.32) and nitrate-N (0.09; Table 2).

Gross primary production and ecosystem respiration
Reliable model performance statistics (Grace et al. 2015) were achieved for at least two full days in 45 of the 49 sites where dissolved oxygen sensors were deployed. Gross primary production estimates ranged from 0.03 to 6.56 g O 2 Ám À2 Áday À1 and ER estimates ranged from 0.36 to 6.87 g O 2 Ám À2 Áday À1 (Figs. 2g, h and 4c, d); however, GPP was generally greater in UGR compared to CC (Fig. 2g). Gross primary production increased with watershed area (Fig. 2g) in both UGR and CC, but for a given watershed area, GPP was greater on average in UGR. Ecosystem respiration increased with watershed area in UGR but in CC, ER peaked at sites with watershed area near 100 km 2 with lower estimated ER in sites with smaller and larger watershed areas (Fig. 2h). The ratio of GPP to ER (GPP:ER) increased with watershed area in both basins, but for a given watershed area was Fig. 4. Spatial patterns in (a) biofilm chlorophyll a, (b) biofilm ash-free dry mass, (c) gross primary production rates, (d) ecosystem respiration rates, (e) juvenile salmonid densities, and (f) the ratio of gross primary production to ecosystem respiration. greater on average in UGR (Fig. 2j). Surprisingly, Chl a explained little variation in GPP in both UGR (r 2 = 0.01) and CC (r 2 = 0.11).
The top three ranked models predicting GPP explained between 71% and 73% of variation across this stream network (Table 1). The proportion of this variation attributed to fixed effects was much greater for GPP than for either of the standing-stock estimates (Chl a and AFMD).
Fixed effects accounted for 0.7-0.73 of the variation, while the proportion of variation explained by autocorrelation was minimal (<0.02). Relative variable importance averaged across all models indicated that temperature (0.37) and solar access (0.36) were the most important fixed-effect predictors of GPP, followed by nitrate-N (0.16) and phosphate-P (0.11; Table 2).
In contrast to models predicting GPP, models were poor at predicting ER. The top three ranked models predicting ER only explained between 0.17 and 0.21 of the variance, and the proportion of this variance attributed to fixed effects was less than 0.10 (Table 1), indicating little predictive power of ER using the four fixed effects used in this study. Phosphate-P was the most important fixed effect based on relative variable importance averaged across all models (0.38) followed by temperature (0.31), nitrate-N (0.23), and solar access (0.09).
Metabolism data (GPP and ER) were collected continuously at three stations to identify any Notes: RMSPE is the root mean squared prediction error. LOOCV is the squared residuals between observed and predicted values through leave-one-out cross-validation. Correlation composition is the proportion the explained variance partitioned into fixed effects, spatial autocorrelation, and the residual nugget. Bold values indicate p-values < 0.05. temporal trends occurring throughout the 6week duration of the study. At these three stations (1 in CC and 2 in UGR), we observed some temporal variation, but variation within sites was generally small and changes over time were less than variation among sites (Appendix S1: Fig. S2). Consequently, we concluded that any variation in measurements due to temporal changes over our 6-week sampling interval was likely overshadowed by spatial differences.

Juvenile salmonid spatial patterns
Juvenile salmonid density ranged from 0.001 individuals/m 2 to 1.492 individuals/m 2 across the 43 sites in the two basins (Fig. 4e). In both UGR and CC, juvenile densities were greatest in headwater sections where GPP was very low, and there was a negative relationship between GPP and salmonid density across UGR (r 2 = 0.35; P = 0.001) and CC (r 2 = 0.10; P = 0.28). Similar to other studies in nearby basins , salmonid density was negatively correlated with temperature in UGR (r 2 = 0.74, P < 0.001). However, in CC this relationship was not evident (r 2 = 0.01, P = 0.680).

DISCUSSION
High-resolution sampling throughout the river network, combined with multiple analytical approaches allowed us to quantify spatial patterns of nutrient concentrations, periphyton standing stocks, GPP, ER, and the distribution of biota that are of management concern. Models were effective in explaining 70-75% of the variation in periphyton standing stocks and GPP using a mixed-effects modeling approach that accounted for spatial autocorrelation (e.g., SSN), which allowed for accurate prediction of these variables at unsampled locations where explanatory covariates were also collected (Appendix S1: Figs. S3 and S4). Temperature, light availability, GPP and GPP:ER increased with watershed area, and while there were outliers, these trends were generally consistent with predictions outlined in the River Continuum Concept (Vannote et al. 1980). However, differences in rates of GPP between sub-basins for a given watershed area suggest that local factors within these watersheds control the nature of these relationships (Minshall et al. 1983(Minshall et al. , 1985. Salmonids were most abundant in cool stream sections where rates of GPP were low, a result that contrasts observations from a nearby basin where salmonid density was positively correlated with GPP (Saunders et al. 2018). The opposing relationships between these studies potentially reflect different landscape filters that regulate the distribution of salmonids between these basins (Poff 1997). In the sub-basins evaluated in our study, temperature has been shown to influence spawning locations and limit the distribution of juvenile salmonids , which may restrict juvenile salmonids to cool, oligotrophic areas where low rates of primary production may be limiting energy flow to the food web.

Nutrients
Nutrient concentrations were not strong predictors of biofilm standing stocks or GPP in our study sites. Nutrient concentrations within a river are typically poorly correlated with GPP (see Bernhardt et al. 2017). However, experimental nutrient additions to stream ecosystems have resulted in enhanced GPP, ER, and secondary production (Peterson et al. 1993, Slavik et al. 2004, Cross et al. 2006, providing evidence that nutrient supply, rather than concentration, often limits GPP, ER, and bottom-up drivers of secondary production. Results from our nutrient diffusing substrate experiment demonstrate that periphyton Chl a accrual in UGR was primarily limited by nitrogen and secondarily co-limited by nitrogen and phosphorous, while CC was primarily limited by nitrogen alone. Detecting a nutrient response in the NDS experiment but not the statistical modeling approach could be the result of different periphyton communities colonizing artificial substrates vs. natural substrates or differences in grazing rates (Cattaneo and Amireault 1992). Alternatively, we suggest that low nutrient concentrations across both sub-basins and a small range in nutrient concentrations throughout the watershed created conditions in which nutrients would fail to emerge as strong correlates in a linear mixed-model analysis.
While there was limited evidence for a relationship between nutrient concentrations and GPP at individual points, spatial patterns suggested interactions between nutrients concentrations, productivity, and associated nutrient demand. During periods of low flow, nutrient demand is expected to be highest and nutrient supply lowest, increasing the potential for biotic controls on nutrient concentrations (Wollheim et al. 2018). Rates of production during these low-flow periods may increase demand and further drive spatial patterns . For example, nitrogen uptake (Tank et al. 2017a) and nutrient retention efficiency (Sabater et al. 2000), which affect longitudinal patterns of nutrient concentrations, were greater in streams with open canopies where autotrophic demand would be higher, compared to more shaded streams. Light availability and GPP increased with watershed area in both UGR and CC, which may indicate greater cumulative demand with downstream distance that resulted in depleted nutrients. Both nitrate-N and phosphate-P concentrations were elevated in shaded headwaters of UGR but decreased with watershed area. Conversely, nitrate-N concentrations in CC decreased with watershed area, while phosphate-P concentrations remained relatively uniform. This suggests that demand for both nitrate and phosphate is high throughout UGR, but in CC demand is high for nitrate, but not phosphate. These spatial nutrient patterns along with the NDS results suggest that during summer low-flow periods, instream demand for nutrients in UGR and CC exceeds supply and therefore exerts control on nutrient concentrations (Wollheim et al. 2018).

Patterns and predictors of GPP and ER
Due to historic difficulty quantifying stream GPP, researchers have measured periphyton Chl a and AFDM as proxies for GPP. However, periphyton Chl a and AFDM were poorly correlated with reach-scale GPP in our study. Although this contrasts with studies that report positive associations between Chl a and GPP (Morin et al. 1999, Dodds 2006, Rodriguez-Castillo et al. 2018, it is consistent with other studies (Velasco et al. 2003, Izagirre et al. 2008. We quantified periphyton Chl a on benthic substrates, but this assessment did not include aquatic macrophytes. Macrophytes were generally uncommon throughout the basin but were observed in some headwaters (as bryophytes) and warm, mainstem sites (as vascular plants as well as filamentous algae), which could lead to discrepancies if macrophytes contribute substantially to whole-reach GPP (Kaenel et al. 2000). Additionally, standing stocks reflect accrual of periphyton after grazing by organisms, and spatial differences in grazing rates may result in a decoupling of standing stocks and GPP. Regardless of the drivers of this decoupling, our results suggest that using standing stocks as a proxy for GPP could lead to inaccurate conclusions about local-and network-scale primary production across our study region.
Solar access (potential light exposure) and temperature emerged as the best predictors of GPP in this system. Metrics of solar radiation potential have been positively correlated with GPP in a number of studies encompassing a wide geographic area (Mulholland et al. 2001, Bernot et al. 2010, Finlay 2011, Hoellein et al. 2013, Tank et al. 2017a, Rodriguez-Castillo et al. 2018, Saunders et al. 2018. As streams widen, solar radiation reaching stream channels is expected to increase due to the decreased ability of riparian vegetation to shade streams (Vannote et al. 1980). However, for a given stream size, we found that solar radiation was variable, which may be attributed to historic and ongoing landuse within these sub-basins ). This heterogeneity resulted in a decoupling of the relationship between solar radiation and temperature and allowed us to include both these covariates as predictors of GPP, improving model predictions of GPP at the network scale.
Gross primary production and ER were well correlated in UGR (GPP explained 72% of the variation in ER), but not in CC (GPP only explained 5% of the variation in ER). Ecosystem respiration represents the combined carbon consumption by autotrophs and heterotrophs, and ER and GPP may be coupled through autotrophic respiration (Hall and Beaulieu 2013). GPP:ER was greater in UGR (mean = 0.67) compared to CC (mean = 0.17), increasing the potential for autotrophic respiration to influence ER in UGR. Alternatively, GPP an ER in UGR may be linked but not truly coupled if carbon fixed by autotrophs is released as dissolved organic carbon, which then acts as a key source of organic matter fueling adjacent heterotrophs (Hotchkiss and Hall 2015). Lastly, GPP and ER may be largely functioning independently but are controlled by similar factors in UGR (Hall and Beaulieu 2013). For example, increasing temperature with watershed area in UGR may have increased rates of both GPP and ER independently. Ultimately, the decoupling of GPP and ER in CC may explain why factors associated with GPP across both basins did not emerge as strong predictors of ER across both basins.
Fixed-effect covariates-temperature, solar access, and nutrient concentrations-explained very little variation in ER across these two subbasins. Further, the combination of fixed effects and autocovariance was only able to explain approximately 21% of the variation in ER. This contrasts with Rodriguez-Castillo et al. (2018) in which SSN models that included GPP and temperature as fixed effects explained up to 67% of the variation in ER across a stream network. In the UGR and CC, other factors such as the availability of dissolved, fine, and particulate organic carbon may have been more associated with ER rates in CC (Roberts et al. 2007). Identification of additional covariates that drive spatial structure and stream metabolism at the network scale will allow ecologists to gain a deeper mechanistic understanding and predictive ability.

Evaluation of RCC predictions
The high-resolution analysis in the present study allows us to evaluate how key abiotic factors expected to change throughout a river network may affect ecosystem processes. The River Continuum Concept (RCC, Vannote et al. 1980) proposes increasing light availability as streams widen to drive increased GPP and GPP:ER until streams become large enough where light attenuation and turbidity reduce benthic primary production. Our results generally support RCC predictions for low-to midorder streams, as solar radiation, GPP, and GPP:ER increased with drainage area across both sub-basins. These spatial patterns are consistent with a growing body of research demonstrating increasing GPP and GPP:ER as a function of stream size within midorder streams (Meyer and Edwards 1990, McTammany et al. 2003, Ortiz-Zayas et al. 2005, Finlay 2011, Hoellein et al. 2013, Rodriguez-Castillo et al. 2018). However, our results also support predictions that local conditions within a watershed, whether natural (e.g., meadows) or anthropogenic (e.g., land-use practices), may dictate the relationship between watershed area and metabolism within a network leading to outliers in GPP:ER that deviate from RCC predictions (Minshall et al. 1983(Minshall et al. , 1985. For example, GPP and GPP:ER were greater for a given watershed area in UGR, where temperature and solar radiation were greater, compared to CC, which was cooler and more shaded on average. Finlay et al. (2011) found a non-linear relationship between watershed area and light availability with an abrupt transition occurring near 100 km 2 in which light and GPP rapidly increased. We did not detect abrupt changes in light as a function of watershed area in UGR or CC, potentially because riparian communities and land-use legacies in these sub-basins have led to reduced riparian cover ). These effects would be greatest in the smallest streams where canopies can close over streams resulting in little light penetration. It is unclear whether restoring vegetation and shading within these sub-basins would manifest in abrupt transitions in light and GPP, as observed in Finlay et al. (2011). Ultimately, when considering many points throughout a network, we simultaneously find support for the RCC in the general trends of light, GPP, and GPP:ER, but also support for the idea that there may be deviations from the RCC model at multiple points along a river network and between watersheds.

Spatial patterns of GPP and salmonids
In UGR and CC, juvenile Chinook and steelhead densities (individuals/m 2 ) in 2016 were greatest in areas that corresponded to low rates of GPP, resulting in a negative relationship between GPP and salmonid density. This contrasts with Saunders et al. (2018) in which GPP was positively correlated with juvenile salmonid density (individuals/m) in the nearby John Day River basin in NE Oregon. However, Saunders et al. (2018) explicitly sampled within a geographic extent exhibiting temperatures suitable to support salmonids. Sampling within a broader geographic extent within the same basin as Saunders et al. (2018), previous studies evaluating fewer sites (n = 7), but a wider range of temperatures, reported negative relationships between salmonid density and stream temperature as well as proxies for primary production (periphyton biomass) and prey availability (invertebrate density; Li et al. 1994. Thus, across broad geographic extents, temperature, including thermal refugia (Ebersole et al. 2003), may exert greater control on salmonid spawning and rearing distributions (Justice et al. 2017), but within habitat exhibiting suitable temperatures, spatial patterns in GPP may influence the energy available for higher trophic levels and the distribution of juvenile salmonids (Saunders et al. 2018).
Fish and overall river management in the basins evaluated in this study are primarily focused on juvenile salmonids, but constraining our assessment of management needs to juvenile salmonid density, rather than total fish density or biomass, limits our ability to consider relationships between GPP and higher trophic levels. In cooler, headwater sections of UGR and CC, fish species composition is dominated by salmonids but warm water species increase with distance downstream (Jonasson et al. 2015). The inclusion of total fish biomass as a response variable may have resulted in a positive relationship with GPP, as has been observed in other studies evaluating metrics of productivity or invertebrate prey biomass and top predator biomass (Hawkins et al. 1983, Kaylor andWarren 2017). Nonetheless, river sections with high salmonid densities but low primary production may provide managers with targeted areas where actions can be taken to enhance primary production or food availability (Naiman et al. 2012).

SSN usage and implications
The inclusion of spatial autocorrelation can improve prediction of response variables at unsampled locations across a stream network (Isaak et al. 2014). However, a trade-off is that covariates used to predict response variables need to be available at both sampled and unsampled locations, often limiting covariates to broad geographic descriptors (e.g., elevation, drainage area). While this approach can increase predictive power throughout a network, mechanistic linkage may be lost. We were able to maintain mechanistic linkage by utilizing an extensive monitoring program-the Columbia Habitat Monitoring Program (CHaMP)-to sample a subset of sites and then predict response metrics at unsampled sites where mechanistic covariates were explicitly measured (e.g., solar access) rather than proxies to represent variables of interest (e.g., stream order as a proxy for light reaching streams). As a result of including these sites, we were able to use the most commonly associated explanatory variables of GPP (i.e., nutrients, light, and temperature) to predict response variables at unsampled locations (Appendix S1: Figs. S3 and S4), and these covariates explained far more variation in GPP than autocorrelation. Using broad geographic covariates would have likely decreased the proportion of variation explained by covariates and increased the proportion explained by autocorrelation. Data from CHaMP surveys are available across other sub-basins of the Columbia River, and there are a number of other large-scale monitoring efforts across this region (e.g., the Aquatic and Riparian Effectiveness Monitoring Program), making this kind of analysis broadly available.
Our results indicate that at unsampled sites, ecosystem processes can be predicted along with scaled estimates of error, which increases spatial resolution within a basin.

CONCLUSIONS
In this study, we explored linkages among ecosystem properties and metabolic processes, and we compared those patterns to the distribution of aquatic biota throughout a river network. Resolving linkages between abiotic conditions and biotic responses in GPP and ER is one incremental step in the larger task of understanding the complex patterns of metabolic regimes (Bernhardt et al. 2017) and the role of spatial patterns in rivers. Further, in our comparison of stream metabolism patterns with those of ESA-listed fish species, we take a first step in meeting the need for a spatially explicit understanding of the river network that considers the relationships between food webs, nutrients, and aquatic biota to better inform our understanding of streams at the network scale and ultimately stream management (Fausch et al. 2002, Naiman et al. 2012, Saunders et al. 2018. Our approach provides a framework that is applicable to large portions of the Columbia River basin as well as other areas where habitat monitoring is conducted at a comparable spatial scale.