Linking stream ecosystem integrity to catchment and reach conditions in an intensively managed forest landscape

. Forests are vital to maintaining headwater stream integrity in forested biomes, which ensures the delivery of aquatic ecosystem services downstream. Forest harvesting, however, can alter land – water linkages and compromise stream integrity. Historically, the main effects of forestry on streams have been documented by studies that used relatively few (mainly abiotic) indicators and which focused on single harvesting events. However, forest management is expected to intensify in the future to meet increasing global wood demand and it is likely that our present understanding does not adequately capture the cumulative effects that streams will be subjected to under intensive forest management. To address this, we assessed the effects of varying forest management intensities on the integrity of 15 forest headwater streams in northwestern New Brunswick, Canada. We used a comprehensive approach to link multiple biotic and abiotic indicators of stream ecosystem integrity to reach- and catchment-level characteristics including forest management (e.g., cumulative harvesting over time, road density, forest condition). Most indicators detected the gradient in forest management intensity with abiotic indicators responding most strongly. Streams in catchments with highest management intensity (especially road density) tended to have higher ﬁ ne inorganic sediment deposition and entrainment, water cations and carbon, dissolved organic matter humi ﬁ cation, and water temperature. These abiotic differences were associated with higher bio ﬁ lm biomass and shredder densities, but lower leaf decomposition. Evidence from our multi-indicator approach elucidated a potential effects pathway of higher inorganic sediment content in bio ﬁ lms of organic matter potentially limiting or altering its use by microbial and benthic macroinvertebrate (BMI) communities and resulting in reduced leaf decomposition rates. Overall, this study shows that current best management practices in an intensively managed watershed (and legacy effects from past management such as older road systems) do not fully protect against an increased delivery of terrestrial materials to streams with resulting habitat and biotic changes, but that they are mostly effective at preventing the impairment of BMI communities.

the delivery of aquatic ecosystem services downstream. Forest harvesting, however, can alter land-water linkages and compromise stream integrity. Historically, the main effects of forestry on streams have been documented by studies that used relatively few (mainly abiotic) indicators and which focused on single harvesting events. However, forest management is expected to intensify in the future to meet increasing global wood demand and it is likely that our present understanding does not adequately capture the cumulative effects that streams will be subjected to under intensive forest management. To address this, we assessed the effects of varying forest management intensities on the integrity of 15 forest headwater streams in northwestern New Brunswick, Canada. We used a comprehensive approach to link multiple biotic and abiotic indicators of stream ecosystem integrity to reach-and catchment-level characteristics including forest management (e.g., cumulative harvesting over time, road density, forest condition). Most indicators detected the gradient in forest management intensity with abiotic indicators responding most strongly. Streams in catchments with highest management intensity (especially road density) tended to have higher fine inorganic sediment deposition and entrainment, water cations and carbon, dissolved organic matter humification, and water temperature. These abiotic differences were associated with higher biofilm biomass and shredder densities, but lower leaf decomposition. Evidence from our multi-indicator approach elucidated a potential effects pathway of higher inorganic sediment content in biofilms of organic matter potentially limiting or altering its use by microbial and benthic macroinvertebrate (BMI) communities and resulting in reduced leaf decomposition rates. Overall, this study shows that current best management practices in an intensively managed watershed (and legacy effects from past management such as older road systems) do not fully protect against an increased delivery of terrestrial materials to streams with resulting habitat and biotic changes, but that they are mostly effective at preventing the impairment of BMI communities.

INTRODUCTION
Among the many ecosystem services provided by forests, aquatic ecosystem services are some of the most valuable. For example, forests filter pollutants before they reach waterbodies and stabilize soil and stream banks to mitigate erosion (Furniss et al. 2010). Ecosystem services such as food provision, recreation, and hydrologic control (e.g., flood control) are also determined by the degree of control that forests exert on energy and water fluxes, such as subsidization of stream food webs, or rain interception, evapotranspiration, and infiltration (Wells et al. 2010). However, because forest streams are intimately linked ecologically to their surrounding landscape (Hynes 1975), disturbances to forests such as harvesting may interfere with ecological functions and jeopardize important aquatic ecosystem services (Pohjanmies et al. 2017). For example, forest harvesting can result in greater hydrological variability (e.g., higher peak flows) for several years after harvesting due to lower infiltration rates and/or reduced precipitation, interception, and evapotranspiration Wondzell 2005, Buttle et al. 2009), affecting the service of flood control. Forest harvesting can increase suspended solids and nitrates (Croke and Hairsine 2006, Kreutzweiser et al. 2008b, Richardson and B eraud 2014, Webster et al. 2015, potentially compromising the provision of clean water; increase stream water temperature and thermal diel fluctuations ; and decrease inputs of organic matter and woody debris to streams (Bilby and Ward 1991, England and Rosemond 2004, Santiago et al. 2011, Burton et al. 2016. This may cause direct and/or indirect changes in the composition and structure of stream biota (Mellina and Hinch 2009, Kreutzweiser et al. 2013, Richardson and B eraud 2014, leading to altered ecosystem functions such as leaf litter decomposition and primary production (Kiffney et al. 2003, Mckie and Malmqvist 2009, Kreutzweiser et al. 2010, Yeung et al. 2017, and, in turn, affect the delivery of aquatic ecosystem services (Balvanera et al. 2006, Woodward 2009).
To minimize impacts, the forest industry has implemented best management practices (BMPs), including riparian buffer zones (i.e., streamside restricted-harvest forest reserves) and guidelines for stream crossings and road construction (Schilling 2009, McDermott et al. 2010, OMNR 2010. Their effectiveness has been assessed (Broadmeadow andNisbet 2004, Cristan et al. 2016), but most have focused on effects of single harvesting events (e.g., Jackson et al. 2007, Witt et al. 2016, Yeung et al. 2017). However, disturbances can happen at different times and locations in the catchment and with different harvesting techniques (e.g., clearcutting or thinning) and may be compounded by the effects of roads, site preparation, herbicide application, etc., resulting in a high potential for cumulative impacts of forest management on aquatic ecosystems (Kreutzweiser et al. 2013, Webster et al. 2015). Yet, comparatively little is known of the cumulative effects of intensive forest management on stream ecosystems, which could be a critical information gap considering that the industry is expected to intensify in the future (Creed et al. 2016). Therefore, we sought to understand whether existing BMPs are effective at protecting streams and their associated aquatic ecosystem services from the cumulative effects associated with this intensification by conducting our study in an intensively managed forest watershed.
Most studies have assessed the effects of forestry on stream ecosystems with only a limited set of biotic and/or abiotic indicators (Janisch et al. 2012, Oliveira et al. 2016, Witt et al. 2016, with few incorporating indicators addressing stream ecosystem function (Cristan et al. 2016). However, each indicator may respond differently to catchment disturbance and interact in unpredictable ways (e.g., synergistic/antagonistically; Piggott et al. 2015, Nõges et al. 2016. To understand the often complex interactions and responses, holistic approaches that incorporate multiple structural and functional components of the ecosystem are needed (Gessner andChauvet 2002, Parr et al. 2016); this could improve our understanding of the equivocal responses to forest harvesting in the literature (Richardson and B eraud 2014).
To address these two knowledge gaps, we used a multi-indicator approach to assess the condition of 15 streams with catchments ranging in forest management intensity in New Brunswick, Canada, in 2014 and 2015. Twelve of the streams were located in the Black Brook (BB) forestry district, because (1) BB is considered one of the most intensively managed forests in Canada and it is third-party certified as sustainably managed under the Sustainable Forestry Initiative (Etheridge et al. 2005), and (2) state-of-the science forest characterization and harvest data were available from high-resolution light detection and ranging (LiDAR) and other advanced remote sensing and GIS techniques, enabling a more refined characterization of landscape and forest characteristics. Additionally, three unmanaged catchments were located in Mount Carleton (MC) Provincial Park, providing a regional range of natural variation in the indicators. The objective of this paper was to understand how different intensities of forest management and the derived catchment and reach variables affect stream ecosystem integrity. This objective was addressed by measuring multiple physical (water chemistry, water temperature, and sediment deposition), structural (algae and Benthic macroinvertebrate [BMI] communities), and functional (leaf litter decomposition) indicators of stream ecological integrity across a gradient of forest management intensity in a twoyear observational study.

Site selection and characterization
The study was conducted at two locations in northwestern New Brunswick, Canada, within the Atlantic Maritime ecozone. The 190,000-ha Black Brook District (BB), owned by J.D. Irving, Ltd. (JDI), is one of the most intensively managed forest lands in Canada (Etheridge et al. 2005), and MC, the reference area, is a Provincial Park~100 km east of BB with no harvesting history (Fig. 1). Black Brook is part of the Sisson Ecodistrict in the Central Uplands Ecoregion, and it is composed mainly of Ordovician-Devonian sedimentary rocks, and deep, loamy soils in areas of low relief and less fertile, shallow, stony soils in higher terrain (Zelazny 2007a). Black Brook is composed of shade-tolerant hardwood stands (25%), mixedwood (18%), softwood-cedar (15%), and softwood (42%) forests. Eighty-eight percent of the softwood forest was comprised of plantations in 2002 including 56% black spruce (Picea mariana), 30% white spruce, 9% Norway spruce (Picea abies), and 3% pine (Pinus sp.; Etheridge et al. 2006). Mount Carleton is part of the Ganong Ecodistrict in the Highlands Ecoregion, and it is underlain by Devonian felsic volcanic rocks, which result in coarse-textured, poor soils, but are accompanied by fine-textured particles derived from metasedimentary rocks (Zelazny 2007b). Balsam fir and black spruce dominate the vegetation in this ecodistrict, but in mid-slopes white spruce, white birch (Betula papyrifera), yellow birch, and red maple (Acer rubrum) are common as well (Zelazny 2007b).
Twelve low-order streams (1st-3rd) were selected in the northern end of BB ( Fig. 1), which ranged not only in harvesting intensity (18-100% of the catchment harvested in the last 10 yr), but also in road density (21-89 m of road per ha), stream crossings (0-4), catchment size (77-389 ha), and forest composition (deciduous/ mixed/coniferous dominated; Appendix S1: Table S1). Mount Carleton was selected as the reference area despite being part of a different ecoregion because there was an absence of unmanaged catchments closer to BB. The reference stream MC2 had a beaver dam~1 km upstream of the sampling reach, which presented an opportunity to compare anthropogenic (harvesting in BB) to natural (beaver dam in MC2) disturbance.
Catchments in BB were delineated using an ArcGIS (ESRI 2011) flow accumulation grid based on a 5-m LiDAR-derived digital elevation model (DEM) obtained from JDI, whereas catchments in MC were delineated by using the provincial 20-m DEM. Detailed characterization of each BB catchment was achieved using LiDAR-derived forest structural metrics, photo-interpreted composition data from high-resolution digital stereo imagery (Forest Resource Inventory or FRI), and topographic information derived from the LiDARbased DEM (see Catchment and reach-level explanatory variables for more details). High-resolution LiDAR was not available for MC catchments, with implications for analyses (see Statistical analysis).
Sixty meter long stream reaches were selected for sampling immediately upstream of road crossings. However, for eight out of 15 streams, roads intersected stream channels further upstream (at least 300 m) of the sampling reach. Road crossings were therefore included as a potential stressor contributing to the cumulative effects from whole-catchment forest management practices (see Landscape features).

Measurement of response variables
Sediment deposition.-Seven sediment traps (50-mL polypropylene centrifuge tubes placed in bricks; Kreutzweiser et al. 2009) were deployed on the stream bottom of depositional areas along the study reach. After 23-24 d in 2014 (October-November) and 25-30 d in 2015 (September-October), caps were replaced and tubes withdrawn and stored frozen.
Contents of the thawed tube were poured through a 250-lm sieve and the filtrate retained. The filtrate was filtered through pre-ashed 1.2-lm GF/C Whatman filters (Whatman, Little Chalfont, UK), and filters dried at 60°C for 48 h and ashed at 500°C for 2 h to calculate ash-free dry mass (AFDM). The size class used for this study (1.2-250 lm) is representative of clay and silt particles known to be transported from roads or harvested sites and which pose biological problems (e.g., clogging gills; Waters 1995).
Water chemistry and dissolved organic matter (DOM) quality.-Water samples were collected at the downstream end of each sampling reach in September, October, and November of 2014 and in August, September, and October of 2015. Samples were kept refrigerated and in the dark until analyzed in the laboratory.
Water chemistry variables were measured at the Great Lakes Forestry Centre (Sault Ste. Marie, Ontario, Canada) following standard methods (Hazlett et al. 2008) for pH, conductivity, alkalinity, cations and anions (Ca, K, Mg, Na, SO 4 , Cl, SiO 2 ), nutrients (NO 2 +NO 3 , NH 4 , total N, dissolved organic carbon, DOC; dissolved inorganic carbon, DIC; reactive P, total P), and metals (Al, Fe, Mn, Zn, Cd, Cu, Ni, Pb). Principal components analysis (PCA) was conducted with the mean value of water chemistry variables for each year, and the first and second principal components (described in Results: Water chemistry and DOM quality) were used to represent the collective water chemistry variables in subsequent statistical analyses.
Water samples for DOM quality were filtered at 0.2 lm and characterized at Laurentian ❖ www.esajournals.org University (Sudbury, Ontario, Canada) using Cary Eclipse (Varian Instruments, Walnut Creek, California, USA) and Cary 60 UV-Vis spectrophotometers (Agilent Technologies, Santa Clara, California, USA). Three-dimensional fluorescence scans were run at 5-nm excitation steps from 250 to 450 nm, and emissions were read at 2-nm steps from 300 to 600 nm. The generated excitation-emission matrices were then corrected and adjusted, and variables describing optical properties of DOM were calculated: HIX, an indicator of the humification degree of DOM (calculated following Zsolnay et al. 1999); SUVA, an indicator of DOM aromaticity (Weishaar et al. 2003); E2E3, negatively related to DOM molecular weight (Peuravuori and Pihlaja 1997); the Fluorescence Index, an indicator of DOM origin (terrestrial vs. microbial;McKnight et al. 2001); and the freshness index (BA; Parlanti et al. 2000).
Water temperature.-Water temperature was continuously measured during the sampling season (31 July-10 November in 2014; 15 June-22 October in 2015) with water temperature and level data loggers (Onset HOBO Data Loggers or Solinst Leveloggers), and daily maximum, minimum, and mean monthly, summer, and fall average temperatures were calculated. The relative trends across streams were comparable for the different temperature variables (data not shown), so only averaged daily maximum August temperature (1-31 August) is presented herein.
Biofilm and algal biomass on tiles.-To estimate the biomass of periphyton communities and total biofilm (algae, bacteria, fungal biomass) on rock surfaces in streams, one row of 10 unglazed clay tiles (4.7 9 4.7 cm) was deployed in each of three different riffles within the 60-m reach. The 10 tiles were glued to a band, which was secured to the streambed using rebar. Five of the tiles in each row were pooled to measure total biofilm biomass, and the other five were pooled to quantify chlorophyll a content as an indicator of periphyton biomass. In 2014, the incubation started at the end of July/beginning of August and lasted 48-51 d. In 2015, tiles were deployed in mid-July and incubated for 24-25 d. The incubation was shortened in 2015 because in 2014 grazers (notably Glossosoma) were found on some of the tiles; the shortened period ensured sufficient time for biofilm growth, while minimizing grazers colonization of the tiles. Upon retrieval, tiles were scraped with scalpels and the slurry washed into Whirl-Pak bags. Samples were kept in the dark and cool in the field and frozen at the end of the day.
In the laboratory, samples were thawed and filtered through pre-combusted and pre-weighed GF/ C Whatman filters. Filters for chlorophyll a analysis were frozen, and filters for total biofilm biomass were dried for 48 h at 60°C and ashed at 500°C for 2 h to calculate AFDM. Chlorophyll a was extracted by submerging the filters in 90% ethanol in glass tubes placed in an 80°C water bath for 7 min, and chlorophyll a concentration of the extract was measured using a Trilogy fluorometer (Turner Designs, Sunnyvale, California, USA). An autotrophic index was calculated by dividing chlorophyll a concentration by total biofilm biomass.
Leaf decomposition and benthic macroinvertebrate communities.-Senescent speckled alder leaves (Alnus incana ((L.) Moench)) were pre-leached for 48 h and dried at 30°C for 48 h. Leaves were weighed into 4.0 AE 0.1 g groups and placed in mesh bags with a 5 9 10 mm mesh size. Leaf packs were tied to bricks placed on the stream bottom or tied to rebar driven into the substrate where water levels were too low to use bricks. Six leaf packs were distributed within each 60-m reach at each site and incubated for 31-33 and 33-35 d in September 2014 and 2015, respectively. Leaf packs were discarded from one of the streams (BB05) in 2014 due to low water levels during the incubation period. Leaf packs were retrieved and the contents preserved in 37% formaldehyde (~10% of the volume).
In the laboratory, containers were emptied into a 250-lm sieve and leaves washed individually with distilled water. Using a dissecting microscope, invertebrates were picked from the material, identified to genus (with the exception of Chironomidae and Ceratopogonidae as well as non-Insecta taxa-all identified to family), and classified according to their functional feeding group using Merritt et al. (2008). These data were used to calculate abundance, richness, Margalef's richness, percentage of EPT (Ephemeroptera + Plecoptera + Trichoptera), percentage of Chironomidae, and the percentage of each functional feeding group. Residual leaf material was dried at 60°C for 48 h and ashed at 500°C for 2 h; AFDM was calculated by subtracting ash mass from dry mass. Percent AFDM lost was the difference between final and initial AFDM (initial was ❖ www.esajournals.org calculated to be 95.2% of the starting dry mass during a preliminary study), and percentage of lost AFDM per degree-day was calculated by dividing percentage lost AFDM by accumulated degree-days during the incubation period. Ash mass, that is, the leaf inorganic mass, was used as an indicator of the degree of very fine sediment entrainment in biofilms on leaf material.

Catchment and reach-level explanatory variables
Explanatory variables were classified into five groups under two spatial scales: catchment and reach variables. Catchment variables included those related to harvest, catchment forest condition (structure and composition), and landscape characteristics. Reach variables included those related to riparian forest condition and stream morphology. In addition to being used as indicator variables, abiotic endpoints (sediment loading, water chemistry, DOM quality, and water temperature) were treated as explanatory variables for a separate set of models in the analyses of biotic indicators (leaf decomposition, BMI community, and biofilm). See Appendix S1: Table S2 for the list of variables available within each category.
Harvest levels in catchments.-Harvest variables were calculated from the GIS information on stands harvested each year since 1986 (earliest year with harvest data available). The area of each catchment harvested each year by different harvesting methods was calculated and then summarized into metrics of cumulative area harvested in the last 5, 10, 20, and 30 yr. Harvesting method was either clearcut (usually 80% or more of the trees removed) or partial harvesting (usually 35-50% of the trees removed, including commercial thinning of spruce plantations and selection-based harvesting in deciduous stands). Total cumulative harvest area was the sum of clearcut and partial harvests over each of the four time periods. Because six of the 12 BB catchments underwent some harvesting between the sampling periods in 2014 and 2015, the values of harvesting variables differ between years.
Landscape features.-Landscape features that could potentially affect stream ecosystems were quantified. Catchment area upstream of the sampling sites was obtained in ArcGIS following watershed delineation from the 5-m DEM and using a flow accumulation threshold of 0.1 ha.
Stream length and density (unit length/watershed area) were calculated from stream shapefiles. Road variables were calculated from the road shapefiles from BB. Road crossing frequency was the number of road crossings divided by total stream length, and road density was the sum of all road lengths divided by the area of the study catchment. Note that roads in BB are gravel roads and that salt or dust settlers are not used. Mean elevation and slope of the catchment were calculated from the DEM. Percent effective variable source area (% effVSA), which describes the hydrological connectivity between the catchment and the stream, was calculated on a hydrologically conditioned DEM, where the D8 algorithm was used to calculate contributing areas and the stream network (a contributing area of 0.25 ha was used as the channelization threshold) following the method used by Mengistu et al. (2014).
Catchment forest condition.-Forest composition (species assemblages) and structure (age, height) variables were derived from FRI data. A percent of total catchment area covered by each tree species was calculated for the overstory (>2 m height) and understory (<2 m) layers, and forest stands were classified as deciduous, coniferous, or mixed forest depending on which overstory and understory tree species predominated. The average height, size class, crown closure, and developmental stage of the forest in each catchment were also calculated from the FRI. In BB, six additional forest structure metrics were calculated from LiDAR for each catchment by averaging the structural metric values based on 6 pulses/m 2 from each 5 9 5 m cell in each catchment: P90, CrC2, CrC10, S2, S10, and VCI (see Appendix S1: Table S2 for details).
Riparian forest condition. -From the FRI data, percentage of deciduous/coniferous/mixed forest and riparian height were obtained for a 30 m wide band along each side of the stream adjacent to and extending upstream of the sampling reach. Percent canopy openness was measured by using photographs taken every 10 m within the 80-m sampling reach with a Canon EOS 50D camera with a 185°SuperFisheye (5.6 mm F/5.6) lens and a self-leveling mount (R egent Instruments, Ville de Quebec, Quebec, Canada). Photographs were processed using WinSCANOPY 2009a for Canopy Analysis (R egent Instruments).
Stream morphology.-Stream morphology was characterized by visually surveying eight 10 m long sections and averaging for each 80-m reach for flow structure (% riffle/run/pool), substrate composition (% bedrock/boulder/cobble/gravel/ sand/silt/clay), width and averaged depth, and number of large woody debris (>10 cm diameter) structures.

Statistical analysis
Box plots were constructed for each stream indicator to look for differences among streams and between years, followed by a linear mixed model ANOVA testing for between-year differences (year as a fixed factor and stream as random). Log-transformation was performed when necessary to meet the assumptions of normality and homogeneity of variance. The inter-annual congruence of each stream indicator was further assessed by conducting a Pearson's productmoment correlation analysis between 2014 and 2015 results.
For model construction and selection purposes, only BB streams were included because (1) several variables derived from high-quality forest and landscape characterization data were only available for BB (Appendix S1: Table S2), and these data may better separate the effects of harvesting from those of forest condition or landscape characteristics; and (2) the geological and topographical differences between BB and MC resulting from these areas being part of different ecoregions could confound the comparison between reference and harvested sites. Therefore, MC sites were only used to determine a regional range of natural variation in the stream indicators and to compare indicator values in BB to those in nearby unmanaged systems using box plots.
Linear regression models were built with stream indicators as response variables and catchment and reach variables as explanatory variables (EVs; Table 1). As a first step during EV selection for the regression analyses, a Pearson's correlation analysis was conducted among EVs within a category. When several variables within a category were correlated (as for harvest, catchment forest condition, and stream morphology categories), PCAs were performed for each category with all the variables in it (on centered and scaled data), and the PCs that captured the main gradients within the category were selected for regression analyses. For the harvest category, two separate PCAs were done with 2014 and 2015 data because harvesting occurred in some of the catchments between years. The variance explained by each PC as well as the variables correlated with the selected PCs is summarized in Table 1; PCA results for each EV category can be found in Appendix S2: Figs. S1-S4. Landscape variables did not show strong correlations among each other and the PCA did not identify strong gradients, so the two variables that were most related to the research questions and that can potentially influence stream condition were selected, that is, road density (which was correlated with number of stream crossings) and percentage of effVSA. There were only five variables within the riparian forest condition category, so instead of conducting a PCA, the two most representative variables were selected: percentage of riparian conifer (negatively correlated with three other riparian variables: percentage of riparian mixed, percentage of deciduous forest, and height) and canopy openness. Stream morphology was not included in the regression analyses for water chemistry, DOM quality, and water temperature under the assumption that this category of EVs would not influence the response variables. All the variables used in the regression analyses are shown in Table 1.
A second set of regression models was built to assess how abiotic endpoints were associated with biotic indicators and included the following: leaf inorganic mass (representing fine inorganic sediments; see Sediment deposition), water chemistry PC1 (representing cation/carbon concentrations; see Results: Water chemistry and DOM quality), water chemistry PC2 (representing DOC and metal concentrations; see Results: Water chemistry and DOM quality), DOM humification (HIX) and aromaticity (SUVA), and maximum August temperature (Table 1).
Model selection was conducted using ordinary least-squares regression and Akaike Information Criterion adjusted for small sample sizes (AIC c ; MuMIn package in R version 1.15.6, Bart on 2016). All potential models were first constructed from the subset of nine (or eight) catchment and reach EVs (including the null model), and the three best models (lowest AIC c score) were selected for each indicator (Anderson 2008); the same was done with the subset of six abiotic ❖ www.esajournals.org stream EVs for the biotic indicators. DAIC c , Akaike weights (x i ), and R 2 were calculated to evaluate the relative importance and the variance explained by each model.
Regression analyses were complimented with hierarchical partitioning analysis (hier.part package in R version 1.0-4; Walsh and Mac Nally 2013), which assessed the relative importance of each of the nine (or eight) EVs and separated the independent and joint contributions of each EV to the variability of the response variable (Chevan and Sutherland 1991). Hierarchical partitioning results were complimented with simple linear regressions, and the coefficient of determination values (R 2 ) and slope coefficients (only when regression models were statistically significant at a = 0.1) are reported. For the biotic indicators, analyses were conducted twice, first with the nine catchment and reach variables and then with the six abiotic indicators.
Leaf pack BMI community composition dissimilarity across streams was visualized in two Notes: DOC, dissolved organic carbon; DOM, dissolved organic matter; VSA, variable source area; LWD, large woody debris. For principal components (PCs), the variance explained by each PC and the variables positively (+) and negatively (À) correlated with each PC have been specified.
† Abiotic stream indicators were treated as explanatory variables only in the case of regression models explaining biotic indicators.
‡ Stream morphology was not included as an explanatory variable in the regression models explaining water chemistry, DOM quality, and water temperature.
§ The sign (+/À) of the loading of the variables in the principal component has been inverted for consistency and readability.
dimensions using non-metric multidimensional scaling (NMDS) graphs. To assess how the ordination of these streams based on their BMI communities related to environmental gradients, EVs were fitted to NMDS ordinations as linear vectors using the envfit function in R (based on 999 random permutations). All these tests were performed using vegan (Oksanen et al. 2017) in R (version 3.3.2., R Core Team 2016).

Sediment deposition
The deposition of fine inorganic sediments (range: 0.021-0.684 g) in traps and the entrainment in leaf pack biofilms (range: 0.11-0.59 g) were significantly lower in 2015 than in 2014 (F 1, 174 = 154.3 and F 1, 158 = 38.8, P < 0.001), whereas the organic percentage in traps (range: 20-64%) was lower in 2014 than 2015 (F 1, 177 = 19.8, P < 0.001; Fig. 2). Yet, inter-annual correlations were high (r = 0.60, 0.87, and 0.91, respectively). Reference streams, overall, had less fine inorganic sediment deposition and entrainment than most BB streams, but a greater proportion of organic fine sediments. Both indicators of inorganic sediment load in the streams were highly correlated (r = 0.87-0.91), so only model selection results for leaf biofilm entrainment are presented.
According to AIC c model selection, models containing only road density were the best to explain inorganic sediment load in both years, and the association was positive (+; Table 2). However, models also including stream morphology PC1 (À, which is positively related to percentage of riffle, percentage of cobble), percentage of riparian coniferous (+), and/or percentage of effVSA (+) closely followed the best model.
Of all nine variables considered, road density (+) had the highest independent effect on fine inorganic sediments in both years (I = 27.2% and 36.6%) according to hierarchical partitioning, as well as the variable explaining the highest proportion of the variance according to linear regression (R 2 = 0.74 and 0.73), followed by percentage of riparian conifer (+) and stream morphology PC1 (À; Appendix S3: Fig. S1).

Water chemistry and DOM quality
Trends in water chemistry across streams (represented by PC1 and PC2) remained relatively constant between years (r = 0.99 and 0.97, respectively), and including or excluding reference streams in the PCA yielded very similar results (the PCAs that include reference streams are shown and described). PC1 (~55% of the variance captured) was strongly and positively related to pH, conductivity, alkalinity, Ca, and DIC and negatively to SiO 2 ; hence, PC1 represented cation and carbon dynamics (Fig. 3). PC2 (~21% of the variance captured) was positively related to DOC, Cl, Al, and Fe, representing DOC and associated metal dynamics. Reference streams had the lowest PC1 values and intermediate PC2 values.
Among all the linear regression models built with the eight EVs to explain the trends in cation/carbon dynamics in BB, the model including harvest PC2 (À, related to cumulative harvest), percentage of riparian conifer (+), and canopy openness (À) had the strongest support in 2014 according to AIC c , but it was closely followed by the model excluding canopy openness (Table 2). In 2015, the best supported model included forest condition PC1 (+, related to mature, deciduous forests), road density (+), and percentage of riparian conifer (+), followed by the model with percentage of riparian coniferous as the only EV. Percent riparian conifer (+) had the highest independent effect on cation/carbon concentrations (I = 38%) and explained the largest amount of variance (R 2 = 0.74) in 2014, followed by road density (+) and harvest PC2 (À) (Appendix S3: Fig. S2). In 2015, both percentage of riparian conifer and road density had equally high independent effects (I = 32.2% and 30.9%) and R 2 s (0.75 and 0.63, respectively).
For DOC/metals dynamics in BB, the model including road density (+) and harvest PC1 (À, positively related to partial harvest in both years and negatively to recent clearcut in 2014) had the strongest support in both years, but it was closely followed by the model with only harvest PC1 in 2014 (Table 2). Harvest PC1 had the highest independent effect on DOC/metals dynamics in both years (À, I = 34.9% and 28.6%) and explained the largest amount of variance (R 2 = 0.59 and 0.49), followed by road density (+) and canopy openness (+; Appendix S3: Fig. S2). The clearcut PC was positively related to DOC and metal concentrations, whereas the partial/total harvest PC was negatively related.
Among the optical metrics describing DOM quality, only metrics describing DOM humification (HIX) and aromaticity (SUVA and E2E3) showed temporal consistency within years, so only models for HIX and SUVA are presented. Dissolved organic matter humification was best explained by road density (+) in 2014, but it was closely followed by the model with both percentage of effVSA (+) and road density (Table 2). In 2015, the model including road density (+), clearcut PC2 (À), forest condition PC2 (+, related to percentage of mixed forest and S2), and percentage of riparian conifer (+) was the most supported, but the model with only road density had also considerable support. Road density had the highest independent effect on DOM humification in both years (I = 33.4% and 38.3%), as well as the highest R 2 (0.63 and 0.60), followed  Notes: DOC, dissolved organic carbon; DOM, dissolved organic matter; VSA, variable source area. The catchment and reach explanatory variables included in the three best regression models (M1, M2, and M3) are presented, as well as their corresponding DAIC c , weight (%), and R 2 (R 2 for simple and adjusted R 2 for multiple regression models). The +/À symbol indicates the sign of the coefficient for that variable within that regression model. In bold are the best models (lowest AIC c value). by percentage of effVSA (+) and percentage of riparian conifer (+; Appendix S3: Fig. S2).
Finally, DOM aromaticity was best explained by percentage of riparian conifer (À) or road density (À) in 2015, but the null model was the likeliest one in 2014. Harvest PC1 (positively related to partial harvest, negatively to recent clearcut) had the highest I in 2014 (34.8%), but the linear regression (LR) was not significant (Appendix S3: Fig. S2). In 2015, percentage of riparian conifer had the highest I and R 2 (À, 30.6% and 0.34, respectively), followed by road density (À).

Water temperature
Maximum August water temperatures (range: 8.9-15.2°C) were strongly correlated between years (r = 0.93) and were best explained by canopy openness (+) and percentage of riparian conifer (+) ( Table 2). Percent riparian conifer had the highest independent effect on temperatures in BB streams in both years (I = 26.3% and 20.6%, respectively), but these values were similar for canopy openness (+) and harvest PC1 (À, positively related to partial harvest and negatively to recent clearcut) in 2014, and for canopy openness and road density (+) in 2015 (Appendix S3: Fig. S3). However, percentage of riparian conifer explained less variance than canopy openness or road density in both years according to the linear regressions.

Biofilm and algal biomass
Total biofilm biomass values (range: 0.2-2.0 g/m 2 ) were significantly higher in 2015 than in 2014 (F 1,43 = 6.7, P = 0.01), but there were no significant differences in algal biomass (range: 0.3-24.9 mg/m 2 ; F 1,43 = 1.4, P = 0.24); both endpoints showed inter-annual congruence (r = 0.68 and 0.59, respectively; Fig. 4). Reference streams had total biofilm biomass and algal biomass values similar to those at the low end of the range for BB streams. Autotrophic index values (data not shown) followed similar trends as algal biomass, so only the latter was used in subsequent models.
Road density (+) was included in all but one of the six best models explaining biofilm biomass (Table 3). In addition, the best model in 2014 included percentage of effVSA (+) and canopy openness (+), and the best model in 2015 included percentage of riparian conifer (À). In 2015, road density was also the variable with the highest R 2 (0.52) and independent effect (I = 27.3%) on biofilm biomass, but percentage of effVSA had a stronger independent effect (although a lower R 2 ) in 2014 (Appendix S3: Fig. S4). The abiotic variable with the strongest independent effect on biofilm biomass was DOM humification (+) in 2014 and fine inorganic sediments (+) in 2015 (Table 4); both variables explained 30% of the variance in 2015 linear regressions, but DOM humification explained considerably more than sediments in 2014 (79% vs. 46%; Appendix S3: Fig. S6). The model including only DOM humification explained as much of the variance as the best model built based on the nine catchment and reach EVs in 2014, but this model explained considerably less variance in 2015 and was only slightly more likely than the null model.
In 2014, the best model explaining algal biomass included canopy openness (+) and harvest PC1 (+, positively related to partial harvest and negatively to clearcut), and it was twice as likely as the model with only canopy openness (Table 3). However, in 2015 the model including canopy openness was the most likely one, but it was only 10% more likely than the null model. Canopy openness (+) was the variable with the highest independent effect on algal biomass in both years (I = 51.4% and 39.7%), as well as the only variable producing a significant simple linear regression model (Appendix S3: Fig. S4). None of the models including abiotic variables were significant, and the null model had stronger support (Table 4).

Leaf decomposition and benthic macroinvertebrate communities in leaf packs
Leaf decomposition.-Average decomposition rates across sites (range: 0.06-0.15% AFDM per degree-day) were significantly greater in 2015 compared to 2014 (F 1, 157 = 4.3, P = 0.04), but 2014 and 2015 values were strongly correlated across streams (r = 0.70; Fig. 4). Two reference streams (MC3 and MC4) had decomposition values in the mid-range of those in BB, whereas MC2 values were similar to those at the low end of the range.
The best model explaining leaf decomposition in BB included harvest PC2 (+, cumulative harvest) and percentage of riparian conifer (À) in 2014, and it was twice as likely as the second-best model, which included stream morphology PC1 (+, representing percentage of riffle, percentage of cobble, and width) as the only EV (Table 3). In 2015, the best supported model included road density (À), canopy openness (+), and stream morphology PC1 (+), and it was twice as likely as the next best model excluding stream morphology PC1. Stream morphology PC1 (+) had the greatest independent effect on leaf decomposition in both years (+, I = 19.4% and 21.3%, respectively), but was closely followed by percentage of riparian conifer (À) and harvest PC2 (+) in 2014, and by canopy openness (+),percentage of riparian conifer (À), and road density (À) in 2015 (Appendix S3: Fig. S5). The most strongly supported regression model based on the six abiotic indicators explained as much of the variance as the most supported models based on the nine catchment and reach EVs, but it only included water chemistry PC1 (À, cations/carbon; Table 4). Water chemistry PC1 had the highest independent effect on leaf decomposition and R 2 but was closely followed by sediments (À; Appendix S3: Fig. S6).   Benthic macroinvertebrate community. -Benthic macroinvertebrate abundance (range: 36-480 individuals/leaf pack) and richness (range: 10-23 genera/leaf pack) were higher in 2015 than 2014 (F 1, 158 = 70.9 and 22.4, respectively, P < 0.001) and showed an inter-annual correlation of r = 0.30 and 0.58, respectively (Appendix S4: Fig. S1). The proportion of EPT (range: 22.3-69.7%), chironomids (range: 19.6-74.7%), and functional feeding groups (Appendix S4: Fig. S2) also varied across streams. In some streams, shredders were the most abundant group (e.g., BB01 and BB02-~70%, dominated by Leuctra and Capniids), whereas in other streams collectorgatherers (e.g., BB06 and BB08-~40-50%, Baetis and Ephemerella), scrapers (BB05-41%, Epeorus), or predators (BB04 in 2014-50%, Dicranota and Rhyacophila) dominated. In reference streams MC3 and MC4, metrics describing the BMI community in the leaf packs tended to fall within the range of BB stream values. MC2 had the lowest percentage of EPT and highest percentage of chironomids of all 15 streams; percentage of shredders at this site was closer to the upper range for BB sites.
Because abundance and richness were correlated (r = 0.70-0.78), only abundance modeling results are presented. In 2014, the null model had more support than the next best supported model explaining BMI abundance on leaf packs in BB (which included canopy openness (+) or effVSA (À) as the only EV; Table 3); in 2015, however, the model with canopy openness was the most likely (+), followed by those including harvest PC1 (+, mostly representing partial harvest) and road density (À) in addition to canopy openness (DAIC c >1.79 though). In 2014, percentage of effVSA had the highest independent effect on BMI abundance (I = 28.7%), followed by forest condition PC1 (related to mature, deciduous forests) and canopy openness, although none of the simple linear regressions including each of these three EVs were significant (Appendix S3: Fig. S5). In 2015, canopy openness (+, I = 52.7%) had a higher independent effect than percentage of effVSA (À). For the regression models based on the six abiotic indicators in BB, the best models included temperature (+) and DOM humification (À) or temperature and inorganic sediments (À; Table 4). Temperature was the variable with the highest independent effect on BMI abundance in both years, but in 2015, unlike in 2014, the linear regression was non-significant (Appendix S3: Fig. S6).
The proportion of shredder BMIs in BB leaf packs was best explained by percentage of effVSA (+) and percentage of riparian conifer (+) in 2014, but the model with only stream morphology PC1 (À, representing percentage of riffle, percentage of cobble, and width) had considerably more support than the remaining models in 2015 (Table 3). Percent effVSA (+) had the highest independent effect on the proportion of BB shredders in 2014 (I = 24.3%), followed by percentage of riparian conifer (+), stream morphology (À), forest condition PC2 (À, related to percentage of mixed forest and S2), and road density (+; Appendix S3: Fig. S4). In 2015, stream morphology had the highest I (36.5%), followed by percentage of riparian conifer, and simple linear regressions were only significant for these two variables. For the abiotic indicators, the 2014 model including DOM humification (+) and water chemistry PC2 (À, related to DOC and metals) had the greatest support; in 2015, the model including only DOM aromaticity was more likely (À; Table 4). In 2014, DOM humification (+) had the highest independent effect on percentage of shredder, but it was closely followed by inorganic sediments (+) and water cations/carbon (+) (simple linear regressions with each variable explained 60% of the variance); but, in 2015, DOM aromaticity (À) had the highest I and R 2 (Appendix S3: Fig. S6).
Reference streams showed a distinct BMI community composition in leaf packs in both years relative to harvested sites, with differences being driven mainly by the presence of Taenionema, Diplectrona, and Isoperla in reference streams (Fig. 5). Several EVs were correlated with the NMDS ordination of streams based on their leaf pack BMI community, but relationships were stronger with reach EVs than with catchment EVs for BB sites. Among reach variables, stream morphology PC1, inorganic sediments, DOM humification, and water cations/carbon were significantly correlated with the BMI ordination in both years. Road density and forest condition PC2 (related to overmature/young forests) were the catchment variables significantly correlated with the BB BMI ordinations. With respect to the direction of these correlations, most variables identified a similar gradient with the ordination; that is, the ordination of streams based on their BMI community was correlated with the ordination of streams from high to low inorganic sediments, water cations/carbon, DOM humification, percentage of run/silt, and road density (Fig. 5).

Sediment deposition
Roads and harvesting in forest management are typically the most important contributors to elevated sediment levels in streams (Reid and Dunne 1984, Croke and Hairsine 2006, Kara et al. 2014, Al-Chokhachy et al. 2016, Davies et al. 2016, and roads are almost always greater contributors than harvesting itself (Croke andHairsine 2006, Webster et al. 2015). Here, we show that this also occurs in an intensively managed landscape under best forest management practices. Density of BB roads (which are gravel) was the best predictor of the amount and inorganic/ organic ratio of fine inorganic sediments deposited in traps and entrained in leaf biofilms, whereas the degree of harvesting had a considerably weaker effect (only entrainment in leaf biofilms presented in this article). Sediment deposition increased at the road-intensive sites even in 2015 when overall fine sediment deposition was much lower. The lower deposition in that year was probably due to lower precipitation during trap deployment (105 vs. 67 mm in 2014 and 2015, respectively), confirming that sediment transport to streams is influenced by rainfall intensity (Croke and Hairsine 2006). Regardless of the source and hydrologic conditions, the current BMPs (in conjunction with legacy effects from past management such as an older road system) in BB, an intensive forest management setting, are insufficient to prevent sediment increases in streams.

Water chemistry and DOM quality
The first water chemistry PC (conductivity, alkalinity, pH, Ca, and DIC) was strongly and positively related to forest management intensity, mainly to road density. Although water cation/ carbon concentrations did not show a strong relationship with the harvest PCs, they showed a strong individual correlation (r = 0.64-0.79) with recent harvest (<10 yr). This agrees with previous studies that have shown increased conductivity and cations being delivered to streams from roads, as well as increased levels of cations after logging due to a reduction in plant uptake, increased decomposition and mineralization, elevated water fluxes, etc. (Danehy et al. 2007, Kreutzweiser et al. 2008b, Richardson and B eraud 2014. It is interesting to note that positive correlations between harvest and cations/ carbon were only observed for recent harvest (<10 yr), whereas the accumulated harvesting over longer periods of time (<30 yr) had a negative effect, possibly due to regenerating and rapidly growing vegetation absorbing these dissolved compounds. Fine inorganic sediments and water chemistry PC1 were strongly correlated, indicating higher delivery of both waterborne fine and dissolved materials in more intensively managed catchments. As with sediments, percentage of coniferous forest (especially in the riparian forest) was strongly and positively related to water chemistry PC1, which suggests either that these forests were less effective at intercepting waterborne materials or that the greater management intensity in catchments with coniferous riparian forests (usually catchments with conifer plantations) delivered higher rates of waterborne materials.
Dissolved organic carbon, nitrogen, Fe, Al, and Cl in BB stream waters tended to increase with recent clearcut (<5, <10 yr) but decreased with partial harvesting, suggesting that the complete removal of trees had a greater effect on stream water chemistry. Increases in nutrients and metals in streams after intensive logging and the amelioration of nutrient exports to some extent by partial or stem-only harvesting have been observed previously (Feller 2005, Kreutzweiser et al. 2008b. Fe and Al concentrations in BB streams were also strongly related to roads (and sediments), showing that these metals were probably delivered to streams attached to sediment particles from roads. Although elevated at some sites, all the water chemistry parameters were below levels deemed unsafe for aquatic life and drinking water (CCME 1999).
Of the DOM quality variables in the current study, DOM humification (HIX) was the most temporally consistent and strongly related to road density. Streams with more humic DOM also had the greatest DOC, fine sediments, DIC, and Ca, suggesting that sites with high road density generally delivered more waterborne materials such as terrestrial DOM characterized by more humic compounds (Kalbitz et al. 1999). That percentage of effVSA was also positively related to HIX suggests that both anthropogenically (via roads) and naturally (via effVSAs) enhanced hydrological connectivity contributed to the delivery of more humic terrestrial compounds. In 2015, DOM aromaticity (SUVA) was negatively related to HIX and road density, which is surprising considering that terrestrial DOM is characterized by more aromatic DOM and higher SUVA values (Weishaar et al. 2003), but others have also reported an apparent discrepancy between HIX and SUVA (Catal an et al. 2017). Because humic acids contain a greater proportion of aromatic groups than fulvic acids (Findlay and Sinsabaugh 2003), it could be hypothesized that the observed increase in humic substances is mainly driven by the more labile fulvic acids. Regardless, neither were related to harvesting, and this agrees with other studies (Burrows et al. 2013, Cawley et al. 2014, De Wit et al. 2014. Because DOM quality controls its reactivity, influences ecosystem functions, and affects microbial food webs (Weishaar et al. 2003, Docherty et al. 2006, Emilson et al. 2017, these forest management-related changes could have biological implications in streams, such as altered biofilm production as discussed in Discussion: Biofilm and algal biomass.

Temperature
While riparian forest removal typically leads to increased and more variable water temperatures, especially in the summer , the effects of upland harvesting on stream temperatures are less predictable, with some suggesting that riparian buffers are effective at reducing such impacts (Wilkerson et al. 2006, Clinton 2011 and others not (Kiffney et al. 2003, Witt et al. 2016). In the current study, even though all the streams had a no-harvest riparian buffer (usually about 30 m wide), percentage of recent clearcut (<5 and <10 yr) was strongly and positively related (especially in 2014) to maximum summer temperature. The warming of groundwater in areas of upland clearing can increase stream water temperature , and although we did not monitor groundwater temperature, our results support this hypothesis. This is further supported by the fact that only recent clearcut, and not partial harvesting, contributed to warmer stream temperatures. Our results also agree with the 5-to 10-yr postlogging recovery time for stream temperature proposed by , since we only detected the effects of clearcut in the most recent 5-and 10-yr harvest periods. Road density was positively related to water temperature in BB as well, which could result from the bare roads promoting the warming of shallow groundwater, or runoff from warmed road surfaces during rain events.
Maximum temperatures in the BB streams never exceeded 16°C and were thus well below the threshold of 24°C set to protect coldwater fish species such as brook trout (EPA 1986). However, the observed temperature differences among streams apparently affected other organisms in this study, with temperature being positively related to BMI abundance and richness in leaf packs.

Biofilm and algal biomass
Total biofilm biomass in these streams responded positively to management intensity, agreeing with previous observations of elevated biofilm biomass at logged sites even when 30-m buffers were retained (Kiffney et al. 2003), as well as in thinned forests (Danehy et al. 2007). In the current study, roads had the strongest positive effect on biofilm biomass, followed by recent clearcut, and both may have contributed to the strong positive relationship between biofilm biomass and DOM humification and DOC/metal concentrations. The hypothetical DOC-related increase in biofilm at intensively managed sites would have occurred despite its more humic nature and may be because microbial communities in forested headwater streams are adapted to the inputs of terrestrial-and more humic-DOC (Kreutzweiser and Capell 2003, Burrows et al. 2013, Emilson et al. 2017. If our hypothesis is true and the increased DOM humification is mainly driven by the less aromatic/more labile fulvic acids, it would explain why biofilms responded positively to increases in humification. If not, these results would add evidence in favor of DOM quantity overriding the importance of quality in terms of stimulating ecosystem productivity (Marcarelli et al. 2011, Roiha et al. 2012.
Even though sediments are known to negatively impact biofilms (Izagirre et al. 2009, Jones et al. 2012, we saw increases in total biofilm biomass at sites with higher sediment deposition. This was probably driven by the greater inputs of DOM supporting heterotrophic microbes (as well as entrained FPOM), since algal biomass in BB was not related to forest management intensity in this study and others (Davies et al. 2016) and detritivores (bacteria and fungi) tend to competitively exclude primary producers under low light and nutrient availability and high allochthonous C supply (Daufresne and Loreau 2001, Mindl et al. 2005, Danger et al. 2007. Algal biomass was mainly related to canopy openness, indicating that primary production may be light, but not nutrient, limited in the BB streams. BB07 and BB10, the streams with the highest algal biomass values, had reaches of narrow buffer zones and sections that were affected by recent windthrow (personal observation), which resulted in more open canopies. Increased riparian canopy openness from blowdown increased light levels and algal biomass in streams, creating effects similar to riparian harvesting (Lecerf et al. 2012) and decreased buffer width (Wilkerson et al. 2010).

Leaf decomposition and benthic macroinvertebrate communities
Leaf decomposition. -Leaf decomposition was negatively affected by management intensity (namely roads) in BB despite the associated warmer temperatures, with decomposition being approximately 15% lower in more intensively managed streams. This could be the result of the enhanced delivery of waterborne particulate and dissolved materials associated with roads, as evidenced by the strong negative effect of water cation/carbon concentrations and fine inorganic sediments on decomposition. Inundation of leaves by sediments lowers their decomposition rates Waide 1982, Lecerf andRichardson 2010), and although the leaf packs in this study did not show signs of burial in sediments (personal observation), fine sediments were entrained in leaf biofilms, and that may have resulted in their reduced palatability and decomposition (Danger et al. 2012). Interestingly, the negative relationship between conductivity/ Ca values and decomposition observed in BB contrasts with the positive relationship observed elsewhere, since Ca tends to be a limiting nutrient for decomposing microbes (Egglishaw 1968, Lecerf andRichardson 2010). This could be a result of lower microbial activity (Emilson et al. 2017), lower decomposition attributable to invertebrates (McKie et al. 2006), and/or adverse effects of sediments overriding the benefits of higher Ca levels (as well as that of warmer water temperatures) to microbial communities.
Although the effect of harvesting on leaf decomposition was weaker than that of roads in the current study, here and elsewhere decomposition can be negatively affected by harvesting when riparian buffer zones are retained (Kreutzweiser et al. 2008a, Lecerf andRichardson 2010). Yeung et al. (2017) proposed that the recovery time of this crucial ecosystem function ranges between 8 and 15 yr post-harvest; BB streams fall within the range, since negative correlations between harvest and decomposition were only significant for harvesting in the past 5 and 10 yr. In fact, the accumulated harvesting over longer periods of time (<30 yr) tended to have a positive effect on leaf decomposition.
Benthic macroinvertebrate community.-Although both abiotic indicators and stream function were affected at sites with higher forest management, the BMI communities in leaf packs did not indicate biological impairment. In fact, the most abiotically impacted streams (based on sediments and water chemistry) had among the highest percentage of EPT values, which contrasts with other studies showing decreased richness and/or abundance in leaf pack BMIs (Kreutzweiser et al. 2008a, Lecerf andRichardson 2010). Higher abundance and richness of BMIs in BB were positively related to canopy openness and negatively to effVSAs, which may be due to the higher nutritional quality of algae produced at warmer, more exposed sites in comparison with the terrestrial materials delivered through effVSAs (Guo et al. 2016). While these indicators were not directly related to management intensity, BMI composition was related to physicochemical characteristics (sediments, DOM quality, water chemistry, and temperature) affected by management activities (Kreutzweiser et al. 2008b).
Relative shredder abundance on leaf packs was as low as 8% in some streams, a surprising result considering that these are forested headwater streams (Vanote et al. 1980), but it supports the idea that leaves are important as both a food source and habitat to non-shredder BMIs (Richardson 1992). Increased percentage of shredders has been attributed to increased detrital resources after harvesting (Jackson et al. 2007, Medhurst et al. 2010, but trends in the opposite direction have also been reported (Smith et al. 2009, Richardson andB eraud 2014). In BB, it seems that shredders had a competitive advantage over the other functional feeding groups in streams with highly managed catchments and the associated higher inputs of terrestrial materials (e.g., 70% shredders in BB01 and BB02), with all the feeding groups but shredders declining in relative abundance as sediment deposition increased. In addition, the increase in percentage of shredders with increasing percentage of effVSA also suggests that shredders respond positively to the delivery of terrestrial materials (via the increased catchment-stream hydrological connectivity in the case of effVSAs).
Increases in shredder abundance, however, did not translate into higher decomposition of leaves; in fact, relative and absolute shredder abundance were negatively correlated with decomposition in 2014 (r = À0.64 and À0.61) and relative abundance in 2015 (r = À0.37). Riipinen et al. (2009) proposed that leaf decomposition is mainly related to the size and biomass of shredders, rather than shredder abundance or richness. In our study, the increase in shredder proportion in the most impacted streams was driven by small-bodied stoneflies (Leuctridae, Nemouridae, and Capniidae), and these families are known for shredding less than the larger caddisflies (Dangles and Guerold 2001) that were found at the less impacted BB streams. In fact, decomposition positively correlated with the absolute abundance of >1-mm Lepidostoma caddisfly shredders in 2015 (r = 0.49), but not with <1-mm Capniids and Nemourids, and negatively with Leuctra (r = À0.55), supporting this hypothesis.

Harvested (BB) vs. reference (MC)
Reference streams tended to have lower deposition of inorganic sediments and a higher organic fraction in sediment traps and different water chemistry (higher SiO 2 and lower Ca, Mg, conductivity, alkalinity, DIC, SO 4 , and nitrogen levels) from the BB sites. These differences are likely attributable, at least in part, to the absence of soil disturbance (harvesting) and roads in the unmanaged MC catchments. However, differences in water chemistry could also be attributed to the different geology underlying these areas as felsic volcanic rocks (characteristic of MC) are known to be rich in elements such as Si, Na, and K, which were higher at these sites. Sedimentary rocks (characteristic of BB) tend to be richer in dolomites and calcium carbonates, explaining the higher pH, Mg, Ca, and C levels in these streams (Feller 2005). Regardless of the relative influences of forest management or geology, these abiotic differences did not translate into significant biotic differences between harvested and reference streams. Reference streams had values for leaf decomposition, BMI abundance and richness, and total biofilm biomass that fell within the range observed for harvested sites.
Among reference streams, MC2 tended to have the most extreme values for most of the indicators, which may have reflected influences of a beaver pond located about 800 m above the sampling reach (Naiman et al. 1986). MC2 had higher water temperatures, sediment deposition, water DOC/Fe/Mg, DOM alkalinity/humification, and percentage of shredder, and lower water SiO 2 , decomposition rates, and BMI richness and abundance than the other two reference streams. Therefore, MC2 shared many of the same characteristics as the most intensively managed catchments in BB. Although this study was not designed to assess the effects of beaver dams on stream ecosystems, these results indicate that natural (beaver dam) and anthropogenic (forest management) disturbances could be shifting stream ecosystems in similar directions and via similar mechanisms. A common pathway of effects could potentially be through the release of terrestrial materials into the water. This release would be facilitated by the inundation of soils and flooding of vegetation in the case of beaver dams (Catal an et al. 2017) and by soil disturbance and road networks in the case of forest management.

CONCLUSIONS AND MANAGEMENT IMPLICATIONS
The suite of indicators measured in this study reflected the gradient in forest management intensity, especially road density, in BB. As a general trend, streams in the most intensively managed catchments had higher fine sediment deposition and entrainment, conductivity, cations, DIC, DOM humification, warmer water, more biofilm biomass, lower leaf decomposition, and higher shredder and lower grazer densities (Fig. 6). The main challenge was to discern between causality and correlation among variables, since harvesting intensity and roads were correlated with many other catchment and reach variables. Catchments in lower elevations tended to have more intensive forest management (i.e., degree of harvesting and road density), more coniferous-dominated riparian and catchment forests, and lower riffle/run ratios (Fig. 6). However, differences in stream indicators appeared to be related more to forest management activity (especially road density) than to forest composition, with road density affecting the greatest number of stream indicators in BB (Fig. 7). Overall, biotic indicators were more directly influenced by reach than catchment variables; reach variables that were both related (e.g., sediments and water chemistry) and unrelated (e.g., canopy openness) to management intensity were strongly associated with biotic indicators. In addition, many indicator responses (e.g., increased fine sediment deposition and conductivity) can be linked to forest management impacts from previous literature thereby inferring causality (Croke and Hairsine 2006, Kreutzweiser et al. 2008b, Richardson and B eraud 2014. Abiotic indicators, namely fine inorganic sediments and water chemistry, were strongly correlated with forest management intensity. Sediments were more strongly related to roads, whereas water and DOM quality were related to harvesting intensity and forest composition in addition to roads. Therefore, it appears that the higher degree of recent harvesting in low-elevation catchments (mainly conifer plantations) resulted in higher levels of ground disturbance overall, enhancing the biogeochemical processing in forest soils through temperature (from canopy opening) and organic matter (from slash) increases and promoting higher transport of water and waterborne materials to streams (Kreutzweiser et al. 2008b). This delivery was facilitated by the higher road density and stream crossings, which potentially contributed to greater runoff volumes and peak flows (Buttle et al. 2018) and resulted in more fine sediments, cations, and DOC reaching the streams in the most intensively managed catchments.
The measurable differences in habitat quality among streams related to intensity of forest management did not translate into impaired BMI communities. However, at the more sedimentinfluenced sites, there tended to be reduced leaf decomposition, which is a critical invertebrate-mediated ecological function in forested headwater streams (Gessner and Chauvet 2002). This underscores the importance of looking at both structural and functional biotic indicators to assess the effects of catchment disturbance on stream ecosystems (Christensen andBartuska 1996, Giller 2005), and the importance of measuring a suite of abiotic indicators to better understand the biological responses and effects pathways. For example, contrasting conclusions could be drawn from this study when looking at indicators of management intensity individually, ranging from site impairment (based on leaf decomposition), to no effect (based on leaf pack BMI abundance), to enhanced productivity (based on biofilm biomass). Evidence from our multi-indicator approach elucidated a potential effects pathway of higher inorganic sediment content in biofilms of organic matter potentially limiting or altering its use by microbial and BMI communities and resulting in reduced leaf decomposition rates. This was only possible by including a suite of indicators representing different components of the ecosystem and lends Fig. 6. Diagram summarizing the relationships among catchment and reach explanatory variables (gray and blue boxes, respectively), and abiotic and biotic stream indicators (yellow and green circles, respectively). Arrows between box/circles indicate hypothesized relationships based on regression analyses, and the up/down arrows after variables indicate positive/negative relations.
itself to an improved understanding of the linkages between disturbance and effects.
Biotic responses to forest harvesting can be positive or negative or show no changes (Richardson and B eraud 2014). In the current study, these indicators were related not only to area harvested, but also to roads, forest condition, and stream morphology, which could explain the breadth of responses in the literature as sites would vary in these variables. In addition, we demonstrated that assessing a gradient in forest management intensity rather than a discrete category (presence/absence of harvesting) or categories improved our understanding of the directionality and thresholds of the change in stream indicators. For example, leaf decomposition in BB streams was higher, lower, or similar to values in reference streams, and captured all three types of responses previously reported (Kreutzweiser et al. 2008a, Mckie and Malmqvist 2009, Lecerf and Richardson 2010. Because we considered a broad gradient in management intensity, we could show that decomposition responded positively to lower intensity and negatively to higher intensity in comparison with reference streams. In summary, results indicate that even in an intensively managed forest landscape such as BB, contemporary management practices (which include legacy effects from past management such as the road system) had little impact on Fig. 7. Strength of the effect of nine catchment and reach explanatory variables and six abiotic stream explanatory variables (columns) on 12 stream indicators measured in this study (rows) in 12 harvested Black Brook streams (BB). Dark gray indicates that an explanatory variable showed the strongest independent effect on an indicator; light gray indicates that a variable was significantly related (based on simple linear regressions, a = 0.1) to an indicator, but that the effect was less strong than that of the variable in dark gray. 14 or 15 indicates that the regression was statistically significant or that the effect was the strongest only in 2014 or 2015, respectively. forest stream biotic communities. Most of the harvesting-related effects were strongest for recent harvesting (<10 yr), suggesting that stream ecosystems in BB have recovered from the impacts of older harvesting. However, the changes in physicochemical habitat quality, especially the fine sediment deposition resulting from roads, and the reduced rates of leaf litter decomposition at the most intensively managed streams warrant attention. Increased protection of streams and their ecosystem services should be directed at reducing ground disturbance and minimizing roads and sediment inputs in catchments with high degrees of management activity.