Point-scale organic-matter decomposition in streambeds is weakly associated with reach-scale respiration

Stream and river ecosystems play a central role in the movement and decomposition of particulate organic matter, serving as a conduit between terrestrial hillslopes and coastal environments. Microbe-catalyzed decomposition generates simpler organic molecules that fuel respiration, often in the sediments of these ecosystems. However, the degree of connection between sediment-associated respiration (ER_sed) and organic-matter-decomposition potential remains poorly understood. It is also unclear whether organic-matter-decomposition potential is more closely associated with ER_sed, whole-ecosystem respiration (ER_tot), or water-column respiration (ER_wc). We examined the link between particulate organic-matter-decomposition potential – using cellulose-based cotton strips as a standardized substrate – and all three components of respiration across 48 sites in the environmentally diverse Yakima River Basin (Washington State, USA). We hypothesized that decomposition within sediments would be most strongly related to ER_sed, but decomposition rates were more closely associated with ER_tot, less so with ER_sed and not at all with ER_wc. This suggests that point-scale particulate organic-matter-decomposition potential within stream/river sediments is more closely associated with integrated system respiration rather than with processes confined to sediments or the water column alone though these relationships were weak overall. Further, across the basin, decomposition rates nearly spanned the previously reported global range for streams and rivers and were best explained by total dissolved nitrogen (TDN), sediment grain size, and aridity of the upstream drainage area. These results highlight the strong influence of land cover and basin-scale biophysical variation on sediment-associated decomposition processes and indicate that mechanistic models of organic-matter decomposition in streams and rivers should account for coupled sediment–water–land interactions.

1 Introduction

Stream networks are major components of the global carbon cycle (Cole et al., 2007; Drake et al., 2018; Talluto et al., 2024). Whole-stream metabolism is often studied as an integrated outcome of processes occurring across the continuum from the air-water interface down through sediments that are below the stream itself (Battin et al., 2023; Tank et al., 2010). The sediments include the interface between the streambed and water column (i.e., the benthic zone) and the spatial domain below this interface, referred to as the hyporheic zone (Boulton et al., 1998; Krause et al., 2011; Wondzell, 2011). The benthic zone can be highly productive with significant primary producer (e.g., algae) and heterotrophic microbial biomass (Allan et al., 2021). The hyporheic zone is also highly biogeochemically active due, in part, to surface water flowing through it and mixing with groundwater to stimulate heterotrophic microbial activity (Boano et al., 2014; Lewandowski et al., 2019; Zarnetske et al., 2011). Processes occurring on and around sediments across the benthic-to-hyporheic continuum are often jointly responsible for the bulk of biogeochemical activity in stream systems (Burrows et al., 2017; Fellows et al., 2001; Garayburu-Caruso et al., 2025; Naegeli and Uehlinger, 1997), with important exceptions in large rivers (Roley et al., 2023).

Metabolic processes within integrated stream systems are linked to the building of and breaking down of organic matter (Hall and Hotchkiss, 2017a; Odum, 1956). Streams are commonly net heterotrophic, meaning they mineralize more organic-matter-associated carbon than they fix by in-stream primary production (Battin et al., 2023; Bernhardt et al., 2022). This emphasizes the importance of understanding organic-matter decomposition in streams and its connection to respiration rates, which are ultimately linked to rates of elemental cycling. Organic-matter decomposition is commonly measured in stream systems by quantifying the breakdown of specific substrates (Benfield et al., 2017; Woodward et al., 2012). Cellulose-based cotton strips are an increasingly common model substrate for such studies as they enable broad comparisons across streams (and other types of ecosystems(Colas et al., 2019; Filbee-Dexter et al., 2022; Tiegs et al., 2024; Vyšná et al., 2014). Like other standardized substrates used in decomposition studies (e.g., litter bags), cotton-strips quantify organic-matter-decomposition potential, that is, the capacity of the ecosystem to decompose organic matter, rather than the actual rate of native organic-matter decomposition in situ. A standard approach is to place them in the field for a known amount of time, retrieve them, and measure the loss in tensile strength, as a measure of the degree of decomposition. This approach has revealed many factors that impact decomposition in streams such as temperature, land use, aqueous chemistry, sediment texture, stream flow, location within the stream network, and canopy cover, among others (Griffiths and Tiegs, 2016; Tiegs et al., 2024).

Studies have examined the relationship between organic-matter decomposition and whole-stream respiration (Mancuso et al., 2023; Pingram et al., 2020; Young and Collier, 2009), but have not specifically tied organic-matter decomposition within sediments to sediment-associated respiration. This leads to an open question and the focus of our study: To what degree are point-scale rates of organic-matter-decomposition potential within sediments linked to reach-scale respiration associated with the whole-stream ecosystem (ER_tot), the sediments (ER_sed), or the water column (ER_wc)? ER_tot represents reach-scale aerobic respiration from autotrophs and heterotrophs across benthic, planktonic, and hyporheic zones. ER_sed comprises reach-scale sediment-associated respiration from benthic and streambed sediments, rooted and submerged plants, and hyporheic zones that are hydrologically connected to the active channel; ER_wc constitutes reach-scale planktonic respiration occurring only in the water column.

We specifically tested the hypothesis that reach-scale ER_sed rates are strongly linked to point-scale measurements of organic-matter decomposition within streambed sediments. More specifically, we tested the prediction that cotton-strip-decay rates are best explained by ER_sed, with little additional variation in decay rates explained by ER_tot or ER_wc. To test our hypothesis and associated prediction we used field deployments across the Yakima River Basin (YRB). The YRB is an environmentally diverse basin in southeastern Washington State (USA) that is ∼16 000 km² and with a stream network that culminates in the 7th-order Yakima River. To generate the data needed to test our hypothesis, we used a combination of sensors and cotton strips across 48 sites in the YRB that collectively spanned a continuum from small mountainous streams in coniferous forests with little human impact to a large lowland river in an arid environment surrounded by significant agricultural land use. The resulting patterns help to fill a fundamental knowledge gap in our understanding of how point-scale organic-matter-decomposition potential relates to reach-scale respiration and can be used to inform models that aim to mechanistically integrate biogeochemical processes within and across stream networks.

2 Methods

To evaluate the linkages between organic-matter decomposition and stream ecosystem respiration we took advantage of a prior study (Garayburu-Caruso et al., 2025) that separated ER_tot, ER_sed, and ER_wc across environmentally divergent locations in the YRB (Fig. 1). Garayburu-Caruso et al. (2025) used dissolved oxygen (DO) sensors, dark bottle incubations, and the single-station method (Odum, 1956) to estimate these three components of respiration. In addition, that study deployed sensors to log water temperature, used here to calculate cumulative degree days as described below. These and associated contextual data were downloaded from existing data packages (Delgado et al., 2023; Forbes et al., 2023; Garayburu-Caruso et al., 2023), and methods are described in detail by Garayburu-Caruso et al. (2025). Data and computer scripts used in the current study are consolidated in Stegen et al. (2025). In brief, DO timeseries were analyzed via StreamMetabolizer (Appling et al., 2018) to estimate ER_tot. The length of the reach that a given DO sensor integrates varies across streams and has been estimated to be roughly three times the turnover length of O₂ (Hall and Hotchkiss, 2017b). As such, integrated reach lengths likely varied across sites due to differences in reaeration and discharge. Field sites were located in separate, non-overlapping reaches, helping to minimize potential spatial autocorrelation between neighboring sites in the estimated metabolic rates. To estimate ER_wc, 2 L opaque bottles containing a DO sensor were filled with stream water and incubated in situ; the rate of DO drawdown was used as the estimate of ER_wc. The difference between ER_tot and ER_wc was used as an estimate of ER_sed, which represents all respiration in the stream system that is not directly occurring in the water column. ER_sed therefore includes respiration in the hyporheic zone, the streambed surface, and rooted/submerged plants. We note, however, that ER estimates derived from open diel O₂ methods capture all oxygen-consuming processes, not solely aerobic respiration. These estimates could also include abiotic oxidation of dissolved Fe(II) or the oxidation of other end products from anaerobic metabolism (Demars et al., 2015; Piatka et al., 2021), other oxidation processes such as nitrification and photooxidation of organic matter (Demars et al., 2015; Estapa and Mayer, 2010), or O₂ inputs from groundwater (Hall and Tank, 2005). Because ER_sed was calculated by the difference between ER_tot and ER_wc, it may be sensitive to these non-respiratory O₂ sinks. However, biological metabolism is generally considered the dominant O₂-consuming pathway in streams and we therefore interpret ER_sed primarily as sediment-associated biological O₂ consumption while acknowledging that other contributions may be captured.

Figure 1 Biophysical clusters, sampling locations, and example conditions across the YRB. (a) The inset map shows the location of the YRB within the contiguous United States, with black indicating Washington State and red indicating the YRB. The YRB is shown with multiple colors, which correspond to biophysical clusters, as presented in Laan et al. (2025) and summarized briefly in the legend. Black circles are locations where decay rates were estimated. (b) Photos provide examples of the breadth of conditions studied across the YRB, post-incubation states of cotton strips, and deployment of the cotton strips.

To take advantage of the ecosystem respiration study, cotton strips made of Artist fabric (following the protocol of Tiegs et al., 2013) were deployed at the same time as the DO sensors. They were deployed upstream of the DO sensors to capture the upstream reach that influenced the DO sensor readings and prevent disturbance during sensor maintenance. This design prioritized among-reach coverage across the basin (48 sites) and provided point-scale estimates of decomposition potential at a single location within each reach.

The cotton strips were deployed in the shallow hyporheic zone for 35 continuous days at 48 sites across the YRB (Fig. 1). At each site, deployment locations were selected to be representative of the reach in terms of sediment size, flow velocity, and substrate composition. In larger river systems (e.g., 7th-order reaches), deployments were limited to wadeable areas near channel margins. Deployment and retrieval days varied across four days, but all strips were deployed for 35 d: deployments were from 25–28 July 2022 and retrievals were from 29 August–1 September 2022. Cotton strips were cut from bolts of 12-ounce, heavy-weight cotton fabric composed of 95 % cellulose (Style 548; Fredrix, Lawrenceville, GA, USA). Each strip was 27 threads wide and cut to 8.0 cm by 2.5 cm (after Tiegs et al., 2013). Each cotton strip was laid flat in a stainless-steel mesh cage (10.8 cm×4.5 cm, RSV Jumbo Mesh Herb Infuser) to minimize physical damage and feeding by macroinvertebrates, thereby emphasizing microbe-based decomposition. At each site, four cages with one strip each were attached to the underside of clay bricks ($\mathrm{20}\mathrm{cm}\times \mathrm{10}\mathrm{cm}\times \mathrm{5.5}\mathrm{cm}$) with stainless steel wire. The brick/cage/strip setup was nestled into streambed sediment such that the cages/strips were within the sediments and the brick was at the sediment/water interface. The four cages were next to each other. This setup kept the cotton strips out of direct light and within the sediments while allowing water to flow past the cotton strips.

After the 35 d incubation period, cotton strips were carefully removed from cages and gently brushed with gloved hands in stream water for approximately 10 s to remove large debris. Cleaned cotton strips were placed in 50 mL conical centrifuge tubes with 70 % ethanol. Tubes were capped and rolled approximately 10 times before ethanol was removed. After completing this step, clean 70 % ethanol was added to the 50 mL tube to minimize further microbial-based decomposition. Cotton strips were transported in the ethanol filled tubes on blue ice to Pacific Northwest National Laboratory in Richland, WA. At the laboratory, the ethanol was removed from the tubes and cotton strips were air dried overnight prior to further drying in an oven at 40 °C for at least 24 h. After drying, cotton strips were stored in air-tight containers with desiccant.

Dried cotton strips were shipped to Oakland University, Rochester, MI for tensile-strength analysis following the protocol in Tiegs et al. (2013). A tensiometer was used to estimate tensile strength (Mark 10 MG100 with a Chatillon TCM 201 with roller jaws). The tensiometer pulled each cotton strip at a rate of 2 cm min⁻¹. Some of the cotton strips were completely degraded such that there was no material to measure, while other cages only contained fragments that were too small to measure tensile strength. These non-detect strips accounted for 5.4 % of all deployed strips. In both cases, a limit of detection was assigned as the lowest tensile strength calculated in Tiegs et al. (2019) divided by 2, resulting in a final value of 0.05. This was done to avoid statistical artifacts that can arise when simply introducing a value of 0 (i.e., the natural log transformation in Eq. 1 would be undefined). Removing these strips entirely would bias the analyses away from conditions with very fast decomposition.

Tensiometer data were converted into decay rates using tensile loss calculated via Eq. (1) (as in Mancuso et al., 2022).

$$\begin{array}{}\text{(1)}& K={\displaystyle \frac{-\ln (T_{\mathrm{s}}/T_{\text{sc}})}{\text{Time}}}\end{array}$$

Here, K is the decomposition rate, T_s is the post-incubation tensile strength of the deployed cotton strips, and T_sc is the mean tensile strength of control strips that were not incubated in the field. The time variable was calculated as either the number of chronological deployment days (i.e., 35 d) or the number of degree days. Using degree days as the time variable accounts for variation in temperature across field sites and was estimated separately for each site as the sum of mean daily river temperature over the incubation period. We use K_cd and K_dd to represent decay rate per chronological day or per degree day, respectively. The values of K_cd and K_dd were estimated for each individual cotton strip and then replicates were averaged to provide a single site-level value for K_cd and K_dd.

We examined both K_cd and K_dd to evaluate whether the connection between decomposition and reach-scale respiration rates depends on accounting for temperature variation across the study basin. This is particularly relevant in the YRB because our field sites ranged from colder headwater streams to warmer low-gradient rivers. To test our hypothesis, we conducted both univariate and multivariate regression-based analyses. We used ordinary least squares regression to examine the strength of univariate correlations between K_cd or K_dd and ER_tot, ER_wc, or ER_sed. We complemented this univariate analysis with multiple regression analysis to find an optimized model to explain variation in either K_cd or K_dd.

Further, to explore how other system variables may explain variation in decomposition rates, a LASSO (Least Absolute Shrinkage and Selection Operator) regression model was built using physical, chemical and environmental variables (Table S1 in the Supplement) as inputs, and K_cd or K_dd as the response variables. Variables were cube-root transformed and z-score normalized to reduce the impact of high leverage points in the regression analysis and to equally weight all variables. The LASSO regression was performed over 100 iterations, each with a different random seed using the cv.glmnet function in the glmnet R package (Friedman et al., 2010). β coefficients were normalized to the maximum β coefficient for each iteration, then averaged over the 100 iterations for the final reported value. Both the raw and normalized mean β coefficient and standard deviation are reported in addition to the R² (Table 2).

3 Results and Discussion

3.1 Decomposition in the Yakima River Basin spans globally reported rates

Both K_cd and K_dd exhibited a wide range of values (Fig. 2), effectively spanning the theoretical maximum of what could have been observed with our deployment setup. This is evidenced by some cotton strips being completely consumed prior to retrieval (5.4 % of all strips, K_cd and K_dd maximized), while others were largely intact (K_cd and K_dd≈ 0). This variation is surprising given the relatively small spatial domain sampled by this study, and emphasizes that environmental heterogeneity can surpass the effects of spatial extent (Mancuso et al., 2022). The environments studied here ranged from pristine locations in the mountainous headwaters of the YRB to lowland locations with heavy agricultural influences (Fig. 1; Laan et al., 2025). This emphasizes the value of environmentally diverse watersheds like the YRB as useful testbeds to study variation in decomposition rates within single hydrologically connected basins.

Figure 2 Decay rate distributions and relationships to each other and temperature. Kernal density functions for (a) K_cd and (b) K_dd from a global streams and rivers dataset and from the YRB. (c) Scatterplot correlating K_cd to K_dd. Histograms summarizing the distribution of each rate are provided on the outer boundaries. (d) K_cd related to temperature summed across the deployment period; summed temperature was used to calculate K_dd. Blue lines represent Lowess spline fits as regression analysis was not required for interpretation.

Comparing decay rates from the YRB to a global dataset from >500 streams and rivers (Tiegs et al., 2024) showed that YRB rates spanned nearly the entire global range and had substantial overlap with the bulk of the global distributions (Fig. 2a and b). Dominant peaks in the YRB rate distributions were shifted slightly towards faster rates, relative to primary peaks in the global dataset (Fig. 2a and b). This shift towards faster rates and the wide range in rates are likely due to several factors linked to rates being estimated in later summer with high temperatures, and established biological communities due to several months since high-flow disturbances (Collier et al., 2013b; Grimm and Fisher, 1989; Mancuso et al., 2023). The two decay rates were also closely correlated with each other, though the relationship weakened towards locations with faster decay rates (Fig. 2c). This suggests a weak influence of temperature in the YRB; a strong influence of temperature should lead to a weak relationship between K_cd and K_dd. These results contrast previously reported influences of temperature on particulate organic-matter decomposition (Benbi et al., 2014; Griffiths and Tiegs, 2016), and likely reflect dominance of other influential factors such as variation in nutrient concentrations (Rosemond et al., 2015). Temperature-driven decomposition is also expected to lead to a strong relationship between K_cd and summed temperature, but we observed a very weak relationship despite ∼4-fold variation in summed temperature (Fig. 2d). This range in temperature among sites would likely be smaller in other seasons, and we do not therefore expect a strong influence of temperature to emerge in the YRB by conducting the study in other seasons. These results emphasize the need to understand factors governing variation in decay rates across the YRB. This is especially true given that rate distributions within this one basin span nearly all globally observed decay rates.

Across the YRB there appears to be potential for some spatial organization for both K_cd and K_dd (Fig. 3). Visual inspection of the maps suggests that the spatial organization may be stronger for K_cd than for K_dd. To evaluate this possibility more rigorously, we regressed each decay rate against upstream drainage area (Fig. 4). In this case, drainage area is meant to reflect position within the YRB. We used drainage area in preference to stream order because it is a continuous variable directly tied to the spatial domain a given stream integrates, whereas stream order is categorical and primarily reflects stream network topology. Associated regressions were significant (p<0.05) with both decay rates increasing with drainage area (Fig. 4). The relationship with drainage area was stronger, in terms of R², for K_cd. Both relationships were, however, relatively weak with R² values of 0.22 and 0.14 for K_cd and K_dd, respectively (Fig. 4). Nonetheless, the existence of a significant relationship after controlling for temperature (i.e., for K_dd) indicates that spatially structured factors other than temperature influence decay rates. This is not surprising as studies using cotton strips have found several factors that influence decomposition, such as nutrient concentrations, turbidity, and many others (Collier et al., 2013b; Pingram et al., 2020; Tiegs et al., 2024). Before exploring a broad suite of potential explanatory variables, we tested our hypothesis that decomposition rates will be better explained by sediment-associated respiration (ER_sed) than by respiration in the water column (ER_wc) or by respiration of the integrated stream system (ER_tot).

Figure 3 Spatial distribution of decay rates across the YRB. Each map shows elevation profiles and either K_cd (a) or K_dd (b). Colored circles are field locations where rates were estimated. The color of each circle is related to decay rate as indicated in the legends, and circle size is scaled to decay rate to further facilitate visual interpretability.

Figure 4 Decay rates increase with drainage area. K_cd (A) and K_dd (B) regressed against upstream drainage area and fit with an ordinary least squares linear regression. Associated regression models are shown as solid black lines and statistics are provided on each panel. Data were cube-root transformed to improve normality before analysis.

3.2 Respiration in sediments explains little variation in decomposition

Contrary to our hypothesis, we found that both decay rates were most strongly connected to ER_tot, less so with ER_sed, and not at all with ER_wc (Fig. 5). Univariate models using ER_tot were better than multivariate models using ER_sed, ER_wc, and their interaction. This is evidenced by univariate models using ER_tot having AIC values more than two units lower than multivariate models containing ER_sed and ER_wc (Table 1). Multivariate models were not used with ER_tot because it contains ER_sed and ER_wc. Garayburu-Caruso et al. (2025) showed that ER_sed accounted for the majority of ER_tot across most of the YRB, with 88 % of locations showing ER_sed contributions exceeding 50 % of ER_tot. The strong correspondence between these two metrics might explain the similar R² values observed. We note, however, that the difference in explained variance between ER_tot (R²=0.29 for K_cd; R²=0.16 for K_dd) and ER_sed (R²=0.22 for K_cd; R²=0.13 for K_dd) was modest, and the stronger association with ER_tot should be interpreted cautiously.

Figure 5 Decay rates related to each of the three aspects of stream ecosystem respiration. Both K_cd (A–C) and K_dd (D–F) show strongest relationships with ER_tot (A, D), weaker relationships with ER_sed (B, E), and non-significant relationships with ER_wc (C, F). Ordinary least squares linear regression models are shown with solid black lines and associated statistics are provided on each panel. All variables were cube-root transformed to improve normality prior to regression analysis.

Table 1 Comparison of univariate and multivariate regression models explaining variation in K_cd and K_dd. Model structures are indicated, along with change in AIC relative to the best model. ER_tot was not used in multivariate models because ${\text{ER}}_{\text{tot}}={\text{ER}}_{\text{sed}}+{\text{ER}}_{\text{wc}}$. Regression statistics for the univariate models are provided in Fig. 5; only the best performing univariate models, in terms of R², are shown. The models with ER_sed and ER_wc, but not the interaction term, are effectively the same as the ER_tot model because ${\text{ER}}_{\text{tot}}={\text{ER}}_{\text{sed}}+{\text{ER}}_{\text{wc}}$. They are included as a point of reference for the model that also includes the ER_sed×ER_wc interaction term.

To interpret these results, we note that cotton strips were deployed a few centimeters into the riverbed sediments, beneath a brick set flush with the streambed surface. We conceptualized this deployment as the shallow hyporheic zone, which was a reason we hypothesized that decay rates would be most strongly connected to ER_sed. However, the results indicate that point-scale decomposition of particulate organic matter within the shallow hyporheic zone is linked to respiratory processes occurring in both the sediment and water column. Therefore, our deployment depth might have reflected sediment-water interface processes in addition to shallow hyporheic zone processes, leading to stronger associations of cotton strip decomposition with ER_tot. We also note that ER_sed estimates may not exclusively reflect biological respiration, as non-respiratory O₂-consuming processes can contribute to these estimates (see Methods). We propose that if our deployment configuration was complemented with a simultaneous deployment that enabled growth of benthic algal biofilms on the cotton strips, the combined decomposition from both deployments could, in some systems, capture more of the processes that contribute to ER_sed. This would provide a more complete view of sediment-associated biogeochemical function, potentially leading to a stronger correlation between decomposition rates and ER_sed. While this remains to be tested, the underlying idea is that primary producers support a large portion of respiration associated with riverbed sediments, which is supported by recent analyses showing a strong link between ER_sed and gross primary production across the YRB (Garayburu-Caruso et al., 2025).

3.3 Slower rates of K_cd and K_dd are in streams with coarse sediments, set within wet forests

Given that our hypothesis was rejected and that ER_tot explained only about 29 % and 16 % of K_cd and K_dd, respectively (Table 1 and Fig. 5), we used a discovery-based approach to explore other system variables that may explain further variation in decomposition rates. LASSO-based modeling indicated that total dissolved nitrogen (TDN) and the median grain size of sediment texture (D50) were most important for explaining K_dd, while TDN and aridity were most important for explaining K_cd (Table 2). Other variables were retained in the LASSO models (Table S1), but we interpreted TDN, D50, and aridity as the most important because they consistently had the largest normalized coefficients. This interpretation is based on these variables having mean normalized coefficients above 0.5 – in terms of absolute value – meaning they were at least 50 % as important as the most important variable in the 100 LASSO model runs. Further, the LASSO coefficients for these variables had a coefficient of variation less than 0.5, meaning that across the 100 LASSO model runs the values of their normalized coefficients were relatively stable (Table S1). The LASSO modeling also confirmed a relatively weak influence of temperature, evidenced by relatively small and highly variable β coefficients for summed temperature in the K_cd model (Table S1); temperature was not used in the K_dd model.

Table 2 Regression coefficients from LASSO models explaining variation in K_dd and K_cd. Explanatory variables were cube-root transformed to reduce influence from high leverage data points and standardized as z-scores to enable direct comparison of the regression coefficients. LASSO models were fit over 100 seeds. Regression coefficients (β) and R² values were averaged across the 100 seeds. Normalized regression coefficients were calculated by dividing each β coefficient by the maximum β coefficient in each seed. Standard deviation (sd) is reported for each variables' coefficient over the 100 seeds and used to calculate the coefficient of variation (cv). Variables shown have an absolute value of mean normalized β of >0.5 and cv<0.5 to emphasize variables that were consistently important across seeds. Results for all variables, both normalized and not normalized, can be found in Table S1.

Both decomposition rates increased with higher TDN concentrations, while K_dd decreased with larger D50 and K_cd decreased with higher aridity index values (Table 2). To more deeply interpret these relationships, we examined Pearson-based univariate correlations between these three explanatory variables and other variables included in the LASSO models. This is important because of strong collinearity among some explanatory variables (Fig. S1 in the Supplement). In this case, variables identified as being the most important may be acting as proxies for one or more additional variables. We found that TDN was most strongly correlated with percent agricultural land cover of the upstream drainage area and sulfate concentrations in the stream water (Fig. S1). Increases in both decomposition rates with TDN may, therefore, reflect agricultural inputs of nutrients that increase overall microbial activity of the stream ecosystems we studied. D50 was most strongly correlated with the aridity index, which was most strongly correlated with percent forest cover; the correlation between D50 and aridity is unlikely to reflect a causal connection, while aridity and forest cover most likely are causally linked. If the relationship between decomposition and D50 is causal, it could be mediated by the total surface area available for microbial attachment. This would, however, influence decomposition only to the extent that microbial biomass in adjacent sediments impacts microbial activity on the cotton strips. Coarser sediments have much less surface area, potentially limiting overall microbial activity.

Considering the directionality of the univariate relationships, in context of the LASSO outcomes, suggests slower decomposition – for both rates – in streams with relatively coarse sediments and set within relatively wet forests. This contrasts with Clapcott and Barmuta (2010) who found faster decomposition in coarser sediments. The discrepancy is likely because we excluded macroinvertebrates while they did not, and they interpreted the link to sediment texture as due to greater habitat availability for macroinvertebrates in coarser sediments. Locations with slower decomposition should, thus, primarily be in higher elevation, relatively pristine parts of the YRB, while faster decomposition occurs at lower elevations impacted by agricultural inputs. These results are consistent with previous work showing greater cotton strip decomposition in impaired streams (Young and Collier, 2009), those with little natural land cover (Collier et al., 2013a; Webb et al., 2019), and in streams with higher nutrient concentrations (Ferreira et al., 2015; Pingram et al., 2020; Tiegs et al., 2013). In addition to differences in nutrient concentrations between higher and lower elevation sites, we expect less light penetration to streams in higher elevation sites because of more forest cover and smaller streams. Though not measured here, less light could suppress autotrophic production which may limit heterotrophic respiration (Bernhardt et al., 2022; Mulholland et al., 2001; Young and Huryn, 1999). Lower autotrophic production could therefore slow decomposition relative to high-light conditions by reducing the supply of labile carbon from phototrophs that can prime organic-matter degradation (Danger et al., 2013; Howard-Parker et al., 2020).

3.4 Point-scale decomposition is associated with processes across the sediment-water continuum and land features

Our results collectively indicate that to study shallow hyporheic zone decomposition processes, it is not sufficient to conceptualize organic-matter decomposition by only considering sediment or water column processes; one must examine the integrated system. The implication of our analyses is that point-scale organic-matter decomposition potential was more closely associated with integrated system respiration than with individual respiration components. Our findings suggest that mechanistic models of stream ecosystem respiration may benefit from accounting for sediment processes, water-column processes, and land-cover influences from beyond the stream. Focusing exclusively on the hyporheic zone is insufficient, even in systems for which the hyporheic zone accounts for most reactions (Boano et al., 2014; Burrows et al., 2017; McClain et al., 2003). This is further emphasized by previous work showing that most respiration occurs in the water column of large rivers (Gardner and Doyle, 2018; Roley et al., 2023). Garayburu-Caruso et al. (2025) also show that fractional contributions of ER_sed to ER_tot is often high, but that there is significant variation in those fractional influences across the YRB. This variability is due to ER_wc being fast enough, in some locations, to account for more than 80 % of ER_tot (Garayburu-Caruso et al., 2025). Similarly, Laan et al. (2025) found substantial overlap between the distribution of ER_wc rates from the YRB and ER_tot from across the contiguous United States. The overall picture is that decomposition is the result of integrated processes occurring across the sediment-water continuum and influenced by external factors tied to land cover and land use. We infer that these integrated processes are influenced by biophysical variation across the YRB (Laan et al., 2025), leading to decomposition rates within this single basin that resemble global rate distributions and nearly span the global range of observed rates (Tiegs et al., 2024). Other basins that contain only one ecoregion or have homogeneous land cover may be expected to have a narrower range of decomposition rates (Webb et al., 2019). We note that our decomposition estimates represent point-scale conditions at wadeable locations and should be interpreted as such rather than as reach-integrated measures of decomposition. Capturing within-reach spatial variability in decomposition across habitats and channel depths would complement the among-reach approach used here, and may lead to stronger relationships between decomposition rates and reach-scale ER. Nonetheless, models applied to any stream network that aim to predict spatiotemporal variation in decomposition rates would be well served by considering processes throughout the integrated watershed system.