Rapid Changes in Ground Vegetation of Mature Boreal Forests—An Analysis of Swedish National Forest Inventory Data

: The boreal forest ﬂoor vegetation is critical for ecosystem functioning and an important part of forest biodiversity. Given the ongoing global change, knowledge on broad-scale changes in the composition and abundance of different plant species and species groups is hence important for both forest conservation and management. Here, we analyse permanent plot data from the National Forest Inventory (NFI) on changes in the vegetation over a 10-year period in four regions of Sweden. To limit the direct and relatively well-known effects of forest management and associated succession, we only included mature forest stands not inﬂuenced by forestry during the 10 years between inventories, and focused on vegetation change mainly related to other factors. Results show strong decrease among many species and species groups. This includes dominant species such as Vaccinimum myrtillus and Deschampsia ﬂexuosa as well as several forest herbs. The only species increasing are some mosses in the southern regions. Our data do not allow for a causal interpretation of the observed patterns. However, the changes probably result from latent succession in combination with climate change and nitrogen deposition, and with time lags complicating the interpretation of their relative importance. Regardless of the cause, the observed changes are on a magnitude that suggest impacts on ecosystem functioning and hence highlight the need for more experimental work.


Introduction
The boreal forest floor vegetation is critical for ecosystem functioning and constitutes an important part of forest biodiversity. It affects soil properties, nutrient cycling and stand succession and hence overall stand structure and development [1,2]. Given the ongoing global change, knowledge on broad-scale changes in the composition and abundance of different plant species and species groups is hence important for both forest conservation and management.
Plants react and change their abundance in relation to both abiotic and biotic factors at various spatial and temporal scales. The most prominent changes occur during the succession after major disturbance events such as fire and forest harvest. These changes are well-studied in relation to vegetation succession following both natural disturbances and forest management [3][4][5][6], with plant species following relatively predictable change from shade intolerant to more shade tolerant species. In addition to succession, vegetation change is also governed by climate as well as other biotic and abiotic factors, and in the boreal forest is particularly related to nitrogen availability [7][8][9].
Fennoscandian boreal forests cover about 60 MHa [10] and are subjected to intensive sustained-yield forestry mainly based on clearcut harvesting systems [11]. The region is experiencing significant climate change with predictions for 2100 of temperature increase from 2-4 degrees and 10-20% increased precipitation (RCP4.5 [12]). These forests are naturally nitrogen-limited [13] and hence sensitive to elevated anthropogenic nitrogen deposition [7,14]. The basal area of the forests across Sweden has increased with on average 20% during the past three decades, with the largest increase in northern Sweden [15]. All these factors may influence forest-floor vegetation, and it is a challenge to disentangle the relative contribution of the direct effects of forest management from the indirect effects of climate and pollution [16]. Yet, this is critical to devise management strategies that protect both species richness (biodiversity) and overall ecosystem function.
National Forest Inventory (NFI) data encompass a large number of sample plots from the entire forest landscape and as such may provide an opportunity to better understand how different factors drive changes in forest vegetation [17]. Covering large geographical regions and the full range of stand ages, such data allow analysing vegetation change while controlling for specific factors such as stand age and forestry interventions. Hence, NFI data may increase our knowledge on the relative importance of different forest management operations and other factors influencing forest-floor vegetation. Although vegetation change is influenced by numerous interacting factors in a complex manner, NFI data represent an important empirical source for identifying general trends in vegetation change and for generating hypotheses on factors influencing these trends.
In the present study we utilize data from the Swedish NFI on the abundance (percentage cover) of 25 species (vascular plants, mosses, lichens) and 6 major species groups. We analyse the changes occurring during a 10-year period in mature forest stands that did not experience harvesting (final felling or thinning) during the study period. This limits the direct role of succession after harvesting and highlights the importance of other factors. Specifically, we aim to: • Identify if directional changes in percentage cover of species and species groups have occurred in old/mature forests that have not been subject to forest management during a 10-year period. • If such changes are present, to analyse if the direction of change varies among four regions in Sweden, stretching from hemiboreal to northern boreal regions.

•
To explore factors that might have contributed to the identified changes.

Study Area
The study area comprises all of Sweden except the southernmost part (temperate region) and follows the regional delineation used in the Swedish NFI ( Figure 1). About two-thirds of the area is covered by forests (based on the FAO, Food and Agriculture Organization of the United Nations, forest definition; [18]) with a total area of 28 million hectares [15]. Industrial forestry has occurred on most of the area since the early 1800s, with large-scale implementation of clearcut forestry from the middle of the 1900s [19].

NFI Data
The Swedish NFI includes a network of systematically located sampled plots where many variables, including the cover of species in the ground vegetation, are monitored [17]. The NFI design is based on clusters of sample plots in square-formed tracts (4-8 plots per tract depending on region). The plots are evenly spaced and located around the tract perimeter with the length of tract side varying from 300 to 1200 m among regions. For our analysis, we used data on cover of vascular plants, mosses and lichens in 100 m 2 permanent (with repeated inventories) circular plots during one full inventory cycle with the first and second inventory year separated by 10 years. The first year of inventory ranged from 1995 to 2003. The cover estimates include surveyor variation, but this variation is considered relatively small compared to random errors [20]. Hence, we base all our analyses on a pair-wise comparison of plots with a 10-year time interval between inventories. We only included plots on productive forestland (annual growth > 1 m 3 ha −1 yr −1 ). As auxiliary information we extracted data on thinnings taking place within 10 years prior to the first inventory (0/1), stand productivity (m 3 ha −1 yr −1 ), change in basal area (m 2 ha −1 ) over the 10-year period, change in canopy cover (%) over the 10-year period and stand maturity class at the first inventory (4 classes of old/mature forest; Table 1). The stand maturity classes included in the study have a minimum requirement of dominant and subdominant trees with a diameter at breast height (DBH) larger than 20 cm and are based on stand age in relation to lowest allowed and lowest recommended harvest age. The Swedish NFI includes inventory of 70 forest-floor species, or combinations of species not identified individually in the field [21]. Several of these species are relatively rare. For a species to be included in our analysis it had to occur in at least 20 sample plots in a region and be present in at least three of four regions. Although fulfilling the inclusion criteria, three species typical for early successional stages were excluded (i.e., Epilobium angustifolium, Rubus ideaus and Geranium sylvaticum), since they are assumed to be in a transient successional phase also in mature forests. With the applied selection criteria, our analysis is based on 25 species and more than 12,500 individual species by region observations across 1300 study plots in each of the two inventory years.

Statistical Analysis
The change in abundance was analysed by calculating a 95% confidence interval around the mean absolute change in cover between the two inventories for individual species and six major species groups (graminoids, herbs, dwarf shrubs, spore plants, lichens and mosses) within individual regions. For the species group analysis, we included all species surveyed by the NFI. Although the NFI plots are located in clusters (tracts), we assume that observed changes are independent among plots. The distance between plots is several hundreds of meters and even for the most frequent species, on average only two to three plots per tract are included in the analysis. This assumption was also made in several recent studies using vegetation data from the NFI (e.g., [4,6]).
The influence of the four stand variables on the change in abundance was explored by linear regression, using the function "lm" in the statistical software R [22]. The analysis was exploratory since there are a very large number of potential models that can be developed and the purpose was not to predict change per se. Rather, we wanted to exclude factors relating primarily to successional aspects and focus on change that is more likely due to other abiotic factors. Hence, we considered region as a "block" in the analysis, and each of the four stand variables was tested individually for a change in percentage cover for each species (see Appendix A, Table A1). For continuous variables (stand productivity and change in basal area), we only report the degree of explanation (adjusted R 2 ) when the factor significantly influenced percentage cover change. For the categorical factors (thinning prior to first inventory and stand maturity class) we provide for each class the percentage cover change for individual species when these factors were significant.
To explore the responses of individual species we selected a set of Ecological Indicator Values (EIV) [23] that directly link to the main factors influencing the abundance of vascular plants in boreal forests. We compared six EIVs (heat requirement, continentality, light, moisture, soil reaction (pH) and nitrogen) for species showing significant change with the EIVs for the other plant species monitored by the NFI by calculating mean EIVs for each of the groups. If the mean EIV differed by more than one unit, we consider this as an indication that the factor influenced the observed change (see Appendix , Table A2).
An NMDS ordination was performed for the six major species groups to examine the effects of inventory time and region on percentage cover. The functions "metaMDS" and "envfit" in R-package "vegan" [24] were used for the analysis, with a solution converging for a 4-dimensional solution based on a Bray-Curtis dissimilarity matrix (final stress 0.0618). Significance test of inventory time (first or second inventory) and region was based on 1000 permutations. For the full data set of all species, no converging solution was possible to obtain. Further details on ordination analysis are given in Appendix .

Stand Structure
On average, 17.5% of the 1300 sample plots was thinned during the 10-year period prior to the first inventory, ranging from 9.4% in region 2 to 25.5% in region 4. The basal area increased as expected, with about 10% across the regions (Table 2). Canopy cover estimates were only available for a limited number of plots (179) and indicate an increase of about 12% across the regions (Table 2).

Change in Species Abundance
We found significant changes in the abundance of many species in several regions during the 10-year period. Out of 93 analysed species-region combinations, 37 showed a significant decrease in cover, whereas only 4 combinations (all mosses) increased (Table 3), that is, about 44% of the species-region combinations showed a significant change in cover. Albeit the absolute change was relatively small, the relative change was high. For species that decreased, the average relative decrease was 34%, while those that increased changed on average with 22%. The number of species showing significant change ranged from 7 to 11 decreasing species and from 0 to 2 increasing species per region. Table 3. Species showing significant changes (p < 0.05) in cover in individual regions (see Figure 1). Initial cover denotes the average cover for the first inventory year while the absolute change and the relative change refers to the average change occurring after 10 years. N refers to the number of plots included in the analysis. The 95% confidence interval for the absolute change is defined by its two limits: the Lower Confidence Limit (LCL) and the Upper Confidence Limit (UCL).

Factors Correlating with Change in Cover
The explored linear regression models (including stand productivity, presence of thinning during 10 years prior to the first inventory, change in basal area during the 10year period and stand maturity class) were in most cases nonsignificant (see Appendix A, Table A1). However, some of the independent variables influenced cover of the species: Stand productivity-Reindeer lichens (Cladonia grp. Cladina) and other ground lichens tended to decrease their cover less on plots with higher forest productivity, while broadleaved grasses decreased more at high productivity plots. However, the models were very weak and explained 1% or less of the change in cover.
Change in basal area-Anemone nemorosa and broad-leaved grasses decreased more at sites with a higher increase in basal area. However, these models only explained about 2% of the change in cover.
Thinning before the first period-Melampyrum spp., Vaccinium myrtillus and Calluna vulgaris decreased in unthinned plot and increased in thinned plots. By contrast, Carex globularis and Pleurozium schreberi were negatively influenced by thinning (Table 4). Table 4. Relative cover change for species where thinning occurred during 10 years prior to first inventory significantly (p < 0.05) explained the observed variation. Sample size in parenthesis. Stand maturity class-The change in cover of some species (Oxalis acetosella, Anemone nemorosa, Vaccinium myrtillus, and Cladonia grp. Cladina) maintained their cover better in the younger maturity classes, except for Cladonia grp. Cladina, which decreased most in the youngest class. Lycopodiaceae showed varying cover change in the different maturity classes (Table 5). Table 5. Relative change in cover of species where maturity class (stand age proxy) was a significant (p < 0.05) factor explaining the observed variation. The sample size in given in parenthesis. For definition of maturity classes see Table 1.

Species 33 34 41 42
Oxalis acetosella The mean value of three EIVs for the vascular plants showing decreasing cover, differed with more than one unit compared to the other species monitored in the NFI plots (Appendix , Table A2); Heat EIV was 1.3 units lower, Soil pH EIV was 1.8 units lower and Nitrogen EIV was 2 units lower. Hence, the decreasing species occurred in colder climates and less productive forests compared to the other forest species.

Change in Species Groups
When combining the species into six functional groups, most of the groups showed a significant decrease in their cover. Only mosses in region 3 increased (relative increase of 9%). The change was not significant for seven of the 24 species group-region combinations ( Table 6). The average relative cover change for decreasing species groups was 23%, ranging from 4% to 44%. Herbs significantly decreased in all regions, while graminoids and lichens decreased in three regions. Table 6. Change in cover of six species groups in four regions (see Figure 1). Initial cover denotes the average cover for the first inventory year, while the absolute change and the relative change refer to the average change occurring after 10 years. N refers to the number of plots included in the analysis. The 95% confidence interval for the absolute change is defined by its two limits: the Lower Confidence Limit (LCL) and the Upper Confidence Limit (UCL). NS denotes nonsignificant (p > 0.05) change. The NMDS ordination of the six groups mirrors the general pattern of significant change in most species groups ( Figure 2). The ordination highlights that the composition differs among regions and that it is change over time (both factors significant, p < 0.01). However, only a minor fraction of the total variation in species composition was explained. The ordination suggests that the composition becomes more similar and moves towards the composition in region 1, with a stronger relative dominance of mosses and dwarf shrubs. Additional details on the ordination results are presented in Appendix .

Discussion
Our analysis of the change in cover of ground vegetation in mature boreal forests shows a rapid and significant decrease for many of the studied species and most species groups. We recognize that our data do not allow for a formal cause-and-effect analysis of the reasons behind these changes, but the results highlight the need to explore potential factors that contribute to our observations. Our data are large and for most species and species groups include several hundreds of randomly distributed plots throughout boreal and hemiboreal Sweden, adding strong statistical support for the observed patterns.
Previous studies using NFI data from Sweden and Finland have demonstrated clear effects of forest management on forest-floor vegetation [4,6,25,26]. Their results can mostly be attributed to the direct effects of forest management (clearcut harvesting and thinning) and the subsequent succession in forest-floor vegetation. Changes in light, nutrient and moisture availability are obvious after forestry interventions, and vegetation dynamics reflect these changes in abiotic conditions. Two recent studies [25,26] acknowledged that climate change may also play a role in observed changes, but the impact is likely minor compared to the effects of the conducted forest management operations. In our case, we have as far as possible removed the direct effects of management activities by focusing on mature forest stands that have not been subjected to major anthropogenic disturbances.

Succession
Although the study plots have not been recently subject to major management interventions, successional transitions may still potentially explain at least part of the observed changes. Several studies from Fennoscandia have described the pattern of vegetation dynamics after disturbance by using comparable data. For example, in a study from northern Finland [27], succession after disturbance was studied in both managed (after clearcutting) and natural forests (after forest fire), ranging from recently disturbed sites to overmature stands (>100 years). In the managed forest stands, Vaccinium myrtillus increased in cover from mature to overmature forests. Also in Swedish forests, the cover of V. myrtillus increased up to a forest age of about 120 years [28]. This was not the case in our study. Here the cover of V. myrtillus decreased in the two northern regions and remained rather stable in the two southern regions. Hence, our results are more in line with recent studies in Sweden [5,6] showing a general decrease of V. myrtillus and other dwarf shrubs, which could be attributed to increasing stand density across Swedish forests (see [15]). These contrasting observations add to a somewhat puzzling variation in the abundance of this keystone species in Fennoscandian forests [4]. Probably there is an interaction between stand density and stand age in these studies. The importance of stand density is supported by our data, since V. myrtillus decreased less in stands thinned prior to the first inventory period and in younger stand maturity-classes.
It is well-established that light availability, and hence canopy cover, in combination with precipitation strongly affect lichen growth [29,30]. In a study comparing lichen growth rate across sites with different canopy cover, it was shown that canopy cover over 60% reduced growth of two common terricolous lichens (Cladonia stellaris and Cetraria islandica) [31]. This suggests that increasing canopy cover could result in decreased abundance of ground lichens also in our data set. During the studied 10-year period, the canopy cover increased with on average about 6%, from just below to just above the 60% threshold. It should also be noted that the basal area in mature forests in our dataset is higher than in the cited study [31].
However, if our observed changes would be mainly associated with a successional trajectory from younger to older forest, with increasing basal area and canopy cover, we would expect to see consistent effects of the changes in basal area during the study period as well effects of thinning events prior to the 10-year inventory period. Our analyses clearly indicate that this is not the case. We found only minor effects of the change in basal for a few species, that is, for Anemone nemorosa, broad-leaved grasses and for three species that increased in thinned stands and decreased in unthinned stands (Table 4). This suggests that the observed changes in cover of major plant groups are not only driven by successional change and that other factors may contribute to the changes in vegetation.

Climate
A significant increase in average temperature and precipitation was observed throughout Sweden from the early 1990s. Compared to the reference period 1961-1990, the average temperature in Sweden was 1.0 degrees higher during the inventory period 1995-2014, with a 9% increase in precipitation [32]. The long-term trend with a warmer climate is clear, and the response of the vegetation may well be associated with this trend. Given that our observed changes cannot be fully explained by succession, climate change should at least be considered a contributing factor. This is supported by the EIV analysis, since the decreasing species tended to be more associated with colder environments than the other vascular plants monitored by the NFI.
A warmer climate extends the vegetation period and may shift the competitive balance among species, with some winners and some losers. It is well-established that soil warming has significant effects on ground vegetation in tundra habitats [33]. However, it is likely that closed canopy forests mediate the effects of soil temperature on ground vegetation [34]. In a long-term experimental study of soil warming, Hedwall et al. [35] found a strong initial vegetation response in young forests, but the long-term effects were less clear, except that moss cover increased as the forest became older. The soil temperature during winter is strongly connected to the insulating effect of snow cover [36]. Milder winters in the south will give a less consistent snow cover, which may increase soil temperature. Hence, we hypothesize that increased temperature may contribute to the observed increase in moss cover in southern Sweden. By contrast, a decrease in snow cover in northern regions may result in stronger ground frost. This has been shown to have dramatic effects on ground vegetation, with significant decrease in the cover of understory vegetation and mosses [37], in accordance with our observations. An increase in precipitation should favour species and species groups typical of more moist conditions. However, our results contradict this since plants typical for moist forests (Equisetum sylvaticum, Rubus chamaemorus, Carex globularis and Sphagnum spp.) also decreased in cover. Neither do the EIVs provide any indication that species associated with moist habitats were favoured.

Nitrogen Deposition
Deposition of atmospheric nitrogen is a factor known to have significant impact on boreal forest ground vegetation [38]. There is a clear gradient through Sweden with higher deposition in the southwest and lower towards the north [39,40]. Nitrogen deposition has decreased since the 1980s across Europe [41]. However, the change is not as evident in Sweden, besides an indicative trend of decreasing nitrogen deposition in the southwestern part (i.e., including parts of region 4) since around year 2000, and for long time series in northern regions [40]. Yet, the EIV analysis suggests that species associated with more nutrient-poor environments (low soil pH and less nitrogen demanding) are overrepresented among species showing decreasing cover.
In experimental studies it has been shown that narrow-leaved grasses (dominated by Deschampsia flexuosa) and herbs increase with increasing nitrogen deposition while dwarf shrubs and lichens decrease [14,42]. In comparison with our results, the effect of nitrogen corresponds to the decline in dwarf shrubs and lichens, but not for narrow-leaved grasses and herbs, since all these species groups decline in mature Swedish forests. It should be noted, however, that even if experimental studies often show clear effects on ground vegetation, the effect of the overall atmospheric deposition of nitrogen is less evident, mainly because experimental doses applied are commonly much higher than background deposition [38]. This complicates the interpretation of our results.

Time Lags in Responses
Most of the studied species and species groups are clonal and long-lived. Hence, although abiotic conditions change, associated changes in cover and occurrence of plants likely include time lags, that is, integrating effects over a longer time period than the studied 10-year period. For instance, major distributional and abundance changes in Finnish ground vegetation are predicted up to mid-late 21st century due to increased temperature [9], and hence highlight trends over longer time periods. When evaluating potential factors behind the observed changes, this implies that increasing stand age and changes in climate and nitrogen deposition occur in parallel time scales, which might require several decades before fully manifested.

Conclusions
We have shown large changes in the ground vegetation in mature Swedish forests over a relatively short time period (10 years). These changes may imply ecosystem-level impacts on the function of boreal forests [43]. At this stage, we cannot attribute the observed changes to a single factor. The changes likely result from a combination of succession, climate change and nitrogen deposition, and where time lags complicate the interpretation. This confirms the conclusion from the review by [16] that it is necessary to consider both previous land management and functional properties of plant communities when analysing and modelling effects of the ongoing global change. Table A1. Effects of stand productivity, basal area change, presence of thinning 10 years prior to first inventory period and stand Maturity class (see Table 1) on the change in cover during the studied 10-year period. Values are p-values from linear models with regions treated as a block in the analysis. Bold numbers indicate a significant effect. N denotes number of plots where the species occurred.

Appendix B
To further explore patterns among the vascular plants showing significantly decreasing cover we utilized the ecological indicator values (EIV) provided by [23]. Although available for a large range of aspects, we limit this comparison to EIVs directly linked to factors well-known to influence boreal forest plants, that is, heat requirement, continentality, light, moisture, soil reaction (pH) and nitrogen. For these six EIVs, an ordinal scale is provided by [23]. Although the scales are ordinal, we have calculated the mean EIV values for the species showing a significant decrease in cover in our study and compared these with the mean values for all other vascular plants included in the Swedish NFI. However, note that broad-leaved grasses cannot be included since the group includes several species with contrasting EIVs.

•
Heat indicator-range from 1 (high alpine) to 13 (climate cultivation zone 1 in southernmost Sweden) • Continentality indicator-range from 1 (hyperoceanic) to 9 (hypercontinental) • Light indicator-range from 1 (deep shade) to 7 (always in full sun) • Moisture-range from 1 (very dry) to 12 (deep permanent water) • Soil pH-range from 1 (strongly acidic, pH < 4.5) to 8 (alkaline, pH > 7.5) • Nitrogen-range from 1 (very N-poor) to 9 (mostly on artificially N-enriched soils)    Figure A1. NMDS ordination (dimensions 3 and 4) of the six different species groups where markers provide the average "sites scores" for the four different regions at the two inventories. Arrows denote the direction of change. Note that, to increase the readability of the figure, the placement of the species group scores is not representing their actual scores since these are outside the current axes. However, their relative position indicates where they occur in the ordination. The difference between regions remains clear also along the NMDS 3 and 4 axes, while the change between periods is negligible.