A method for evaluating sediment-induced macroinvertebrate community composition changes in Idaho streams

BURP data

Benchmarks were developed using IDEQ Beneficial Use Reconnaissance Program (BURP) data. BURP is a bioassessment program that measures multiple stream physical, macroinvertebrate, and fish indicators in wadable 1^st–4^th order streams. IDEQ has collected BURP data in ~200 Idaho stream reaches during July-August each year in most years since 1993. BURP monitoring generally excludes federally designated wilderness areas, non-perennial streams, lake outlets, and reaches recently impacted by beaver activity or wildfire. Generally, BURP uses a targeted non-random sample design where crews select reaches representative of the associated stream assessment unit, the spatial scale used by IDEQ for CWA beneficial use support assessments. However, in several years an additional population of sites was also selected using a probability survey design for statewide data summaries. Data from both targeted and probability sampling were included in analyses, as described further below. BURP does not include repeat sampling of the same reach across years, but in some cases, sampling is spatially clustered across years and a single stream segment (National Hydrography Dataset Plus version 2 comid) was sampled multiple times.

All BURP data used in analyses were collected following protocols described in the Beneficial Use Reconnaissance Program Field Manual for Streams (Idaho Department of Environmental Quality (IDEQ), 2017). During each sample event field crews delineated a reach with length 30 times bankfull width or a minimum of 100 m. Macroinvertebrate samples were a composite of three samples, one from each of three separate riffle habitats within the reach. If three riffle habitats were not present, one or more samples were collected in a run. Samples were collected with a 500 µm mesh Hess sampler and composited and preserved with ≥70% ethanol in the field. Samples were identified to the lowest practical taxonomic level (typically to genus or species) and enumerated by EcoAnalysts Inc. following methods described in Richards et al. (2018). The same laboratory and methods were used across years. Substrate particle size was measured using modified Wolman pebble count methods (Wolman, 1954). Crews established three bankfull width transects, with each transect located one meter upstream of a macroinvertebrate sample location. Pebble counts were always conducted after macroinvertebrate samples to avoid disturbing macroinvertebrates during pebble counts. Within each transect at least 50 particles were selected at approximately equidistant intervals (heel to toe or one pace distance, etc.) and particle intermediate axis was measured and assigned to one of 11 size classes. Particles within both the wetted and non-wetted portion of the bankfull channel were included in counts. All particles within bankfull collected across the three transects were used to calculate percent surface fine sediment within the <2.5 mm size class. We assumed this size class represents sand and smaller particles and is functionally equivalent to surface fines measured as <2 mm. Hereafter we refer to percent surface fine sediment <2.5 mm as ‘SF’, shorthand for ‘sand and fines.’

Macroinvertebrate data were used to calculate fine sediment biotic index (FSBI) (Relyea, Minshall & Danehy, 2012) scores. FSBI is an index for assessing surface fine sediment impacts on macroinvertebrate community composition in northwest US streams. Relyea, Minshall & Danehy (2012) used macroinvertebrate data from multiple monitoring programs in the Pacific Northwest, including BURP data, to identify 206 taxa that are common in the Pacific Northwest (occurring in ≥2% of streams in their database). They assigned each taxon a sediment sensitivity score based on the percent fines <2 mm associated with the taxa’s 75^th percentile occurrence. FSBI scores for a sample are then calculated as the sum of taxa-level scores (Relyea, Minshall & Danehy, 2012). FSBI is thus based on taxa occurrence only. FSBI scores generally range from 0–350; high scores indicate many sediment-sensitive taxa are present and therefore elevated surface fine sediment levels likely have not impacted macroinvertebrate community composition. Low FSBI scores indicate few sediment-sensitive taxa are present. FSBI has been applied in state-wide stream condition assessments in Washington State (Larson et al., 2019) and for assessing effects of streambed instability on macroinvertebrates (Kusnierz, Holbrook & Feldman, 2015). FSBI was used rather than other macroinvertebrate indicators because FSBI is sediment-specific and is correlated with more generalized macroinvertebrate indicators such as ephemeroptera, plecoptera, and tricoptera (EPT) taxa richness (Larson et al., 2019). In preliminary analyses there also was a stronger stressor-response relationship between SF and FSBI than with EPT taxa richness.

The ‘reference’ and ‘stress’ site definitions and datasets used by IDEQ to develop BURP bioassessment metrics were used for analyses (Jessup, 2011; Jessup & Pappani, 2015; Idaho Department of Environmental Quality (IDEQ), 2016). These reference and stress site definitions are routinely applied by IDEQ for CWA applications in Idaho (Idaho Department of Environmental Quality (IDEQ), 2016). BURP reference and stress sites were defined as the 10% least and 10% most-disturbed 1998–2007 BURP sites within each of three site classes, based on nine landscape indicators of upstream human landscape disturbance (Jessup & Pappani, 2015). Site classes are ecoregion groupings that help explain macroinvertebrate community structure variability within reference sites (Fig. 1). They were defined by applying non-metric multidimensional scaling and principal components analysis to Idaho macroinvertebrate data within BURP reference sites (Jessup, 2011; Jessup & Pappani, 2015). Site classes are described further in supplemental information. In cases where there were multiple BURP reference or stress sites within a stream segment (comid) across years, one reference or stress site per comid was randomly selected for use in analyses. Summary statistics for BURP reference and stress datasets used for analyses are in Table 1 for SF and Table 2 for FSBI.

Reference benchmarks

Reach-specific SF reference benchmarks (SF_ref) were developed based on relationships between SF and bankfull width within each stream order using data from BURP reference sites (Fig. 2). Strahler stream order used for benchmark development was based on National Hydrography Dataset Plus version 2 1:100,000 scale hydrography used by Idaho for CWA applications. Bankfull width and order both had strong correlations with SF (see supplemental information). SF decreased with order within 1^st–4^th order streams, and order appeared to be a reasonable categorical proxy for stream power and several other physical variables correlated with SF. Relyea, Minshall & Danehy (2012) observed surface fines decreased with stream order using a larger Northwest U.S. dataset. Miller et al. (2021) also observed strong relationships between bankfull width and SF in Idaho streams, and divided streams into groupings based on bankfull width to calculate SF reference benchmarks. Quantile regression was used to define SF_ref, the 75^th percentile value expected among reference streams based on reach order and bankfull width (Fig. 2). All BURP sites that met reference criteria, including sites sampled using a targeted (N = 277) and probability survey design (N = 17) were used in regressions. We used quantile regression because previous studies indicated trying to predict site-specific reference fines condition exactly can result in order of magnitude or larger site-specific prediction errors (Hawkins, Olson & Hill, 2010; Grangeon et al., 2023). Local geomorphic features, within-reach processes, and local hill slope processes can have significant impacts on reach bed sediment conditions but can be challenging to capture in models (Hawkins, Olson & Hill, 2010; Snyder et al., 2013; Keestra et al., 2018). Using quantile regression allowed us to acknowledge this local variability. The 75^th percentile was selected to represent an upper bound SF value expected under reference conditions. This approach is conservative (protective) because 25% of reference sites have SF values exceeding this threshold. The ‘quantreg’ R package (version 5.97) was used for regressions (Koenker, 2023). FSBI reference benchmarks (FSBI_ref) were calculated as the 25^th percentile FSBI value within each BURP site class described above.

Stressor-response benchmarks

Logistic regression was used to predict SF benchmarks with a 50% (SR50) and 75% (SR75) probability that FSBI would be worse (less) than FSBI_ref. Unlike reference benchmarks, stressor response benchmarks were not reach-specific. Separate logistic regressions and benchmarks were developed for each stream order/site class combination. Regressions used all BURP data collected 1998–2021 (reference, non-reference, stress collected using either targeted or probability survey design) with SF as the predictor variable and a binary value (1 or 0) indicating if FSBI was less than FSBI_ref as the response variable. Logistic regressions were implemented using the ‘glm’ function in base R. Regression model fit was documented using the model chi square statistic, odds ratio confidence intervals, and Hosmer-Lemeshow (HL) goodness-of-fit test (Hosmer et al., 1997; Hosmer, Lemeshow & Sturdivant, 2013). The chi square statistic quantifies the difference between deviance associated with a null model with an intercept only and model deviance. A large and statistically significant (<0.05) chi square statistic provides evidence to reject the null hypothesis that the model does not perform better than random chance. The model odds ratio indicates how the odds of classifying a reach as having FSBI worse than reference increases for a 1% increase in SF. When lower (2.5%) and upper (97.5%) confidence intervals for model odds ratio values both exceed one, this provides evidence that increasing SF reduces FSBI. The HL goodness-of-fit test describes the level of agreement between observed and predicted outcomes by decile of predicted probability. The HL null hypothesis is that the model provides a good fit. Low HL test p values provide evidence for rejecting the null hypothesis and concluding model fit is poor. Statistics and logistic regression curve plots were used in a weight of evidence approach to evaluate whether developed regression models were adequate for purposes of defining stressor-response benchmarks.

Relative risk

Relative risk (RR) (Van Sickle et al., 2006) was calculated to describe the strength of association between poor SF condition and poor macroinvertebrate condition. RR was the ratio of two conditional probabilities: the probability that a macroinvertebrate indicator benchmark was exceeded when a SF benchmark was exceeded, divided by the probability that a macroinvertebrate indicator benchmark was exceeded when a SF benchmark was not exceeded. RR was calculated for SF_ref and each stressor-response benchmark using data from two probability survey datasets. First, RR was calculating using data from 85 1^st–4^th order stream reaches sampled across Idaho 2013–2016 using BURP methods and a spatially balanced probability survey design (Idaho Department of Environmental Quality (IDEQ), 2018). Site inclusion probabilities were based on stream order and total state stream length within each order. These BURP probability survey data were not used in reference benchmark development but were used for stressor-response benchmarks. Second, RR was calculated using 43 stream reaches sampled across Idaho Bureau of Land Management (BLM) lands 2013–2016 using a probability survey design and BLM Aquatic Inventory and Monitoring (AIM) program protocols (Bureau of Land Management, 2021; Miller et al., 2021). Site inclusion probabilities were based on stream order categories (1^st and 2^nd order, 3^rd and 4^th order, and 5^th + order) and stream linear extent within each category (Miller et al., 2021). AIM collects paired Wolman pebble and macroinvertebrate data within each reach but uses a different sample design and field methods from BURP. Compared to BURP, more transects are used per reach, more macroinvertebrate subsamples are collected, and different macroinvertebrate collection methods are used, among other differences. A table comparing monitoring programs is included in supplemental information. Calculating RR using AIM data tests if associations between SF benchmark exceedance and macroinvertebrate effects persist across different monitoring approaches.

RR was calculated using SF as the stressor indicator and FSBI as the response variable. SF and FSBI benchmark status was used to rate SF and FSBI as either ‘good’ or ‘poor’ for each sampled reach. In addition, RR was calculated using SF as the stressor indicator and other macroinvertebrate indicators as the response variable. For BURP data, IDEQ’s multi-metric macroinvertebrate index (SMI2) was used. For AIM data, an index describing the reach-specific ratio of observed to expected (O/E) macroinvertebrate taxa (Miller et al., 2021) was used. FSBI is an index based on taxa sediment tolerances, whereas SMI2 and O/E are generalized indicators of macroinvertebrate community condition that are not specific to any pollutant. O/E describes the proportion of macroinvertebrate taxa expected under reach-specific reference conditions that were observed in samples. SMI2 reflects the macroinvertebrate community expected under BURP reference condition based on multiple individual macroinvertebrate metrics (Idaho Department of Environmental Quality (IDEQ), 2016). Program-specific SMI2 (≥52–54 depending on site class) and O/E (≥0.63) benchmarks were used to identify good and poor macroinvertebrate condition. For each dataset, RR values and an associated 95% confidence interval were calculated using the ‘spsurvey’ R package (Dumelle et al., 2023). Calculations used site weights (i.e., inverse of site inclusion probability) to generate unbiased RR estimates. RR values with a 95% confidence interval >1 indicate a risk of poor macroinvertebrate indicator condition when the SF benchmark is exceeded. Higher RR values suggest a stronger association between stressor and response indicators.

Application framework

We developed a simple framework for applying SF and FSBI benchmarks to rate the strength of evidence for a sediment-induced change in macroinvertebrate community composition for each sample event with paired data. The framework rates a sediment effect as ‘unlikely’ if both SF_ref and FSBI_ref benchmarks are achieved, as having ‘mixed evidence’ where only one reference benchmark is achieved, and as ‘likely’ if both reference benchmarks are not achieved. The framework was designed to be consistent with the desire for at least two lines of evidence for a sediment effect when evaluating compliance with Idaho’s narrative sediment criterion (Idaho Department of Environmental Quality (IDEQ), 2016).

We evaluated the framework by applying it to BURP, AIM, and PacFish/InFish Biological Opinion (PIBO) monitoring program (Kershner et al., 2004; Roper, Saunders & Ojala, 2019; Saunders et al., 2022) data collected throughout Idaho. PIBO data were not included in RR calculations because the PIBO sample design was not suitable for RR calculations. BURP, AIM, and PIBO sample events where paired Wolman pebble and macroinvertebrate samples were collected were included in analyses. Because some AIM and PIBO sites were sampled multiple times across years, we randomly selected one sample event per AIM and PIBO site to reduce spatial sampling bias among data used for analyses. Datasets used for evaluation included 1998–2021 BURP data (N = 5,215 sites/events), 2013–2022 BLM AIM data (N = 299 sites/events), and 2004–2019 PIBO data (N = 665 sites/events). Notched boxplots were used to examine how macroinvertebrate indicator (FSBI, SMI2, O/E) distributions varied with benchmark status and framework prediction classes. Statistically significant differences in median macroinvertebrate indicator values between condition categories were inferred by examining overlap (or lack thereof) between boxplot notches representing median 95% confidence intervals. All data and R code used for analyses are available online (https://doi.org/10.17605/OSF.IO/3cZRP). DEQ BURP data used were queried from DEQ’s internal BURP database. AIM and PIBO data were provided by USFS and BLM staff.

Application considerations

The benchmarks and framework developed here could be used to assist stressor identification. The CWA requires states to identify stressor(s) impairing aquatic life and develop a restoration plan to reduce stressor levels. When macroinvertebrate samples suggest the community is not within the “natural range of reference streams or conditions” and thus does not meet the definition of full support of aquatic life use in Idaho WQS (IDAPA 58.01.02.010.42), methods to diagnose the specific stressor(s) causing impairment are needed. Our benchmarks and framework could be used to assist stressor identification but should not be used as the only evidence diagnosing sediment impairment. BURP, AIM, and PIBO bioassessment programs all collect paired Wolman pebble and macroinvertebrate data, but also measure many other parameters relevant to assessing sediment-induced macroinvertebrate effects. The parameters measured differ across programs but include various stream physical habitat integrity parameters and generalized macroinvertebrate metrics. These other data should also be considered when diagnosing sediment impairments. Our framework also does not provide information on potential human sources or transport pathways, which are needed to assess compliance with Idaho’s narrative sediment criterion (Idaho Department of Environmental Quality (IDEQ), 2016). We recommend assessors combine benchmark and framework predictions with other lines of evidence in a weight of evidence approach, and explicitly articulate the lines of evidence and rationale for sediment assessment decisions.

When sediment has been confirmed as a stressor, the CWA requires states to define sediment targets to restore support of aquatic life use. Our benchmarks could be used as sediment restoration targets in streams where macroinvertebrate effects through substrate exposure pathways are the primary concern. Our approach does not address other sediment exposure pathways or ecological effects. Idaho WQS include numeric turbidity criteria to protect aquatic life from water column suspended sediment exposure pathways (IDAPA 58.01.02.250.02.3.e, IDAPA 58.01.02.401.02a). IDEQ guidance for selecting sediment targets also includes suggestions for suspended sediment concentration targets and substrate (surface plus subsurface) sediment targets based on streambed core samples for protecting salmonid spawning habitat (Idaho Department of Environmental Quality (IDEQ), 2003). IDEQ sediment target guidance did not include suggestions for surface fine sediment targets. This study represents one potential approach for selecting surface fine sediment targets.

Our benchmarks and framework predictions are specific to a single sample event within a single stream reach. Assessors conducting stressor identification analyses will likely face data from multiple sampled reaches within an area of interest or sample reaches with multiple sample events across years. Appropriate methods for addressing temporal variability or reconciling conflicting outcome predictions within or across reaches are application-specific and beyond the scope of this study. As described above, we recommend considering framework predictions along with other relevant lines of evidence relevant to a sampled reach. What lines of evidence are available will vary geographically, so appropriate decision rules will also. Management goals, such as the biological endpoint (macroinvertebrates, fish, endangered species, etc.) targeted often varies geographically, affecting decision rules needed.

Our SF reference benchmarks are based on channel bankfull width and implicitly assume reach bankfull width is within reference condition. In cases where human activities have increased bankfull width, we still expect the benchmark and framework to identify sediment-induced macroinvertebrate community changes. SF reference benchmarks decrease as bankfull width increases (Fig. 2) because in reference streams SF decreases with stream power. Generally, when human activities increase channel width stream power decreases and substrate fine sediment increases. In such cases, measured SF would likely exceed our SF reference benchmarks. Patterns in Figs. 4–6 suggest relationships between reference benchmark exceedance and macroinvertebrate responses are robust. The benchmarks and framework are not a tool for identifying stream channel geomorphic changes from reference condition but may complement other tools for evaluating anthropogenic channel alterations.