The Impact of Forest Fertilization on the Ecological Quality of Two Hemiboreal Streams

: The present study aimed to detect any changes in concentrations of nutrients and evaluate the impact on the quality of two hemiboreal streams that collect a discharge from two fertilized Scots pine stands. In 2017, nitrogen-containing mineral fertilizer was spread in pine stands on mineral soil located near the ﬁrst stream. In 2018, potassium containing wood ash was spread in pine stands on organic soil near the second stream. From 2017 to 2020, surveys of physico-chemical parameters, diatoms, macrophytes, and macroinvertebrates were performed to determine the possible effects of fertilization on the ecological quality of the streams. A control site upstream of the fertilized forest stand and a treatment site downstream of the fertilized forest stand was monitored at each stream. Water quality indices, chemical parameters of surface water, and indicator species analysis showed no short-term impact of forest soil improvement with wood ash and ammonium nitrate. We found no clear patterns before and after the fertilization events in both streams, although we did observe inter- and intra-annual differences in aquatic biota and stream ecological quality mainly caused by local environmental factors.


Introduction
Forest management practices include drainage, sanitary cuts, pre-commercial thinning, and proper regeneration, including implementing forest soil preparation methods and selecting planting material. In addition, forest soil fertilization-depending on the fertilizer type and site conditions-is recognized as an effective measure to enhance tree growth [1,2]. Forest soil fertilization has been practiced in Sweden since the 1950s, with 190,000 ha per year fertilized with nitrogen-containing fertilizer in the mid-1970s to improve the growth conditions of forests. However, this practice declined significantly in the 1990s due to an increased ecological focus in Swedish forestry [3]. Due to soil acidification, it is vital to focus research on forest soil treatment with wood ash [4] and determine its subsequent effects on the chemistry of soil and water [5][6][7][8]. In the 1960s, fertilization experiments of Norway spruce were conducted in Denmark, but since 1980, experiments mainly focused on counteracting decline of forest. However, forest soil fertilization is no longer practiced in Denmark [9]. On the contrary, it is estimated that annual fertilization of 5000-10,000 ha forest in Norway would contribute to the CO 2 sequestration in 10 years of 0.14-0.27 million tonnes/year [10].
Local studies have reported the positive impact of wood ash on the radial increment of hybrid aspen plantation [11] in the first 1-2 years and additional radial increment in Table 1. Characteristics and treatments of the studied forest stands located at two hemiboreal streams in Latvia.

Study Sites
The first experimental site was located next to the Age (No.1. in Table 1), Limbazi Parish, in the northern part of Latvia. The experiment was carried out in a 93-year-old stand of Scots pine (Pinus sylvestris L.) on drained organic soil ( Figure 1 and Table 1).
According to the Latvian forest site type classification system, the forest type of the stand is Myrtillosa turf. mel.-mesoeutrophic site [39]. The second experimental site was situated next to the Rusinupe (No.2. in Table 1), Vecumnieki Parish, in the central part of Latvia. The experiment was conducted in a 61-year-old stand of Scots pine on drained mineral soil. The forest type of the stand is Myrtillosa-mesotrophic site. The annual mean air temperature during 1961-2010 in the areas of Age and Rusinupe was 6.2 and 5.9 °C, respectively. The annual mean precipitation during 1961-2010 in these same areas was 667 and 704 mm y −1 , respectively [40].   [41]. The main characteristics of the soil are summarized in Table  2. According to the Latvian forest site type classification system, the forest type of the stand is Myrtillosa turf. mel.-mesoeutrophic site [39]. The second experimental site was situated next to the Rusinupe (No.2. in Table 1), Vecumnieki Parish, in the central part of Latvia. The experiment was conducted in a 61-year-old stand of Scots pine on drained mineral soil. The forest type of the stand is Myrtillosa-mesotrophic site. The annual mean air temperature during 1961-2010 in the areas of Age and Rusinupe was 6.2 and 5.9 • C, respectively. The annual mean precipitation during 1961-2010 in these same areas was 667 and 704 mm y −1 , respectively [40].
The dominant soil type in the forest stand No. 1 is Histosol, and Podzol is dominant in the forest stand No. 2 [41]. The main characteristics of the soil are summarized in Table 2.

Experimental Design and Treatments
Samples were collected upstream and downstream of the experiment site to obtain data on control and impact [42]. The Age is a medium-sized rhithral stream with silt, pebbles, Forests 2022, 13, 196 4 of 15 and rocks on the riverbed. The monitoring was conducted before and after the inflow of the tributary, which collects a discharge and a runoff of the forest stands fertilized with WA, into the Age. Rusinupe is a small potamal stream with sand and coarse particulate organic matter on the riverbed. The monitoring points were established upstream and downstream of the forest stand, where NH 4 NO 3 was spread. In Age, a beaver dam next to the inflow of the tributary was observed throughout 2018. In Rusinupe, a beaver dam was detected upstream in the autumn of 2018. From 2017 to 2020, macrophytes, macroinvertebrates, and phytobenthos were monitored once, twice, and three times, respectively, per year.

Surface Water Sample Collection and Analysis
Surface water samples were collected twice a month from 2017 to 2019. At the first study site, samples were collected from the Age tributary before the inflow into the Age. At the second study site (the Rusinupe), samples were collected upstream and downstream of the fertilized forest stand. Chemical analyses of the surface water samples were carried out in an accredited laboratory of the Latvian Environment, Geology and Meteorology Centre.

Sampling and Analysis of Biological Quality Elements
Because the streams contained different types of riverbed substrates, phytobenthos ( Figure S1) were collected from branchwood and pebbles. Slides were prepared according to the European standard NF EN 13946 for the microscopic study of diatoms. Up to 400 frustules per sample were counted and identified at the species level. The Trophic Diatom Index (TDI20 and TDI100), Watanabe's Index (WAT), and Specific Polluosensitivity Index (IPS) were calculated using OMNIDIA [43] and used to assess ecological quality [44].
Macrophyte ( Figure S2) assemblages in a 100 m stretch were described once a year during the vegetation season. The diversity and abundance of the detected species were determined using a nine-point grading scale. Macrophytes were recorded and identified at the species level or the lowest practical taxonomic level.
The Macrophyte Index for Rivers (MIR) was calculated as described by Szoszkiewicz et al. [45] to estimate the ecological quality of streams. To assign the MIR values to the ecological quality classes as required by the European Union (EU) Water Framework Directive, the calculated MIR was standardized using the ecological quality ratio (EQR) [46].
Samples of macrozoobenthos ( Figure S3) were collected twice a year in spring and autumn. Five replicate samples were collected at each site, in proportion to microhabitat occurrence. The macroinvertebrates were primarily identified at the species or genus level. Oligochaetes were not identified further, whereas Dipetra were identified to the family level. Macroinvertebrate indices were calculated using Asterics v. 4.04 [47]. A total of 102 numerically and ecologically suitable metrics were selected for further analysis.
In the Age, ecological quality based on macroinvertebrates was assessed using the Latvian Macroinvertebrate Index (LMI) [48]. In Rusinupe, we applied a method used in Estonia [44] as there is currently no appropriate ecological quality assessment method based on macroinvertebrates for small streams in Latvia.

Statistical Methods
Principal component analysis (PCA)-an ordination analysis-was performed using PC-ORD 5.0 software [49] to determine the relationship between the selected metrics and samples at each stream. Benthic invertebrate taxa and metric data were log-transformed before conducting the PCA. The abundance of benthic invertebrate taxa in Age and Rusinupe was measured via indicator species analysis (ISA) using the PC-ORD 5.0 software [49]. A common aspect of ecological research is identifying species associated with, or indicative of, groups of samples, including studies of environmental management. In-depth analyses are often required to identify species indicative of specific groups. ISA allows statistically rigorous assessments of these indicator species [50]. Taxonomic data were grouped according to three categorical variables, namely site position to the treatment (upstream or downstream), season (spring, autumn), and year (2017 to 2020). We performed the Monte Carlo test of significance of the observed maximum indicator value for species, with 4999 permutations. Indicator values were calculated according to the method described by Dufrene and Legendre [51]. The non-parametric Wilcoxon rank sum test with continuity correction was used to compare the chemical parameters of water samples collected at the control and treatment sites. The test was conducted at a 95% confidence level using the R program [52].

Chemistry of Surface Water
Following the treatment with wood ash, the average pH of water samples from the Age tributary changed from 6.88 ± 0.10 to 6.91 ± 0.06 (p < 0.01). The average K concentration was slightly higher two years after the treatment with wood ash than during the one-year period before fertilization-(from 0.56 ± 0.13 to 0.69 ± 0.08) ( Figure 2).
Forests 2022, 13, x FOR PEER REVIEW 5 of 16 or indicative of, groups of samples, including studies of environmental management. Indepth analyses are often required to identify species indicative of specific groups. ISA allows statistically rigorous assessments of these indicator species [50]. Taxonomic data were grouped according to three categorical variables, namely site position to the treatment (upstream or downstream), season (spring, autumn), and year (2017 to 2020). We performed the Monte Carlo test of significance of the observed maximum indicator value for species, with 4999 permutations. Indicator values were calculated according to the method described by Dufrene and Legendre [51]. The non-parametric Wilcoxon rank sum test with continuity correction was used to compare the chemical parameters of water samples collected at the control and treatment sites. The test was conducted at a 95% confidence level using the R program [52].

Chemistry of Surface Water
Following the treatment with wood ash, the average pH of water samples from the Age tributary changed from 6.88 ± 0.10 to 6.91 ± 0.06 (p < 0.01). The average K concentration was slightly higher two years after the treatment with wood ash than during the oneyear period before fertilization-(from 0.56 ± 0.13 to 0.69 ± 0.08) ( Figure 2). We detected a statistically significant increase in PTOT concentration (from 0.03 ± 0.01 to 0.04 ± 0.00, p < 0.05) and a significant decrease of dissolved organic carbon (DOC) concentration (from 52.96 ± 2.46 to 41.62 ± 2.19, p < 0.01) after the application of wood ash in the Age catchment area. The average pH of water samples from the Age tributary ranged from 7.29 ± 0.05 (upstream) to 7.29 ± 0.06 (downstream). The average NTOT concentration did not indicate any nutrient loading of the applied NH4NO3 from the fertilized forest stand in the Rusinupe catchment area. The concentration of NTOT was slightly lower downstream of the fertilized forest stand (1.12 ± 0.06 control and 1.00 ± 0.06 fertilized). However, a significantly lower (p < 0.05) concentration of K (0.92 ± 0.03 control and 0.83 ± 0.03 fertilized) was detected in surface water downstream of the treatment site in Rusinupe. We detected a statistically significant increase in P TOT concentration (from 0.03 ± 0.01 to 0.04 ± 0.00, p < 0.05) and a significant decrease of dissolved organic carbon (DOC) concentration (from 52.96 ± 2.46 to 41.62 ± 2.19, p < 0.01) after the application of wood ash in the Age catchment area. The average pH of water samples from the Age tributary ranged from 7.29 ± 0.05 (upstream) to 7.29 ± 0.06 (downstream). The average N TOT concentration did not indicate any nutrient loading of the applied NH 4 NO 3 from the fertilized forest stand in the Rusinupe catchment area. The concentration of N TOT was slightly lower downstream of the fertilized forest stand (1.12 ± 0.06 control and 1.00 ± 0.06 fertilized). However, a significantly lower (p < 0.05) concentration of K (0.92 ± 0.03 control and 0.83 ± 0.03 fertilized) was detected in surface water downstream of the treatment site in Rusinupe.

Macrophytes
A total of 44 macrophyte taxa were found in the stretches of Age and Rusinupe. Of these, 39 taxa were identified to the species level and five taxa to the genus level. The number of species varied from nine macrophyte taxa at the control site in Rusinupe to 18 taxa at the experimental site in Age. Vegetation cover varied from 10% to 40% at the Age control site and from 20% to 50% at the Age treatment site. In the Rusinupe, vegetation cover varied from 10% to 15% upstream and from 70% to 80% downstream of the fertilization site. The species most frequently found at both Age sites were Nuphar lutea, Sparganium emersum, and Fontinalis antipyretica. In the Rusinupe, Sparganium emersum, Lemna minor, Cicuta virosa, and Elodea canadensis were the dominant species at the control

Phytobenthos
We detected 124 and 147 diatom taxa in the Age and the Rusinupe, respectively. The most frequently occurring taxa in the Age were Amphora pediculus, Cocconeis placentula, Gomphonema parvulum, Meridion circulare, Planothidium lanceolatum, Platessa conspicua, Rhoicosphenia abbreviate, and Stauroneis anceps. In the Rusinupe, the most frequently occurring taxa were Amphora pediculus, Cocconeis placentula, Navicula cryptocephala, and Meridion circulare. Ecological quality varied from medium to high in both streams and stretches, and the TDI, WAT, and IPS values varied seasonally. At the Age control site, TDI varied from medium to high quality, WAT from medium to high quality, and IPS from good to high quality. At the Age treatment site, TDI fluctuated from medium to good quality, WAT from good to high quality, and IPS good to high quality.
In the Rusinupe, TDI was in the range of medium to high quality both at the control and treatment stretches. The WAT index showed good to high quality at the control site and moderate to high quality at the treatment site. The IPS values indicated good to high quality both at the control and treatment stretches.

Macroinvertebrates
We found a total of 169 macroinvertebrate taxa, of which 133 were in Age and 104 in Rusinupe. Caddisflies (Trichoptera) was the order with the highest species richness, with 30 taxa in the Age and 22 taxa in the Rusinupe. The macroinvertebrate abundance varied annually and seasonally, from 173 to 1685 specimens per sample in the Age and from 414 to 4043 in the Rusinupe. We observed higher macroinvertebrate abundance in spring in the Age and autumn in the Rusinupe.
In the Age, the LMI indicated moderate (in 2017 and 2018) to good (in 2019 and 2020) ecological quality in spring, but moderate ecological quality (in all years) in autumn. In the Rusinupe, the ecological quality also fluctuated from medium to good, reaching the highest values in spring of 2018 at both sites, in addition to in spring and autumn of 2020, at the experimental sites. Macroinvertebrate communities in the Age were more heterogeneous compared to those in the Rusinupe. Bivalvia, Diptera, and Crustacea were the dominant taxa and comprised 50% to 90% of the total number of benthic organisms found in the Rusinupe. In the Age, these three taxa constituted only approximately 30% to 60% of the total number of organisms; moreover, crustaceans were not among the dominating taxa. In the Age, a slightly stronger dominance of Diptera and Trichoptera was observed at the upstream site than at the downstream site. The number of Bivalvia increased upstream of the Rusinupe treatment site but decreased downstream of the site.

Indicator Species Analysis of the Age
Altogether, 10 indicator species showed significant interannual variability (p < 0.05). Most of these met the criteria for indicator species in 2018 (one species each from the Bivalvia, Diptera, Coleoptera, and Lepidoptera taxonomic groups and two Gastropoda species). Of the remaining, one was designated in 2017 (Ephemeroptera), one in 2019 (Odonata), and two in 2020 (Hydrachnidia and Ephemeroptera). Two bivalve species were significant indicator species for the experimental site, and three were for the control site (one Diptera and two free-living caddisfly species). However, seasonal variability showed the largest number of indicator species in spring (eight in total): four Trichoptera species and one each from the Simuliidae, Coleoptera, Plecoptera, and Hydrachnidia taxonomic groups. Four indicator species were designated in autumn: two Trichoptera species and one species each from Ephemeroptera and Plecoptera (Table 3). PCA of benthic invertebrate metrics (7 of 102 metrics were selected having 70% of the fitted range) showed an impact of season on calculated metric value distribution along the PCA axes (axes 1 and 2 explained 57.3% and 31.9% of the total data variance, respectively). From axis 1, the most significant correlation was found for EPT (%) abundance classes in spring and abundance of Trichoptera larvae in spring; however, from axis 2, total abundance was found in autumn.
The first two axes of the phytobenthos, macrophyte, and selected benthic invertebrate metric PCA explained 50.5% of the total data variance (29.6% and 20.9% for axes 1 and 2, respectively) in the Age. Significant correlations with PCA ordination axes were found for all biological quality element (BQE) group metrics. Positive correlations with axis 1 were found for several diatom metrics and the EQR of macrophyte MIR index, whereas negative correlations were found for several macrozoobenthos metrics. The highest positive correlation coefficients with axis 1 were found for IPS in autumn, whereas TDI in autumn and The highest positive correlation coefficient was found for IPS in autumn and high negative correlation coefficients with axis 1 were found for total coverage of macrophytes, abundance classes (%) of Ephemeroptera, Plecoptera, and Trichoptera (in spring), and abundance of Trichoptera (in spring). Furthermore, TDI in autumn and total abundance of benthic invertebrates (in autumn) had the highest positive correlation with axis 2.
The highest positive correlation coefficients with axis 1 were found for total coverage of macrophytes, abundance classes (%) of Ephemeroptera, Plecoptera, and Trichoptera (in spring), and abundance of Trichoptera (in spring); however, correlation coefficients were negative for IPS in autumn. Furthermore, TDI in autumn and total abundance of benthic invertebrates (in autumn) had the highest positive correlation with axis 2. The abundance of Bivalvia, number of Ephemeroptera, Plecoptera, Trichoptera, Coleoptera, Bivalvia, and Odonata (EPTCBO) taxa, and number of macrophyte species also correlated positively with axis 2. These metrics indicated higher abundances in 2018 at both the upstream and downstream sampling sites, forming a separate group from other samples, thus showing interannual differences in communities of biological quality elements. The upstream site also differed from other samples in 2017 and showed a negative correlation with axis 2 (Figure 3). the PCA axes (axes 1 and 2 explained 57.3% and 31.9% of the total data variance, respectively). From axis 1, the most significant correlation was found for EPT (%) abundance classes in spring and abundance of Trichoptera larvae in spring; however, from axis 2, total abundance was found in autumn.
The first two axes of the phytobenthos, macrophyte, and selected benthic invertebrate metric PCA explained 50.5% of the total data variance (29.6% and 20.9% for axes 1 and 2, respectively) in the Age. Significant correlations with PCA ordination axes were found for all biological quality element (BQE) group metrics. Positive correlations with axis 1 were found for several diatom metrics and the EQR of macrophyte MIR index, whereas negative correlations were found for several macrozoobenthos metrics. The highest positive correlation coefficients with axis 1 were found for IPS in autumn, whereas TDI in autumn and total abundance of benthic invertebrates (in autumn) had the highest positive correlation with axis 2.
The highest positive correlation coefficient was found for IPS in autumn and high negative correlation coefficients with axis 1 were found for total coverage of macrophytes, abundance classes (%) of Ephemeroptera, Plecoptera, and Trichoptera (in spring), and abundance of Trichoptera (in spring). Furthermore, TDI in autumn and total abundance of benthic invertebrates (in autumn) had the highest positive correlation with axis 2.
The highest positive correlation coefficients with axis 1 were found for total coverage of macrophytes, abundance classes (%) of Ephemeroptera, Plecoptera, and Trichoptera (in spring), and abundance of Trichoptera (in spring); however, correlation coefficients were negative for IPS in autumn. Furthermore, TDI in autumn and total abundance of benthic invertebrates (in autumn) had the highest positive correlation with axis 2. The abundance of Bivalvia, number of Ephemeroptera, Plecoptera, Trichoptera, Coleoptera, Bivalvia, and Odonata (EPTCBO) taxa, and number of macrophyte species also correlated positively with axis 2. These metrics indicated higher abundances in 2018 at both the upstream and downstream sampling sites, forming a separate group from other samples, thus showing interannual differences in communities of biological quality elements. The upstream site also differed from other samples in 2017 and showed a negative correlation with axis 2 (Figure 3).  Spring: Abund_S, total abundance in spring; EPTpaclS, abundance classes (%) of Ephemeroptera, Plecoptera, and Trichoptera in spring; EPT/TT_S, number of Ephemeroptera, Plecoptera, and Trichoptera taxa to the total number of taxa in spring; EPTCBOTS, number of Ephemeroptera, Plecoptera, Trichoptera, Coleoptera, Bivalvia, and Odonata taxa in spring; BivAb_S, abundance of Bivalvia in spring; TrichAbS, abundance of Trichoptera in spring. Autumn: Abund_A, total abundance in autumn; EPTpaclA, abundance classes (%) of Ephemeroptera, Plecoptera, and Trichoptera in autumn; EPT/TT_A, number of Ephemeroptera, Plecoptera, and Trichoptera taxa to the total number of taxa in autumn; EPTCBOTA, number of Ephemeroptera, Plecoptera, Trichoptera, Coleoptera, Bivalvia, and Odonata taxa in autumn; BivAb_A, abundance of Bivalvia in autumn; TrichAbA, abundance of Trichoptera in autumn.

Indicator Species Analysis of the Rusinupe
Two indicator species showed significant interannual variability (p < 0.05) for 2017 (one from the Diptera genus and one Plecoptera species). Upstream and downstream site differences showed nine indicator species (p < 0.05): one Gastropoda and one Ephemeroptera species for the upstream site, and two Gastropoda, one Isopoda, one Zygoptera, and three Trichoptera species for the downstream site. Rheophile indicator species were characteristic of the downstream site. Seven indicator species showed seasonal differences (p < 0.05): two significant indicator species (one Ephemeroptera and one Trichoptera species) were found for spring; five taxa (one Zygoptera and two Diptera families, in addition to two Trichoptera species) were found for autumn (Table 4). According to the PCA of the metrics in the Rusinupe, axis 1 of selected benthic invertebrate PCA explained 56.96% of the total data variance, whereas axis 2 explained 22.72%. High correlation coefficient values with axis 1 were found for eight metrics, from which the majority were calculated from the abundance of Bivalvia. However, high correlation coefficient values with axis 2 were found for only one metric, calculated from the abundance of Gastropoda in spring.
The first two axes of phytobenthos, macrophyte, and selected benthic invertebrate metric PCA explained 58.6% of the total data variance (40.9% and 17.7% for axes 1 and 2, respectively) in the Rusinupe. Significant correlations with PCA ordination axes were found for all BQE group metrics. The highest positive correlation coefficients with axis 1 were found for abundance of Bivalvia (number of individuals, % from total abundance and abundance classes (%)), abundance of burrowing/boring organisms (%) in spring and autumn, and active filter feeders (%) in spring and autumn. With axis 2, the highest negative correlation coefficient was for abundance of Gastropoda in spring (Figure 4).
The first two axes of phytobenthos, macrophyte, and selected benthic invertebrate metric PCA explained 58.6% of the total data variance (40.9% and 17.7% for axes 1 and 2, respectively) in the Rusinupe. Significant correlations with PCA ordination axes were found for all BQE group metrics. The highest positive correlation coefficients with axis 1 were found for abundance of Bivalvia (number of individuals, % from total abundance and abundance classes (%)), abundance of burrowing/boring organisms (%) in spring and autumn, and active filter feeders (%) in spring and autumn. With axis 2, the highest negative correlation coefficient was for abundance of Gastropoda in spring (Figure 4).

Water Chemistry
The major advantages of forest fertilization include a reduction in greenhouse gas emissions and economic benefits from the additional wood volume increment [53]. However, the possible adverse effects of this practice on the aquatic environment are not fully understood, particularly in terms of fertilization rates and experimental scale, as high variability of enrichment level can be found at a scale of less than 1 km. The variability in the results of nutrient enrichment experiments can be explained by the sedimentary environment, benthic macrofauna, and macrophyte community [54]. In general, our short-term monitoring results did not indicate an apparent decrease of ecological quality according to chemical parameters in the studied streams, thus supporting the conclusions of an earlier study [20]. In the Age tributary, a significant decrease in DOC and an increase in P TOT concentration were observed one year after applying wood ash, whereas the K concentration after applying wood ash was not statistically significant. It should be noted that aquatic chemistry parameters, including P TOT , K, DOC, and N TOT , display strong seasonal changes due to differences in nutrient input sources and pathways to streams, in addition to nutrient consumption by algae and macrophytes [55]. Although P TOT concentration during the first summer was higher than during the second summer or a year before treatment, due to the short study time and absence of replicate sites, it was not possible to exclude the impact of interannual variability. Results from a similar study [6] conducted in a Swedish peatland showed elevated pH, P TOT , and K concentrations for up to two years after applying wood ash at concentrations (3.1 t ha −1 ) similar to those in the present study. A study [20] on the long-term effects of wood ash application in Finland showed increased K concentration but no effect on P TOT , indicating the slow release of P TOT from wood ash [56]. Similarly, the fertilization of organic soil with NH 4 NO 3 did not affect N TOT concentration, although it has been reported that nitrate leaching occurred from acidic, sandy podosol [57]. In the present study, even after applying wood ash, concentrations of P TOT and N TOT in the Age were low and corresponded to high ecological quality both before and after the treatment. In the Rusinupe, N TOT showed good ecological quality, whereas total P TOT showed poor quality at both the control and treatment sites. We assume that the increased phosphorus concentration in the Rusinupe is caused by the loading of P TOT from drained forest soils in the catchment, as no other significant pollution sources were detected.

Biological Quality Elements
The species richness of diatoms observed in the Age and the Rusinupe would be considered relatively high and comparable with the taxa numbers of small streams in Sweden [5]. Of the species characteristic of forest streams, Gomphonema parvulum was common in the Age, Navicula cryptocephala in the Rusinupe, and Meridion circulare in both. Meridion circulare prefers alkaline or neutral pH values and tolerates slightly elevated concentrations of organically bound nitrogen [58]. Gomphonema parvulum is a eutraphenic species that tolerates low oxygen concentrations (>30%). Navicula cryptocephala requires moderate oxygen saturation (>50%), and M. circulare prefers oxygen saturation above 75%, as both are hypereutraphent species [59]. Their ecological requirements and occurrence reveal these species to be generalists [60]. In addition, our results of the diatom indices failed to detect any changes in the ecological quality of the studied streams, supporting the reported IPS values in Finland and Sweden [5,58]. The IPS values obtained at the streams examined in our study were comparable to those of boreal non-urban streams (mean 14.0), as reported by Teittinen et al. [58]. Although we did not estimate the changes in the biomass of primary producers following treatments with wood ash and NH 4 NO 3 , Hensley et al. [26] reported no impact of fertilization on the production of abundant algal and plant biomass.
The assessment of the effects of nutrients on the status of macrophytes in running waters is confounded by the synergistic effects of the other environmental and biotic variables that affect their growth, e.g., light conditions, substrate type, or stream velocity [61]. The most frequent macrophyte species observed in studied streams are tolerant to habitat degradation and other types of impacts, such as organic pollution [62]. In addition, the free-floating macrophyte species are limited by stream velocity. They reach their highest abundances in slow-flowing streams with sandy and soft, silty substrate similar to the treatment site in the Rusinupe.
Benthic invertebrates showed high seasonal and interannual variation in both species composition and abundance, in addition to in quality indices, thus supporting the results carried out in large rivers [63,64]. Moreover, the differences in macroinvertebrate assemblages between the treatment and control sites in the Age and Rusinupe are most likely caused by the natural variability of macroinvertebrate communities. Nevertheless, our results on macroinvertebrate abundance concur with those of a 3-year study in Germany, showing a distinct variation in seasonality. Regarding the macroinvertebrates, long-term studies are recommended to distinguish clear patterns in quantitative data [65]. Another cause of the macroinvertebrate data variability may be the beaver activity we observed at both studied streams in 2018. Beavers are reported to alter the macroinvertebrate communities, and the species richness of ecologically sensitive taxa, such as mayflies, stoneflies, and caddisflies, tends to decline [66,67]. Furthermore, the ISA results of benthic invertebrate communities support the statements mentioned previously. The herbivorous Lepidoptera indicator species Cataclysta lemntata (Table 3), which feeds on free-floating lemnids [68], may be evidence of beaver activity in 2018. The presence of seasonal indicator species in both streams could be suitably explained by the different life histories of aquatic invertebrates [69]. In the Age, the most significant differences were found within sampling year and sampling season rather than within sampling sites. In the Rusinupe, the most significant differences were found within sampling sites and sampling seasons rather than within sampling years. Such differences may be explained by the varying size of streams and the larger impact of local environmental factors [64] on benthic invertebrate communities in a small stream such as that of the Rusinupe. Studied biological quality element metrics showed high seasonal, interannual, and spatial data variability and we could not describe clear negative or positive patterns of diatom, macrophyte, and benthic invertebrate metric results.
Streams are open ecosystems, strongly connected to their riparian zone, and the aquatic biota are connected to the food resources and nutrients received from the terrestrial environment. Separation of the natural environmental factor impact from an impact of anthropogenic factors on the stream ecosystems remains an open question. A study of plant communities [70] emphasized that climate change induced significant alterations in the hydrological characteristics of lowland streams in temperate regions. Moreover, climate change effects on the structural and functional properties of riparian ecosystems remain to be more fully elucidated.

Conclusions
In conclusion, no negative effects of forest fertilization on the studied streams were observed, though we are aware that a complex set of different factors, e.g., site variability, beaver activity, and seasonal or interannual variations, may affect the ecological quality of the lotic ecosystems. Therefore, we suggest that longer-term monitoring of reference conditions is needed to distinguish the impact of forest fertilization from the natural variability of physico-chemical parameters and biological quality elements. Additionally, site-specific hydrology, measurements of forest soil chemistry, and shallow groundwater are required for better understanding of the processes and effects caused by forest fertilizers.

Data Availability Statement:
The data used in this study are available on request from the corresponding author.