Drainage‐induced browning causes both loss and change of benthic biodiversity in headwater streams

Concentrations of dissolved organic carbon (DOC) have increased over the past few decades, causing freshwater browning. Impacts of browning on biodiversity have been little studied, despite many of the individual stressors associated with browning being known to control freshwater communities. We explored the responses of benthic invertebrates along a wide gradient of DOC concentrations (3.6–27 mg L−1) in 63 boreal streams variously impacted by peatland drainage or peat production. DOC was a prime determinant of macroinvertebrate diversity and abundance, with the strongest negative response in algal scrapers. Threshold indicator taxa analysis indicated an abrupt community change at 12–13 mg DOC L−1, with only four taxa increasing, while 13 taxa decreased along the gradient. Our findings of both a gradual loss and abrupt change of biodiversity along a browning gradient provide a benchmark against which changes to stream biodiversity relative to the predicted browning trend can be gauged.

terrestrial productivity (Solomon et al. 2015). Regardless of the mechanisms causing browning, its ecological consequences are less studied, particularly in lotic ecosystems and benthic habitats therein. Lake studies have generally shown a unimodal relationship between browning (i.e., DOC concentrations) and lake productivity: at low concentrations, productivity is nutrient-limited, then reaches a plateau at around 10 mg DOC L À1 . At higher concentrations, lentic production becomes light-limited (e.g., Bergström and Karlsson 2019). Most evidence for this pattern comes from lake pelagic communities but a similar pattern was recently documented also for benthic algae (Fork et al. 2020).
Corresponding studies on primary consumers other than zooplankton are rare. Jonsson et al. (2015) manipulated browning and warming in experimental ponds and detected a change in community structure toward larger taxa. This response was, however, much stronger for warming and any effect of browning was only detectable when water temperature also increased (Jonsson et al. 2015). In contrast, Vasconcelos et al. (2016) showed in another pond experiment that browning overwhelmed any effects of warming by reducing benthic algal production and invertebrate biomass. It thus appears that browning may exert at least some, and sometimes quite strong, control over benthic consumers but this effect is mainly indirect, reflecting the negative effect of browning on benthic algal productivity. It is also possible that algal quality, not quantity, is decisive; the nutritional quality of benthic algae decreases already at the early stages of browning, even if algal quantity shows little response . Browning is associated with several concurrently changing factors and the individual effects of, for example, increased DOC and nutrients (e.g., Bergström and Karlsson 2019) cannot be distinguished without experimentation. Nevertheless, broad-scale field surveys can help us understand and predict biotic responses to browning (see Solomon et al. 2015).
Impacts of browning on freshwater biodiversity have been surprisingly little studied. Kesti et al. (2022) showed that chironomid midges were more abundant in humic lakes and mayflies and stoneflies in clearwater lakes. However, dividing sites into a strict dichotomy may not be optimal if community responses along a stressor gradient are continuous. Whether biological communities change gradually or exhibit an abrupt change at a distinct point of a gradient is a long-standing debate in ecology. Although current research considers these two models of community organization as endpoints of a continuum (Liautaud et al. 2019), sometimes even small changes in the environment can induce a major, potentially irreversible, regime shift in ecosystems (Andersen et al. 2009;King and Baker 2011). Being able to predict if (and when) ecosystems are approaching such critical thresholds should be particularly relevant for processes that occur over a "human timescale" of years and decades (Cooper et al. 2020) and are likely to continue at the same, or accelerating, pace in the future.
We explored the responses of benthic invertebrate communities along a DOC gradient in 63 boreal streams impacted by various degrees of peatland drainage or peat production. We first examined if diversity (taxonomic richness; feeding trait composition) and invertebrate community composition exhibited a gradual change along the browning gradient or an abrupt shift at a certain point of the gradient, indicative of an ecological threshold in community response. We then examined which invertebrate taxa best indicated the presence of such a threshold (if any), therefore serving as indicators of a browning-induced shift in benthic community organization.

Sampling sites
We selected 63 headwater streams (orders 1-2) to represent a gradient of DOC concentrations from 3.6 to 27 mg L À1 (Fig. S1). The study sites constitute headwaters of two neighboring watersheds in north-eastern Finland, River Iijoki and River Koutajoki. The study area is characterized by extensive peatlands and mixed forests. Each study site consisted of a continuous 50-m long riffle section where all sampling was conducted.

Benthic invertebrate sampling
Benthic macroinvertebrates were sampled in late September 2019 (53 sites) or 2020 (10 sites) with a 2-min kicknet (Ø 500-μm) sample at each site. A sample of this size covers about 1.3 m 2 of the stream bed and is known to capture about 75% of all species present in a riffle (Mykrä et al. 2006). Invertebrates were identified to species or genus level.

Calculation of community metrics
We calculated three Hill numbers (q1, q2, q3) on the macroinvertebrate data. However, as other Hill numbers closely reflected the patterns in q1 (species richness), we only focus on richness hereafter.
We characterized the feeding habits of invertebrate taxa using a European trait database (Schmidt-Kloiber and Hering 2015). The propensity of a taxon to belong to each of four feeding trait categories (scrapers, shedders, gatherer-collectors, predators) was quantified by an affinity score from 0.1 to 1.0 based on a fuzzy-coding approach (Chevenet et al. 1994). Taxa abundances were then weighted by affinity scores.

Environmental measurements
We collected water samples simultaneously with macroinvertebrate sampling and they were analyzed for DOC (mg L À1 ), water color (mg Pt L À1 ), total phosphorus (μg L À1 ), nitrate-N (NO 3 ; μg L À1 ), and pH following national standards (National Board of Waters 1981) (see Table S1). DOC concentration was measured from filtered (Ø 0.45 μm Whatman GF/F) water samples by infrared spectrometry with Shimadzu TOC-VCPH analyzer (Shimadzu Scientific Instruments).
Water temperature ( C) was recorded at 1-h intervals by HOBO Pendant data loggers for 5 weeks prior to invertebrate sampling. Mean substratum size was determined as the weighted average of particle size classes in 10 randomly distributed 0.25 m 2 plots using a modified Wentworth scale from 1 (clay/silt 0.001-0.07 mm) to 10 (large boulder to bedrock >512 mm). Substrate diversity was quantified by calculating Simpson diversity index based on the relative cover of each Wentworth class in each quadrat. Moss cover (%) was estimated in 20 0.25 m 2 plots distributed randomly across the study section. Current velocity (cm s À1 ; Schiltknecht MiniAir20) and water depth (cm) were measured at three positions along each of five randomly placed transects perpendicular to the flow (Table S1). All data are available in Brüsecke et al. (2022).

Data analyses
We first classified streams into three subgroups based on DOC-concentration: low (< 10 mg L À1 ; n = 19), moderate (10-15 mg L À1 ; n = 19) and high-DOC (> 15 mg L À1 ; n = 25), roughly corresponding to 20%, 75%, and 98% of water color observations in ca. 4500 Finnish streams . Environmental gradients were then visualized with principal component analysis (PCA), separately for habitat and water chemistry variables. Differences among the DOC groups were tested with permutational multivariate ANOVA (PERMANOVA) using adonis function of the vegan R package (Oksanen et al. 2018). PERMANOVA was run with 9999 permutations using the Euclidean distances on zstandardized data with zero mean and unit variance.
Least squares multiple linear regressions were used to examine the relationships between water chemistry, habitat variables, and univariate community metrics (species richness; scraper, gatherer-collector, predator and shredder abundance). Abundance-based metrics were log 10 -transformed prior to analyses. For each environmental predictor, we estimated the variance inflation factor (VIF). The VIF values suggested some degree of multicollinearity (1.5-2.7), but all values were well below 5, which has been suggested as a limit for alarming multicollinearity (Hair et al. 2016). The most parsimonious model was selected based on the empirical variable selection (REVS) approach (Goodenough et al. 2012) using the regsubsets function of the leaps package (Lumley 2020) in R statistical software (R Core Team 2021). The model with the lowest Bayesian information criterion was considered the most parsimonious one. To further compare the relative importance of model predictors, we calculated variable importance (% of R 2 ) using the "lmg" method in the relaimpo R package (Grömping 2006) and rescaled the values to sum up to 100%. Confidence intervals (95%) for variable importance values were obtained by bootstrapping (n = 9999).
Macroinvertebrate community change along the DOC gradient was determined using Threshold Indicator Taxa Analysis (TITAN 2.4.1 R package; Baker et al. 2015). TITAN combines indicator species analysis (Dufrêne and Legendre 1997) with nonparametric change-point analysis (Baker and King 2010) and enables identification of change points in the frequency and abundance of individual taxa. Further, it examines if multiple taxa have synchronous responses along the gradient. TITAN calculates indicator values (IndVal) for all taxa and all possible change points along the gradient of interest and uses permutation tests to assess the level of uncertainty in these scores. Permuted IndVal scores are standardized as Z scores and summed for positive ("Z+") and negative ("ZÀ") response values for each possible change point. Peaks in Z values highlight points around which a disproportionate number of taxa exhibits changes in occurrence and relative abundance, thus representing a potential community threshold (Baker and King 2010). We used bootstrapping (n = 1000) to estimate uncertainty around the change points for each indicator taxa. A taxon was assigned as a positive or negative responder if: (1) ≥ 95% of bootstrapped runs were significantly different from random distribution (at p ≤ 0.05, high reliability), and (2) if the change in frequency and abundance of the taxon was in the same direction for ≥ 95% of the bootstrapped runs (at p < 0.05, high purity).

Environmental conditions
The first two principal components of the instream PCA explained 56% of total variance. The first axis represented a size gradient with depth, width and substrate size increasing along PCA 1 (Fig. 1a). The DOC groups did not differ in instream characteristics (PERMANOVA: F 2,60 = 1.01; p = 0.396).

Benthic invertebrate diversity along the browning gradient
Invertebrate richness varied between eight and 40 taxa per site and was mainly (28.3%) explained by DOC (negative relationship) and water temperature (positive relationship) ( Fig. 2; Table 1). DOC contributed 76% of all explained variation in the multiple regression model (Fig. S2b).

Functional feeding groups
DOC explained 33.6% of variation in scraper abundance (Fig. 3a). Residual variation was positively related to water temperature and stream width (Table 1). DOC was clearly the predominant predictor (75% of explained variation) in the multiple regression model (Fig. S2b).
Gatherer-collector abundance (Fig. 3b) was explained in nearly equal amounts by NO 3 -N, DOC and stream width ( Fig. S2c; Table 1). For predatory invertebrates, the univariate relationship with DOC bordered at significance (Fig. 3c) but when all environmental variables were included, predator abundance was only related to stream width (positive) and depth (negative) (Fig. S2d; Table 1). Shredder abundance showed no relationship to DOC (Fig. 3d).

Discussion
The ongoing browning of freshwaters raises concerns about the ecological consequences of this trend. Our data showed a gradual reduction in invertebrate richness with increasing DOC, and an abrupt shift in benthic macroinvertebrate community composition when DOC concentrations reached 12-13 mg L À1 . Importantly, algal scrapers were the only feeding group with a strong negative response to DOC, likely related to an indirect effect of DOC on invertebrates via their key food source, benthic biofilm, which is of lower nutritional quality at higher DOC levels . The importance of the indirect linkage between invertebrate consumers and browning via food resources is further supported by the observation that other feeding groups (whose food resources should be less affected by browning) showed a weaker (gathering-collectors; predators) or no response (shredders) to increasing DOC. Jyväsjärvi et al. (2022) combined a mesocosm experiment and field survey to show a drastic reduction in algal (particularly diatom) biomass and bioavailability of polyunsaturated fatty acids (PUFA) with water browning. At the same time, the availability of terrestrial-derived long-chain saturated fatty acids poor in PUFAs (Taipale et al. 2016) increased. These observations combined with those of this  study indicate that the loss of benthic biodiversity was likely caused by impairment of the nutritional base of invertebrates as a result of browning. The shift in community composition at 12-13 mg DOC L À1 resulted from the loss or decline of many sensitive species, while only a few taxa became more abundant along the gradient. The taxa that benefited from browning are generalists known to cope in stressful environments, for example, the isopod Asellus aquaticus (see O'Callaghan et al. 2019). Asellus could also have benefited from reduced predation risk at high-DOC sites but most of the invertebrate predators in our data do not feed much on Asellus (e.g., Tikkanen et al. 1997). We did not quantify fish abundances at our study sites but some fish species, particularly burbot, bullhead and minnow, are known to occur commonly in peatland-dominated, dark-colored streams (Sutela et al. 2021).
Most of the sites that exceeded 15 mg DOC L À1 drain excessive peatlands and those exceeding 20 mg DOC L À1 were downstream of heavily ditched peatlands or peat production areas. Thus, the key mechanism causing stream browning in our study area is linked to forestry-related land use, particularly peatland drainage. The same observation has been made for other boreal streams and lakes. For example, using a historical data set from Swedish lakes, Kritzberg (2017) identified transition from agricultural to forestry-related land use as a major factor causing browning (see alsoŠkerlep et al. 2020). More specifically, Nieminen et al. (2021) reported ca. 14 mg L À1 higher concentrations of total organic carbon in streams with drained than undrained catchments in Sweden and Finland. These results support the view that increase in forest biomass due to intensive peatland drainage has contributed strongly to water browning in boreal regions. Peatland drainage can also enhance deposition of fine sediments, thus potentially altering benthic biodiversity  (Ramchunder et al. 2012), but this effect seemed minor in our study as substrate size was not an important driver of biodiversity. Elevated loads of organic matter can also increase aluminum concentration, causing direct toxicity for aquatic organisms (Helmer et al. 1990;Vuori et al. 1998). As we did not measure aluminum in our study streams, we cannot exclude the potential added effect of this contaminant on benthic biodiversity in strongly drained watersheds. The few previous studies that have addressed the impacts of peatland management on benthic biodiversity have generally shown that drainage-induced changes to stream habitat and/or water quality reduce invertebrate diversity in recipient streams (Ramchunder et al. 2012;Brown et al. 2019;Rajakallio et al. 2021). Similarly, in a study of 18 streams in northern Sweden, Jonsson et al. (2017) concluded that DOC and pH were the main drivers of benthic invertebrate richness and community composition. However, a direct comparison between their study and ours is complicated because their DOC gradient (9.8-42.4 mg L À1 ) only included one site where DOC concentration was below the threshold value detected in our study. None of these studies focused directly on browning and were therefore not conducted along an anthropogenically caused browning gradient. Thus, our finding of both a gradual biodiversity loss and an abrupt biodiversity change along a browning gradient are novel, providing a benchmark against which forthcoming changes to stream biodiversity relative to the predicted browning trend can be gauged. Thus far most long-term studies have shown a consistent trend of browning (e.g., Asmala et al. 2019;Lapierre et al. 2021) and some studies have used various climate scenarios to predict that freshwater ecosystems will continue to brownify for many decades to come. For example, Weyhenmeyer and Karlsson (2009) showed that, in the worst-case scenario, absorbance in Swedish lakes and rivers will increase on average by a factor of 1.3 until 2030. Relating this prediction to our data, the loss of invertebrate biodiversity in boreal streams impacted by land drainage seems inevitable. Importantly, the order of species loss is usually not random but depends on species' traits in relation to the type of disturbance (Jonsson et al. 2002). In our case, a few tolerant taxa benefited from browning whereas more vulnerable taxa, many of which are algal scrapers, were greatly reduced. Such a major shift in community composition is unlikely to go without substantial ecosystem-level implications as benthic algae are then largely released from top-down regulation by invertebrate grazers.
Our study used a space-for-time substitution design to assess the community-wide consequences of browning on stream invertebrates. The use of such designs was recently criticized by Stetler et al. (2021), because they rely on the assumption that drivers of ecological change through time are the same as those that drive changes in space. Using continentalscale spatial data and a long-term regional study on nutrient and DOC concentrations, Stetler et al. (2021) showed that DOC and nutrients are coupled spatially but not temporally, indicating that, across decadal time intervals, browning may not be associated with increased nutrient concentrations. In their study region, temporal variation in DOC mainly reflects recovery from acidification whereas spatial variation is more site-specific. The situation may be different in boreal streams where the same factor(s) related to increased forest cover (Kritzberg 2017;Škerlep et al. 2020) and peatland drainage (Nieminen et al. 2021) are causing browning across both time and space. Furthermore, as the forest-rich (and forestry-intensive) countries are largely responding to intensifying bioeconomy by further boosting their wood production, browning of freshwater ecosystems is likely to continue. On a more positive side, restoration of drained peatlands is an increasing activity (e.g., Humpenöder et al. 2020) and while little is known about the consequences of peatland restoration on stream ecosystems (but see Ramchunder et al. 2012) there is evidence that restoration by blocking of drainage ditches initially increases nutrient loads (Koskinen et al. 2017), causing potentially harmful short-term impacts on stream organisms and food webs. On a longer time scale, peatland restoration could be expected to reduce DOC loads but longterm responses of freshwater ecosystems to peatland restoration are largely unknown, thus representing a critical research gap for a sustainable management of streams and their watersheds in peatland-dominated catchments.