Aquatic bryophytes play a key role in sediment-stressed boreal headwater streams

Forestry-related land use can cause increasing instream sedimentation, burying and eradicating stream bryophytes, with severe ecological consequences. However, there is limited understanding of the relative roles and overall importance of the two frequently co-occurring stressors, increased fine sediments and loss of bryophytes, to stream biodiversity and ecosystem functions. By using random forest modeling and partial dependence functions, we studied the relative importance of stream bryophytes and fine sediments to multiple biological endpoints (leaf-decaying fungi, diatom, bryophyte, and benthic macroinvertebrate communities; leaf decomposition) using field survey data from headwater streams. Stream bryophyte abundance and richness were negatively related to fine sediment cover, highlighting the detrimental effect of sedimentation on bryophytes. However, bryophyte abundance was consistently more important a determinant of variation in community composition than was fine sediment cover. Leaf decomposition was influenced by shredder abundance, water temperature and, to a lesser degree, stream size. Our results suggest that the loss of stream bryophytes due to increasing sedimentation, rather than fine sediments per se, seems to be the key factor affecting multiple biological responses. Enhancing the re-establishment of bryophyte stands could partly compensate for the negative impacts of sedimentation on bryophytes and, consequently, on several other components of boreal stream ecosystems.


Introduction
In many stream ecosystems, aquatic macrophytes modify erosion and sedimentation regimes and create habitat for other organisms (Sand-Jensen, 1998;Sand-Jensen & Pedersen, 1999). Macrophytes generally increase habitat heterogeneity and thus support diverse and abundant invertebrate communities (e.g., Taniguchi et al., 2003;Suurkuukka et al., 2014;Turunen et al., 2018). Macrophytes also provide refuge from predation and harsh environmental conditions such as floods or droughts (Rantala et al., 2004;Huttunen et al., 2017). In boreal streams, bryophytes and particularly mosses often dominate macrophyte communities. Bryophytes retain fine and coarse organic matter (Muotka & Laasonen, 2002;Turunen et al., 2018) and thus indirectly fuel detritus-based stream food webs. While it is known that bryophytes provide important habitat for stream invertebrates (Suren, 1993) and buffer invertebrate communities against variability (Huttunen et al., 2017) their importance to other ecosystem components, such as microbial communities and stream ecosystem functions remains poorly understood.
Forestry practices are a major cause for increasing sedimentation of boreal headwater streams (Marttila & Kløve, 2010;Turunen et al., 2017). Particularly drainage ditching, construction of unpaved forest roads, and deforestation have increased soil erosion and altered hydrological regimes (Vuori & Joensuu, 1996;Kreutzweiser & Capell, 2001;Sutherland et al., 2002). In Finland, 55% of peatlands have been drained by ditching to enhance forest growth (Turunen, 2008), increasing fine sediment loads, dissolved organic carbon, and nutrient concentrations in streams (Vuori & Joensuu, 1996;Vuori et al., 1998;Nieminen et al., 2017). The excess of fine sediments clog interstitial spaces and impair hyporheic flow exchange within the stream bed, with adverse ecological consequences on, for example, crevice-dwelling invertebrates and gravel-spawning fishes, such as salmonids (Jones et al., 2012;Sear et al., 2016). The instability and scouring of sediments may hinder the establishment of periphyton and interfere with ecosystem functions (McKie & Malmqvist, 2009;Louhi et al., 2017;Turunen et al., 2018). A particularly harmful consequence of excessive sedimentation is that it reduces bryophyte abundance in streams (Matthaei et al., 2006;Turunen et al., 2017). Indeed, experimental evidence suggests that bryophyte loss may have a stronger effect than sedimentation on stream invertebrates and ecosystem functions (Turunen et al., 2018), but the relative importance of these two stressors in field conditions remains unexplored.
We explored the relationship between aquatic bryophytes and fine sediment cover, and their relative importance to biological communities (leaf-decaying fungi, diatoms, macroinvertebrates) and leaf breakdown rate in boreal headwater streams. Our study streams spanned a gradient from near-pristine to heavily forestry-impacted streams where benthic habitats were almost completely covered by fine sediments. Headwater stream ecosystems are typically strongly influenced by natural variation in geology, topography and groundwater inflow (Winter, 2007;Annala et al., 2014). Therefore, we accounted for the effects of natural variability (e.g., water temperature, substratum structure and water chemistry) on benthic communities and ecosystem processes by using random forest modeling and partial dependency functions to enable the exploration of individual effects that bryophytes and sediments may have.
We tested three key hypotheses about the role of, and interactions between, bryophytes and fine sediments in boreal stream ecosystems. First, we expected that bryophyte cover and richness are negatively impacted by fine sediments. We further expected that bryophyte cover and fine sediments both contribute importantly to stream community composition, with bryophytes having a generally stronger effect of the two (see Turunen et al., 2018), especially for diatoms and macroinvertebrates. Fungal community composition was expected to be more driven by water chemistry than by bryophyte or fine sediment cover, since water chemistry typically controls fungal assemblages more strongly than substrate (Niyogi et al., 2003;Tolkkinen et al., 2015). Finally, bryophyte loss was expected to increase leaf decomposition rate, because the lower detrital resources in streams with low bryophyte cover can cause higher shredder aggregation in leaf bags and thus increase the decomposition rate (Turunen et al., 2018).

Study area
Our study area is in the region of Northern Ostrobothnia in central Finland in the headwaters of the River Iijoki basin (catchment area 14,191 km 2 ) (65°17-65°46 N and 27°35-28°35 E) (Online resource 1). The area is sparsely populated and characterized by extensive boreal coniferous forests (82% of the land area). Forestry, with strongly modified drainage network, poses the main anthropogenic pressure to stream ecosystems in the area whereas other pressures, such as agriculture and urban land use, are negligible (Online resource 1). Streams in the area are naturally circumneutral and typically colored by dissolved organic carbon from mires and peatland forests.
About 30% of the total land area of Finland is covered by peatlands, 55% of which has been drained (Turunen 2008). The peak era of draining activity occurred between the 1960s and 1990s (Paavilainen & Päivänen, 1995). Ditch networks were usually drained directly into stream channels that often were also straightened to enhance water removal. This has caused extensive soil erosion and flushing of sediments to streams (Marttila & Kløve, 2010). Construction of new drainage ditches was forbidden by the Finnish Forest Act in 1997 but existing ditch networks are regularly maintained to ensure drainage effectiveness.

Site selection
We selected 27 streams that differ in the intensity of forestry land use (mainly drainage ditching) and fine sediment input. We selected the sites based on map information and discussions with local forestry managers, and verified the suitability of candidate sites by field visits. All the streams were selected to be approximately of same size and have similar slope. The sites have a variable upstream catchment drainage intensity from 0 to 65% (mean 16%) of catchment area drained with ditches. In each stream we selected a 50-100 m long swiftly-flowing riffle habitat where all sampling (see below) was conducted. At impacted sites, the sampling reach was positioned immediately downstream of the main drainage network and, consequently, major sediment sources.

Habitat measurements
Field surveys were conducted between August and October 2013. Substrate structure was visually estimated in fifteen evenly distributed 0.25 m 2 plots as the percentage cover of different substrate size classes. Substrate size was classified by using modified Wentworth scale : fine organic matter = 0, sand (diameter 0.25-2 mm) = 1, fine gravel (2-6 mm) = 2, coarse gravel (6-16 mm) = 3, small pebble (16-32 mm) = 4, large pebble (32-64 mm) = 5, small cobble (64-128 mm) = 6, large cobble (128-256 mm) = 7, boulder (256-400 mm) = 8, large boulder and bed rock ([ 400 mm) = 9. Fine sediment cover for a reach was calculated as the mean proportion of class 0 and 1 across the ten plots. Bryophytes were sampled in fifteen 0.25 m 2 quadrats across a reach. For each quadrate, we recorded each species and their percentage coverage and used mean cover per reach in subsequent data analyses.
Water temperature was measured for 21 days (i.e., during leaf decomposition assays) at 1-h intervals with temperature loggers (iButton; Thermochron, Maxim Integrated, San Jose, USA). Mean water temperatures were calculated for each reach and used in the analyses.
The volume of large woody debris (LWD) was quantified by measuring the length and average diameter of all wood pieces larger than 5 cm in diameter that occurred within the sampling reach. In the case of tree trunks partly outside the channel, only the portion within the bankfull channel was measured. The total volume was divided by the area of the sampling reach to obtain wood volume per square meter. Shading by riparian vegetation was estimated in the middle of the channel at 20 evenly spaced points along the sampling reach. At each point, an estimate of the percentage cover of riparian canopy was obtained by pointing a cylindrical tube (diameter 15 cm) directly upwards and estimating the cover of tree canopy within view. Mean cover was then calculated for the whole sampling reach.
Total phosphorus, water color and suspended solids were measured from water samples collected simultaneously with the biological sampling using standard methods in the FINAS accredited laboratory of Finnish Environment Institute. Water pH was measured in situ with YSI Professional Plus -meter (YSI Inc., Yellowsprings, Ohio, USA).

Biological sampling
Periphytic diatoms were sampled from five randomly selected stones (diameter 10-20 cm) at each reach. Approximately 50 cm 2 of the upper surface of each stone was brushed with tooth brush and the dislodged material was preserved in ethanol. Four hundred diatom valves were counted and identified to species level for each sample, and the relative abundance of species was used in data analyses.
Macroinvertebrates were sampled by taking four 30-s kick samples that covered about 1.2 m 2 of the riffle habitat in each reach. The samples were pooled and preserved in 70% ethanol and all specimens (excluding Chironomidae, Simuliidae and Oligochaeta) were identified to species or genus level.
Leaf decomposition was measured by incubating dried alder leaves (Alnus incana L.) in each reach. We collected alder leaves prior to abscission and air-dried them at ? 23°C for 3 weeks. Dried leaves were enclosed in five bags (8-mm mesh size, six grams in each bag) which were anchored to house bricks and placed in depositional stream areas that naturally accumulate leaf litter. After 21 days of incubation, the leaf bags were collected and transferred to the laboratory, where the remaining leaf material was dried at 60°C for 24 h and weighed. Macroinvertebrates in the bags were collected and the abundance and richness of leaf shredding invertebrates were recorded. The samples were then ashed for 4 h at 550°C to convert dry mass to ash-free dry mass. Leaf mass loss rate (k) was calculated from the exponential decay model M f = M 0 * e -kt , where M f is the final dry mass, M 0 the initial dry mass and t the number of incubation days.
Fungal DNA was extracted from the incubated leaves using a PowerSoil Ò DNA Isolation Kit (MO BIO Laboratories, San Diego, USA) as described in Turunen et al. (2017). The yielded DNA sequences were analyzed using Quantitative Insights Into Microbial Ecology (QIIME) version 1.8.0 (Caporaso et al., 2010). The sequence library was split by samples and quality filtered for each sequence using default settings in QIIME. A total of 1,464,869 and 1,122,516 sequences were retained for bacteria and fungi, respectively. The sequences were clustered as operational taxonomic units (OTUs, 97% similarity) using the Usearch algorithm (Edgar, 2010). Due to the large number of OTUs, those that occurred in less than five sites were excluded from analyses.

Statistical analyses
We first explored the relationship of stream bryophyte cover and species richness to fine sediment cover by linear regression to quantify the influence of sedimentation on bryophytes. The data fulfilled the assumptions of normality and homoscedasticity of linear regression. Bryophyte cover percentages were arcine square root transformed prior to analysis. Significance level was set at P \ 0.05. Collinearity between fine sediments, bryophytes and other environmental variables were explored to ensure that no highly correlated variables (P \ 0.60) were included in the analysis (Online resource 2). We used random forest modeling (Breiman, 2001) to explore the importance and quantify the relationship of each of the biological response variables to bryophyte cover, fine sediments and other environmental variables. Random Forest models were used because they are a powerful tool for detecting and modeling complex relationships between many predictor variables (Cutler et al., 2007). Moreover, the models account for the effects of many predictor variables simultaneously. Therefore, we were able to examine the individual effect of each predictor separately after the effects of all the other predictors were accounted for (Hastie et al., 2001).
As biological response variables, we included community composition of leaf-decomposing fungi, benthic diatoms, and benthic macroinvertebrates, and leaf decomposition rate. To summarize variation in community composition, we ran Non-parametric Multidimensional Scaling (NMDS) by using R-package vegan (Oksanen et al., 2018) and function metaMDS and used the first and second axis scores as response variables in Random Forests. In metaMDS principal components rotate the final configuration so that the variance is maximized on the first axis scores. Diatom and macroinvertebrate NMDS species scores were used to describe the community change related to NMDS axes. For fungi we had only OTU data and the taxonomic composition could therefore not be described.
All predictor variables (Table 1) were first included in the Random Forest models and then the R-package Variable Selection Using Random Forests (VSURF, Genuer et al., 2015) was used to select the most meaningful environmental predictor variables for each biological response. As predictors for leaf decomposition we also included abundance, richness, and evenness of shredders in the leaf bags, and richness and evenness of fungi, as biotic variables that may influence decomposition (McKie & Malmqvist, 2009;Tolkkinen et al., 2013). To determine the importance of selected explanatory variables to each response variable, we calculated percentage increase in mean squared error (MSE) that describes the increase in model error when predictor variable values are randomly permuted relative to model performance with the original data. The larger the increase in error, the more important the variable is for the model. To visualize the relationship between key explanatory variables and response variables we used partial dependency plots (De'Ath 2007) which show the average trend in the response variable as a function of the focal explanatory variable, while keeping all other explanatory variables in the model constant. All analyses were performed using R software (version 3.5.1; R Core Team, 2018).

Results
The stream bryophyte cover varied from 1 to 84 and fine sediment cover from 10 to 100% (Table 1). Mean water temperature ranged from 4.7 to 10.6°C (mean of 7.0°C). Water color (a proxy for dissolved organic carbon concentration) ranged from 15 to 140 mg Pt l -1 (mean 55 mg Pt l -1 ) and total phosphorus (a proxy for trophic status) from 6 to 28 lg l -1 (mean 14 lg l -1 ) ( Table 1).
The most important variables related to leaf decomposition rate were shredder abundance (% increase in MSE = 12.3) (Fig. 2h), mean water temperature (MSE = 9.7) (Fig. 2i) and catchment area (MSE = 8.5) (Fig. 2j), explaining 40.3% of the variability. Decomposition rate increased sharply with increase in shredder abundance, until it reached a plateau at around 100 individuals. Decomposition rate increased also with increasing water temperature and catchment area.  diatom (a, b, c), macroinvertebrate (d, e, f), fungal (g) community structure (NMDS axis 1 scores) and leaf decomposition (h, i, j) to the most important environmental variables

Discussion
We found that bryophyte loss, i.e., reduction of both bryophyte abundance and richness, was related to increasing sedimentation. Despite the significant negative correlation between bryophyte abundance and fine sediment cover, some streams retained substantial bryophyte cover even when high amounts of fine sediment cover were detected, which enabled us to assess the relative importance of the two correlated variables to community structure and leaf decomposition. Our key finding was that bryophyte abundance was consistently more important than the forestryinduced sedimentation in explaining variation in biological community structure, a result that complements previous experimental evidence (Turunen et al., 2018). As expected, bryophyte cover was among the most important variables to diatom and macroinvertebrate community composition. Contrary to our expectations, however, fungal community composition was also strongly influenced by bryophyte cover. Leaf decomposition rate was mostly influenced by variation in shredder abundance and water temperature and, in contrast to our predictions, was not related to bryophyte cover.
Relationships of diatom and macroinvertebrate communities to environmental gradients Diatom species composition changed distinctly along the bryophyte cover gradient. Bryophytes form an attachment substrate for benthic diatoms and often support a distinct epiphytic assemblage from that on stony substrates (Rothfritz et al., 1997;Lim et al., 2001;Knapp & Lowe, 2009). Additionally, large mosses may shade the benthic habitat, which may also influence diatom community composition (Lange et al., 2011). Water color and temperature were also important determinants of diatom communities. In boreal streams, water color is generally correlated dissolved organic carbon (DOC) and iron concentrations in the stream water (Temnerud et al., 2014). Variation in water transparency influences benthic light conditions and thereby diatom communities (Lange et al., 2011). Diatom taxa generally show distinct thermal preferences (Weckström et al., 1997), and temperature is therefore one of the key determinants of diatom species composition.
The main variation in macroinvertebrate community composition was related to variation in water temperature, color, and catchment size. However, the variation in macroinvertebrate NMDS axis 2 scores was related to bryophyte abundance, highlighting that bryophytes do shape invertebrate community structure by providing abundant detrital food and habitat (Turunen et al., 2018). Indeed, many of the species that were associated with sites that had high bryophyte cover were detrivores (e.g., N. avicularis, N. pictetii, Leptophlebia sp., P. meyeri, and M. gelidum). Water temperature had a particularly large effect, accounting for a great share of macroinvertebrate community variability. Some of our study streams are strongly affected by groundwater input, being thus cold and stable environments that support unique invertebrate taxa (Ilmonen & Paasivirta, 2005), many of which are crenophilous (e.g., A. sulcicollis, N. pictetii, Scleroprocta sp., L. nigra) (Ilmonen & Paasivirta, 2005). Water color reflects the dominance of peatlands in stream catchments, and of humic substances in the stream water, which are known to be key determinants of macroinvertebrate assemblages in boreal streams (Aroviita et al., 2009).

Relationships of fungi and leaf decomposition to environmental gradients
Streams with abundant bryophytes typically sustain higher organic matter standing stock (Muotka & Laasonen, 2002;Turunen et al., 2018), which may cause differences in fungal community composition in streams that differ in bryophyte abundance. Fine sediment cover was related to variation in NMDS axis 2 scores, indicating a secondary role to fungal community composition. Contrary to many previous studies (e.g., Suberkropp & Chauvet, 1995;Tolkkinen et al., 2015), natural variation in water chemistry (pH and nutrients) had little influence on fungal communities, likely reflecting the low variation of water chemistry in our study streams. However, our water samples represented a snapshot of water quality during a low flow period (late summer/early autumn) and thus do not necessarily represent the water chemistry driving the species composition.
Leaf decomposition was not related to fine sediment cover or loss of bryophytes. In a recent mesocosm experiment, leaf decomposition rate was reduced in treatments with Fontinalis mosses because mosses supported abundant organic matter standing stocks, resulting in lower shredder aggregation into leaf packs in the presence of mosses (Turunen et al., 2018). In this study, shredder abundance in leaf bags was positively related to decomposition rate, as was also observed by Mustonen et al. (2016) and Turunen et al. (2018). The positive relation of leaf decomposition rate to water temperature is frequently observed in both field (e.g., Martínez et al., 2014) and laboratory (e.g., Ferreira & Chauvet, 2011) experiments. Fungal richness or evenness were both unrelated to decomposition rate (see also Dang et al., 2005;Tolkkinen et al., 2013). Our results thus support the conjecture that shredders contribute relatively more to leaf decomposition than microbes, particularly in streams where water quality is not severely impaired by anthropogenic activities (Graça et al. 2001;Hieber & Gessner, 2002;Taylor & Andrushchenko, 2014).

Conclusions
A major challenge to biodiversity conservation and restoration is that the loss of key species may have cascading ecological effects on all trophic levels and thus fundamentally change ecosystem functioning (Tronstad et al., 2010;Ellis et al., 2011). Our results suggest that bryophytes might have such a role in boreal stream food webs. Increasing sedimentation from forestry practices, especially from drainage, is a key stressor inducing loss of bryophyte abundance and species diversity. However, the consistently stronger relationship of the community composition and leaf decomposition to bryophyte cover than to fine sediments suggests that it is the disappearance of bryophytes that links a large portion of the negative impacts of fine sediments to other communities in boreal streams. This result suggests that restoring sediment-stressed streams by the addition of coarse substrate particles (boulders, wood) may enhance the recolonization of stream bryophytes and thus aid the ecological recovery of stream ecosystems (Turunen et al., 2017). Abundant bryophytes have also been found to decrease temporal community variability and bryophytes may therefore be important also in protecting boreal stream assemblages from stochastic climatic variability (Huttunen et al., 2017). Finally, the strong effect of water temperature on species composition and leaf decomposition rate highlights the potential sensitivity of headwater stream biodiversity to changes in water temperature. There is evidence of temperature increase caused by global warming in boreal spring-fed streams (Jyväsjärvi et al., 2015) and, together with potential hydrological changes (Mustonen et al., 2018), this may threaten the unique biodiversity of these ecosystems.