Introduction 

Norway spruce (Picea abies [L.] Karst.) is one of most important tree species in European temperate and boreal forests and has a long tradition of cultivation ([11], [6]). Besides its ecological importance, Norway spruce is a backbone species for the timber industry in many European countries ([18]). Similar to other economically valuable tree species ([41]), Norway spruce has been favoured and planted beyond its ecological niche. However, it is currently under severe pressure due to global warming and intensified forest disturbances of various forms. Outside its natural distribution in lowlands, this tree species is particularly susceptible to heat and drought spells because of its shallow root system ([6]) and to forest pest infestations such as Eurasian spruce bark beetle (Ips typographus L.), which is arguably one of the most damaging agents in temperate and boreal forest ecosystems ([30], [9], [26], [40]). Special attention should be given to climate change induced transitions towards alternative vegetation states and irreversible alterations in the structure and composition of mountain spruce forests in their natural, high-elevation distribution ranges ([23], [1]). Extensive outbreaks of Ips typographus at higher latitudes ([40]) are another reason for concern.

In the last decade, Slovenia has witnessed a series of back-to-back forest disturbance episodes on an unprecedented scale (an ice storm in 2014, widespread bark beetle outbreaks in 2015 and 2016, windthrows in 2017 and 2018). The largest ever recorded ice storm event in the region occurred in 2014 ([25]), causing massive tree damage (9.32 millions m³) in more than half of Slovenian forests, with a long-lasting effect on forest dynamics in the following years. Catastrophic storms can result in large-scale destruction of forests and cascading effects starting with pervasive infestations of secondary pests or diseases ([36]). While the immediate ice damage was more pronounced for broadleaf tree species, the ice storm was the primary trigger for a dramatic increase in bark beetle outbreaks throughout the country, particularly in Norway spruce stands. Spruce trees were either uprooted or snapped, or the crown was lightly to severely damaged ([7]). The increasing frequency of climate extremes combined with unfavourable forest structure has pushed spruce forests to the edge of their ecological space and has led to the initiation of large-scale forest conversions towards more natural forest vegetation ([18]), a process that has been gaining importance since the early 2000s ([34]). Stand-replacing forest disturbances have caused a transition from mostly local outbreaks to those at regional and national levels ([19]). Disturbance events have been accompanied by extensive post-disturbance management actions (extensive salvage logging and wood transportation) that additionally affected not only stocks of damaged trees but also site conditions (e.g., soils) and the development of young stands (e.g., tree regeneration). These impacts are especially challenging for close-to-nature forest management, which simulates natural processes operating on much smaller spatial scales ([29]).

Site, stand and landscape conditions as well as past land use in forests have a significant influence on their current status, sensitivity and resilience ([13], [31], [24]). Abiotic site conditions, such as a warmer climate or shallow calcareous soils in spruce forests, can contribute to the development of natural disturbances and negative impacts on these forests, of which bark beetle infestations and windthrows are among the most common in Central Europe ([38]). Analysing a decline in the growing stock of a tree species across a wide spectrum of site conditions can provide insight into the climatic, soil, topographic and forest stand factors under which a decline is more likely to occur ([3], [8], [26]). However, many approaches aimed at identifying such factors are based on plausible or intuitive narratives rather than on tested and data-driven concepts. Identified patterns of past spruce decline have generally suggested that certain site conditions (e.g., warmer sites with dry and/or shallow soils) make spruce forests more susceptible to damage resulting from large-scale disturbances. On the other hand, more favourable site conditions reflecting the natural occurrence of primary spruce forests can contribute to lower susceptibility. Ecology-based predictions and derived findings have potentially important forest management implications. However, causal connections show a high degree of complexity as numerous factors and their interactions are involved. Therefore, a complete understanding exceeds the scope of a single study. Given that natural disturbances in forests are characterized by high degree of stochasticity and that disturbance patterns are highly variable ([32]), it is of paramount importance to enrich our knowledge with real-world data-derived and hypothesis-driven results. Studying the dependence of currently prevailing tree species decline on ecological factors can provide a firm basis for improving forest management systems aimed at mitigating the effects of climate change and intensified natural disturbances, particularly for vulnerable forest types such as Norway spruce forests.

Here, in order to identify the main ecological drivers of the recent Norway spruce growing stock decline, we investigated the response of the relative decline along a suite of abiotic site factors (related to geology, climate, soil, terrain topography and forest stand characteristics) using a representative sample of forest compartments across Slovenia. The primary focus of this study was to describe how the effects of the selected predictors associated with climate, soils, topography and stand characteristics might depend on bedrock type and forest naturalness. Naturalness is defined as the degree to which the existing tree species composition reflects the potential natural composition of forest stands ([4], [5]). The relevance of both factors as potential regulators (either mitigating or stimulating a growing stock reduction in forest stands) is well established within the research agenda of forest ecology and management ([7], [3]). However, our understanding of ecological predictors controlling tree species growth and survival remains limited. Moreover, their expected impacts might be blurred by the occurrence of intense forest disturbances with various dependence patterns. A particular disturbance agent may depend more on one factor (e.g., bark beetle outbreaks strongly depend on climatic conditions and the associated physiological status of trees), while the incidence of other disturbances can be attributed to different forest site factors (e.g., stressful abiotic conditions may not be a necessary precondition for windthrows).

In relation to the two main predictors analysed, we predicted that (i) compartments on carbonate bedrock would exhibit a higher spruce decline compared to those on silicate bedrock due to their higher susceptibility to drought stress and (ii) compartments with lower forest naturalness would also exhibit a higher spruce decline given that the tree species composition in such compartments is less adapted to the site conditions.

  Materials and methods 

Study area

Slovenia is located in the transition zone between Central and South-eastern Europe. At the crossroads between the Mediterranean, Alpine and continental climatic regimes, it experiences large climate variation ([35]). Despite its relatively small size, it is comprised of the Alpine, Pre-Alpine, Dinaric, Pre-Dinaric, Sub-Pannonian and Sub-Mediterranean biogeographical regions with their own climatic and forest vegetation specificities. For background information regarding the prevailing natural forest type, see Kutnar et al. ([19]). Covering almost 60% of the land area, forests have played an important role throughout Slovenian history. Although the principles of close-to-nature, sustainable and multifunctional management have been implemented to a greater degree compared to some other European countries, anthropogenic alterations of forest stands have nevertheless been a common practice. Today’s forest status often deviates from that which would have developed without human engagement, with only 22% of Slovenian forests exhibiting preserved tree species composition ([9]). The planting of economically profitable but ecologically unsuitable tree species (poorly adapted to local site conditions) remains one of the main problems in Slovenian forestry. Such practices have resulted in large-scale forest transformations that are comparable to those reported for countries in Central Europe ([32], [18]).

Primary (natural) Norway spruce forests grow abundantly only in the Alps in the north of the country and to a much lesser extent in cold valleys and karst depressions in the Dinaric region in the southern part of the country ([19]). However, in past centuries Norway spruce was extensively planted outside of its natural distribution range because of its high adaptability, relatively low investment costs and valuable timber. Such management-created Norway spruce stands, particularly even-aged spruce monocultures at lower altitudes, are susceptible to damage caused by various abiotic and biotic factors, with minimal resistance against natural hazards ([13], [18]) such as bark beetle infestations ([27]).

Indicators of Norway spruce decline

The proportion of Norway spruce in the growing stock of Slovenian forests has declined significantly over the last 20 years, and this decline has been particularly evident in the last 10 years ([19]). In order to contextualize the spruce decline, we used the Annual Reports of the Slovenia Forest Service (SFS) for the period 2010-2018 ([43]). These reports include the amounts of salvaged trees (m³) in Slovenian forests due to different agents. We analysed selected indicators of spruce decline (temporal dynamics of the percentage of Norway spruce in the growing stock of Slovenian forests, the ratio between salvage and regular logging/felling of spruce) and intensified impacts of forest disturbances (volume of felled coniferous tree species due to different disturbance agents, total area of completely devastated forest stands without any tree canopy cover and a high percentage of bare soil). The proportion of Norway spruce in the total growing stock of Slovenian forests declined from 31.5% in 2010 to 30.5% in 2018 (Fig. 1a). The total salvaged volume of spruce (caused by all disturbance agents) in the period 2010-2018 represented 52.4% of the total harvest. However, in the years before the ice storm disturbance (2010-2013), this proportion was significantly lower (29.5%) compared to the period 2014-2018 after the ice storm (70.8%). Planned (regular) harvests of spruce were thus greatly reduced (Fig. 1b). Salvage logging of coniferous tree species during the period 2010-2018 (Fig. 1c) amounted to a total of 16,238,715 m³. Half (50.2%) of this volume was felled due to insects. The amount of salvage logging in Slovenia is most affected by the occurrence of natural disturbances and insect infestations, especially spruce bark beetle (Ips typographus) attacks. Wind was responsible for 21.8% and the massive ice storm event for 16.0% of salvaged volume (Fig. 1c, evident peak for this agent in 2014). The total volume of salvage logging in the period 2014-2018 (13,901,160 m³) was almost 6-times higher than that in the period 2010-2013 (2,337,560 m³). Large-scale disturbances caused a marked increase in completely devastated forest areas, with a total of 8194 ha in the period 2015-2018. Before the ice storm in 2014, these areas were significantly smaller (407 ha in total for the years 2010, 2011 and 2012 - Fig. 1d).

Fig. 1 - (a) Percentage of Norway spruce in the growing stock of Slovenian forests. (b) The ratio between salvage (sanitary) and regular logging/felling of Norway spruce. (c) Volume of felled coniferous tree species due to different disturbance agents. (d) Area of completely devastated forest stands (large canopy openings and cleared, deforested areas) caused by bark beetle outbreaks (data for 2013 and 2014 not available). All panels refer to the period 2010-2018. All data were extracted from Annual Reports provided by the Slovenia Forest Service ([43]).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Data acquisition and handling

Extensive data on more than 335,000 forest stands and more than 53,700 forest compartments were compiled from Slovenia Forest Service databases for the year 2010 ([44], [45]) and for the year 2018 ([46], [47]). A forest stand is a relatively homogeneous part of the forest in terms of tree species composition, developmental stage, age of the predominant trees, etc. and corresponds to a spatial unit of a typically minimal area of 0.5 ha that is important for forest inventory, planning and management. A forest compartment is a permanent management unit consisting of several different stands, with an average size of 22 ha ([19]). Compartment-level data were derived from stand data, by aggregating stand level data belonging to the same compartment.

From the total pool of forest compartments across Slovenia (more than 53,700 overall), we performed a data sub-setting procedure during which certain arbitrarily defined criteria were used. The analyses included compartments with a minimal area of 5.0 ha, with spruce growing stock in the year 2010 equal to or greater than 20.0 m³ ha^-1 and with a decrease in spruce growing stock (negative relative change) from 2010 to 2018. From the forest stand data, we calculated the average growing stock (m³ ha^-1) of Norway spruce for each compartment. For each compartment, the relative change (%) from 2010 to 2018 in spruce growing stock (GS, m³ ha^-1) was calculated with the formula: [(GS₂₀₁₈ - GS₂₀₁₀) / GS₂₀₁₀] × 100. This resulted in 17.618 compartments in total. This dataset was then divided into two subsets. The first subset (henceforth “subset 1”) included compartments with a relative decline in spruce GS equal to or greater than 10.0%. The second subset (henceforth “subset 2”) included compartments with a relative decline in spruce GS below 10.0%. The division between these two subsets was made in order to analyse their general differences in terms of ecological predictors and other characteristics. Both subsets were additionally reduced by excluding the compartments that were categorized as non-forested (which included fragmented landscape and land use types other than forests) and with an area larger than 50.0 ha. This brought us to the final selection of suitable compartments, i.e., a total of 6355 for subset 1 and a total of 4396 for subset 2, totalling 10.751 compartments. Only subset 1 (including compartments with a spruce growing stock decline > 10%) was subjected to core statistical analyses. The compartments in subset 1 were considered to be more reliable (sufficient drop in growing stock) as they reflect changes caused by sanitary logging and large-scale disturbances. In contrast, compartments with lower spruce decline may include the effects of regular logging interventions as a consequence of small-scale management (these effects were not the focus of this study). Subset 2 (including compartments with spruce growing stock decline < 10%) was used for comparison with subset 1. The total area of compartments was 132,760.5 ha in subset 1 and 94,426.9 ha in subset 2. The mean size of compartments was 20.9 ha in subset 1 and 21.5 ha in subset 2.

Compartment-level data included seven main explanatory variables related to geology, forest naturalness, climate, soil texture, terrain topography and forest stand characteristics. We used two categorical predictors. First, the SFS database contained information on the dominant parent material, with 38 different levels ([45], [47]). To simplify this large number of categories, we implemented an expert-based classification where each parent material was classified into one of the three broad categories for bedrock type: carbonate, mixed and silicate bedrock. This division was made based on selected properties, e.g., CaCO₃ content and the speed and type of bedrock weathering (physical or chemical). The group of carbonate geological substrate was mainly represented by limestones and dolomites. The group of mixed geological substrate included bedrock types such as flysch and carbonate quartz sandstones, whereas the most frequent substrates in the group of silicates were mica schists, quartz sandstones and tonalite. Geology largely defines the physical and chemical properties of forest soils.

Second, information on the degree of forest naturalness was chosen due to its management relevance and close connection with spruce planting outside its ecological niche. The change in tree species composition is a measure of the degree to which the existing forest vegetation deviates from the potential vegetation ([4]). This variable was described with four classes as defined in the SFS databases: preserved (potential natural tree species composition is modified 0-30%), changed (modified 31-70%), strongly changed (modified 71-90%) and altered (modified 91-100%). A higher degree of composition preservation implies higher forest naturalness ([5], [9]).

The other five predictors were numerical (quantitative) variables. Based on the geographical coordinates (latitude and longitude) of the centroid of each compartment, we extracted open-source data for potential annual evapotranspiration (PET, representing atmospheric water demand) and clay content (CLAY) in the soil (depth of 0-100 cm). These two were sourced from the Global Aridity Index and Potential Evapotranspiration Climate Database v2 ([39]) and SoilGrids database ([15]). Potential evapotranspiration reflects both temperature and precipitation conditions (higher PET values indicate warmer and drier sites). Soil clay content is normally associated with soil water storage capacity (a positive correlation between clay content and water storage capacity is expected), and clay content can be used as a proxy for soil moisture. However, the amount of dead water that is not available for the plants is high in clayey soils. Prior to PET, we also tested various climatic variables (e.g., mean annual temperature, amount of annual precipitation, mean temperature of the warmest quarter, precipitation of the driest quarter) derived from the WorldClim database ([12]). However, we ultimately elected PET because of its integrative nature. The R library “raster” ([16]) was used for data extraction from both databases.

Additional relevant predictors of spruce decline were available in the SFS database ([45], [47]). Data on stoniness and rockiness (%), slope aspect (°), slope steepness (°) and Norway spruce proportion (%) in the growing stock were considered. We averaged the percentages of stoniness and rockiness into one variable named “STOROC”. Stoniness and rockiness (defined as the coverage of stones with an average size < 30 cm and rocks with an average size > 30 cm on the forest surface) can be seen as a proxy for soil conditions, with more stony/rocky sites mainly exhibiting shallower soils. The heat load index (HLI) was calculated based on three parameters: latitude, slope aspect and slope steepness, according to equation no. 3 in McCune & Keon ([22]). A higher heat load index (HLI) is characteristic of steep, south-facing slopes, while lower HLI values suggest less sun-exposed, north-facing slopes. The proportion of Norway spruce in the growing stock of the entire compartment at the beginning of the period (year 2010) was also included. In our case, the compartment-level proportion of spruce generally depends on the spatial structure of stands within each compartment (e.g., few spruce-dominated stands within a matrix of prevailing broadleaved forests results in a low spruce proportion). Summarized information about predictor variables is given in Tab. 1. In addition, the mean altitude (m a.s.l.) and total area (ha) of each compartment were taken as supportive factors explaining differences between the two subsets and for interpreting potentially confounding effects in relation to the focal predictors.

Formulation of hypotheses

For each of the selected predictor variables, we set and tested a hypothesis regarding how the predictor was likely to affect spruce decline (Tab. 2). These assumptions were derived from a forest ecology-based understanding of site factors or stand characteristics within compartments. It is important to mention that the assumptions usually apply to bark beetle disturbance (eco-physiological aspect), while other types of disturbances (e.g., ice storm, windthrow) might act differently as their impacts are more dependent on the structural conditions of a forest stand. For example, specific microclimatic conditions (depending on the geographic position and local topography, bedrock and soil) can greatly increase spruce mortality due to acute or chronic drought stress, which increases the susceptibility of trees to bark beetle attacks ([26]).

Statistical analysis

Descriptive statistics for each predictor were calculated separately for subset 1 and subset 2. As the data did not fulfil parametric assumptions (normality and homoscedasticity), differences between the two subsets were tested with the non-parametric Mann-Whitney test and visualized for the comparison of data distribution.

Our dependent (response) variable was the relative decrease in spruce growing stock, expressed in % and calculated with the above listed formula. Using the data from subset 1, we fitted semiparametric models because we expected non-linear relations between the dependent variable and predictor variables. For this purpose, Generalized Additive Models (GAMs - [14]) were implemented using the “mgcv” R package ([42]). Predictors can be estimated in nonparametric fashion, and the particular relationship between the two variables can be assessed from a smoothed drawing line. To avoid a high degree of model complexity, each GAM contained one of the numerical predictors (PET, CLAY, STOROC, HLI or SPRUCE10) and one categorical predictor (bedrock type or forest naturalness). This yielded five models for each categorical predictor and 10 different GAMs in total. Fitting multiple GAMs allowed to explore in detail the response of spruce decline, whereas one “global” GAM with multiple predictors would only provide general trends. In the results section, we report the F-value as an indication of relationship strength between the response and predictor variables. Explained deviance (%) as an overall measure of the explanatory power of the model is also provided. In addition to relative spruce decline, we also modelled the response of absolute decline (expressed in m³) using the same method as that used for the relative decline.

Before running the statistical models, we checked the distribution of the response variable and all five quantitative predictors. Predictors were also tested for possible multicollinearity, but low correlation coefficients prevented this (see Fig. S1 in Supplementary material). When significant departures from normality were detected, appropriate data transformations were performed for both response and explanatory variables. For relative spruce decline (%) a square-root transformation was the appropriate choice and for absolute decline (m³) the log₁₀ transformation was used. To further improve the interpretability of model outputs, standardization of predictors was implemented, resulting in scaled and centred values with mean = 0 and standard deviation = 1. For each quantitative predictor, we tested differences between groups defined for bedrock type (Fig. S2 in Supplementary material) and forest naturalness (Fig. S3). The non-parametric Kruskal-Wallis test with Bonferroni correction for multiple comparisons was used due to data violating the assumptions of normality and homoscedasticity. The most commonly used threshold for p-values (i.e., 0.05) was set in all statistical procedures. Analyses were done in statistical computing software R version 4.1.1 ([28]).

  Results 

Comparison between subsets with low vs. high spruce decline

Forest compartments in subset 1 (including compartments with a relative spruce decline ≥ 10.0%) showed significantly higher potential evapotranspiration (PET), soil clay content (CLAY) and stoniness/rockiness (STOROC) compared to compartments in subset 2 (including compartments with a relative spruce decline < 10.0% - Fig. 2). In contrast, compartments in subset 1 were characterized by a significantly lower proportion of Norway spruce in year 2010 (SPRUCE10), altitude and size (area). The investigated subsets did not differ significantly in terms of head load index (HLI) values (Fig. 2).

Fig. 2 - Comparison of potential evapotranspiration (PET), soil clay content (CLAY), percentage of stoniness/rockiness (STOROC), heat load index (HLI), proportion of spruce in 2010 (SPRUCE10), altitude and area (compartment size) between the two subsets. Subset 1 includes forest compartments with spruce decline > 10% (n = 6355). Subset 2 includes forest compartments with spruce decline <= 10% (n = 4396). Differences were tested with the non-parametric Mann-Whitney test and p-values are given above each panel. The red dot with range within each violin plot represents the mean value ± standard deviation.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Both subsets differed with respect to bedrock type categories and classes of forest naturalness (Fig. S4 in Supplementary material). Almost 70% of all compartments in subset 1 were located on carbonate bedrock, while the share of compartments with carbonate bedrock was lower in subset 2 (59.1%). Silicate bedrock was more frequent in subset 2 (30.3% of all compartments) in comparison to subset 1 (22.5%). Compartments with mixed geology had a similar share in both subsets, i.e., 8.5% in subset 1 and 10.6% in subset 2. More than half of the compartments in subset 1 had preserved tree species composition, whereas this share was lower for subset 2 (40.0%). In other classes, the proportions were as follows: 35.1% (subset 1) and 39.0% (subset 2) for the class of changed composition, 10.5% and 16.6% for the class of strongly changed composition and 2.6% and 4.4% for the class of altered composition (Fig. S4 in Supplementary material).

Spruce decline per bedrock type categories and classes of forest naturalness

Considering only the focal data of subset 1, we investigated the general pattern of spruce decline for different categories of bedrock type and classes of forest naturalness. The overall mean value of relative spruce decline in subset 1 was 38.7%, ranging from 10.0% to 99.8%. There were statistically significant differences between bedrock type categories in terms of spruce decline. Spruce decline was significantly lower in compartments with silicate bedrock compared to those on carbonate bedrock (Fig. 3a). Compartments with mixed geology exhibited intermediate levels of spruce decline but were not significantly different from the category of carbonate bedrock type. Examining the pattern for classes of forest naturalness, we found that compartments with preserved tree species composition experienced significantly higher spruce decline compared to the other three classes (Fig. 3b).

Fig. 3 - Relative spruce decline (data from subset 1) for (a) three bedrock type categories (carbonate, mixed, silicate) and (b) classes of forest naturalness (preserved, changed, strongly changed, altered tree species composition). Different letters above each boxplot indicate statistically significant differences according to the Kruskal-Wallis test with Bonferroni correction for multiple comparisons.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

GAMs with bedrock type as a predictor

Spruce decline generally increased along the potential evapotranspiration (PET) gradient. The relationship was significant for all three bedrock type categories, i.e., carbonate, mixed and silicate bedrock (Fig. 4, first row; deviance explained = 7.94%). Forest compartments with higher clay content in the soil exhibited higher spruce decline and this pattern was detected for all three categories (Fig. 4, second row; deviance explained = 6.46%). The effect of stoniness/rockiness (STOROC) on spruce decline was significant only in case of carbonate bedrock, for which it showed an increase at the beginning of the gradient followed by a steady and slowly decreasing trend (Fig. 4, third row; deviance explained = 5.46%). Heat load index (HLI) proved significant in the category of mixed bedrock. An increasing trend of spruce decline along the HLI gradient was observed (Fig. 4, fourth row; deviance explained = 4.87%). Among all quantitative predictors, the strongest relationship between the response variable was found for the proportion of spruce in year 2010 (SPRUCE10 - Fig. 4, fifth row; deviance explained = 10.70%). Compartments with a lower share of spruce generally exhibited a higher relative decline of spruce growing stock and this was true for all three bedrock type categories. Across these categories, the dependent variable behaved in a similar fashion with some minor exceptions (e.g., STOROC and HLI). Overall, the strongest predictor was SPRUCE10, while HLI had the lowest explanatory power.

Fig. 4 - Estimated response of relative spruce decline (%) to selected ecological predictors modelled with Generalized Additive Models: (PET) potential evapotranspiration; (CLAY) clay content in the soil (depth 0-100 cm); (STOROC) stoniness/rockiness; (HLI) heat load index and (SPRUCE10) proportion of spruce in the growing stock in year 2010. Bedrock type was included in each model as a categorical predictor. The y-axis in each plot represents the partial effect of selected explanatory variable on the spruce decline. The green ribbon around each curve denotes the 95% confidence interval. The F-value along with accompanying statistical significance (significant results in bold; *** p < 0.001, * p < 0.05) of smooth terms is reported at the top of each panel.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Absolute spruce decline (log₁₀(m³) - Fig. S5 in Supplementary material) behaved quite differently compared to the relative decline described in Fig. 4. For example, while the relative decline evidently increased along the soil clay content gradient (Fig. 4), the absolute decline showed a less clear pattern (unimodal or even declining trend along the same gradient). Furthermore, the volume of felled spruce trees increased with the proportion of spruce in the year 2010, whereas the opposite trend was observed for relative decline. More importantly, additional analysis with respect to the absolute decline revealed that compartments with different geology showed various patterns of spruce decline, which can be seen in the shape of the curves produced by the GAMs. This further demonstrates that the effect of ecological predictors on spruce decline depends on bedrock type and suggests the presence of interactions between qualitative and quantitative predictors.

GAMs with the degree of forest naturalness

With respect to forest naturalness (FN) as a qualitative predictor, the patterns observed were similar to those for bedrock type. Higher potential evapotranspiration (PET) caused higher spruce decline, although for some classes of FN preservation (e.g., changed, strongly changed) this trend was less clear (Fig. 5, first row; deviance explained = 7.63%). We came to the same conclusion regarding the relationship between relative spruce decline and CLAY (i.e., increasing trend along the gradient) as it was significant for all four categories (Fig. 5, second row; deviance explained = 7.34%). On the other hand, the dependent variable responded quite differently among FN classes along the STOROC gradient as it showed a less clear trend for the preserved and changed class but an increasing trend for strongly changed and altered compartments (Fig. 5, third row; deviance explained = 6.41%). HLI proved to be a significant predictor only in the class of changed tree species composition, with a slight increase in spruce decline along the gradient (Fig. 5, fourth row; deviance explained = 4.43%). The proportion of spruce in the year 2010 was once again identified as the most powerful predictor (Fig. 5, fifth row; deviance explained = 9.42%). Its effect was significant for all four FN classes and a more or less evident decreasing trend can be recognized along the gradient. In addition, the dependency of the predictor’s effects on the degree of forest naturalness can be recognized, suggesting interactions between qualitative and quantitative predictors.

Fig. 5 - Estimated response of relative spruce decline (%) to selected ecological predictors modelled with Generalized Additive Models: (PET) potential evapotranspiration; (CLAY) clay content in the soil (depth 0-100 cm); (STOROC) stoniness/rockiness; (HLI) heat load index and (SPRUCE10) proportion of spruce in the growing stock in year 2010. The degree of forest naturalness (FN) was included in each model as a categorical predictor. The y-axis in each plot represents the partial effect of selected explanatory variable on the spruce decline. The green ribbon around each curve denotes the 95% confidence interval. The F-value along with accompanying statistical significance (significant results in bold; *** p < 0.001, ** p < 0.01, * p < 0.05) of smooth terms is reported at the top of each panel.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Similar to the relative spruce decline described in Fig. 5, the absolute decline (Fig. S6 in Supplementary material) showed a varying response for different classes of FN. In some cases, the overall trend varied greatly between classes. The most evident exception to this observation was the predictor SPRUCE10, which consistently affected the absolute decline in a similar way (i.e., increasing spruce decline along the gradient).

  Discussion 

Geology and forest naturalness as modulating factors

In this paper, we quantified the recent (2010-2018) Norway spruce (Picea abies) decline in Slovenia caused by abrupt disturbance events of widespread proportions, such as bark beetle outbreaks, ice storm damage and windthrows. We accounted for a wide array of environmental variables that could potentially influence spruce decline, including those related to geology, climate, forest soil, terrain topography and stand features. A key focus was given to the effects of two qualitative explanatory variables with a contrasting nature: bedrock type and forest naturalness (preservation of tree species composition). While the latter has been frequently studied and can be greatly modified by forest management ([9]), researchers have rarely addressed how given geological conditions (natural factors beyond the control of anthropogenic interventions) might influence forest decline.

One general observation derived from analyses was that relative spruce decline significantly differed with respect to bedrock. On average, the highest decline was observed for the carbonate bedrock type. This agreed with our hypothesis that forest compartments on carbonate substrate would exhibit higher spruce decline due to the less favourable soil conditions (lower water storage capacity). Compartments on silicate bedrock had the lowest decline. Although the western part of Slovenia receives higher annual precipitation, these areas are more susceptible to spruce decline due to high prevalence of shallow forest soils on limestone with low water availability and frequent rainfall deficits during summer. Higher temperatures weaken Norway spruce trees, particularly in areas outside the species’ natural range ([20]). Abiotic stress induced by drought may push trees beyond critical physiological thresholds for survival ([2]). Warming temperatures are also expected to facilitate the population development of lethal forest insects such as bark beetles ([37], [8]). In the last few decades, two main outbreaks of bark beetles have occurred in Slovenia. The epidemic period from 2003 to 2008 was initiated by a severe drought in 2003 and the period from 2014 onwards was triggered by an ice storm ([9]). Hlásny et al. ([18]) reported a transition from wind- to drought-driven bark beetle outbreak dynamics in the Czech Republic. Marini et al. ([21]) suggested that drought and climate warming can induce epidemic eruptions even in the absence of abrupt pulses of widespread tree mortality caused by discrete disturbance events.

After the ice storm in 2014, the Dinaric, Pre-Dinaric, Alpine and Pre-Alpine regions became epicentres of spruce decline ([19]). Some of these areas experienced an almost complete depletion of spruce, particularly at lower altitudes. Most of the plantations with a high proportion of Norway spruce trees are at a lower altitude, which means that trees are also often stressed by higher temperatures. On a positive note, shifts in forest composition and structure via host tree removal can counterbalance the enhancing effects of accelerated beetle population development and increased drought-induced susceptibility of spruce to beetles in the long run ([33]). The most prominent reduction of spruce growing stock was, at least to some degree, spatially concentrated (western part of the country, i.e., Dinaric Mountains, mid-altitude hilly landscapes between the Alps and Dinarides) rather than scattered across the country because of the joint effects of two negative factors that were spatially and temporarily synchronized: ice storm damage and drought-prone sites on carbonate bedrock. High-elevation areas on carbonate with high annual rainfall (between 2000 and 2500 mm per year or even more; e.g., calcareous SE Alps) experienced higher spruce decline compared to high-elevation areas on silicate bedrock in the eastern part with less annual rainfall (between 1500 and 1600 mm per year or even less; e.g., Pohorje Mountains). Owing to carbonate bedrock with low soil water storage capacity, the sites are more vulnerable to drought stress compared to sites on silicate bedrock, which can maintain sufficient water supply throughout the summer. Our analysis did not include the period 2019-2021, but from our personal observations and available data from SFS Annual Reports, we can claim that bark beetle infestations of spruce forests also started at higher altitudes because of a rise in temperatures, prolonged summer droughts and forest disturbances (e.g., windthrows). Increasing rates of tree mortality linked to a warmer and drier climate in high-altitude environments have been illustrated by several articles ([2]).

Historical management practices oriented towards economic benefits of forest stands have largely changed the composition and structure of today’s forests. A large proportion of Norway spruce was the main cause of the change in tree species composition of Slovene forests. There was a clear association between the degree of forest naturalness and the proportion of Norway spruce (see Fig. S3 in Supplementary material, bottom left panel). A higher proportion of spruce was characteristic of compartments with a greater alteration of tree species composition. A similar finding was confirmed by De Groot et al. ([9]). The share of spruce in forest compartments was also the decisive factor for the degree of damage associated with forest dieback processes in Sudety and Beskidy Mountains ([3]). We initially hypothesized that compartments with a greater alteration of tree species composition would show a higher spruce decline (Tab. 2). However, the results did not support this assumption as compartments with preserved composition exhibited the highest relative spruce decline among all classes.

A greater alteration of tree species composition is inversely related to forest naturalness ([5]), although naturalness does not automatically mean greater tree species richness as it depends on the potential natural vegetation. Forest ecology theory predicts that species-rich forests with more preserved species composition should show higher degree of temporal stability. Resistance and resilience to the negative effects of natural disturbances can be greatly elevated with mixed and local site-adapted tree species composition. In contrast, changed or completely altered tree species composition are often predisposing factors to negative disturbance impacts. One possible explanation for the discrepancy observed with our data analysis is that intense large-scale forest disturbances are able to override the buffering potential of favourable conditions seen in more preserved tree composition. Another explanation is that small, isolated spruce stands within more preserved, predominantly deciduous forest compartments are probably at even higher risk of significant relative change or even complete decline than larger stands of spruce in altered, spruce-dominated compartments. Because of the responsive and well-organised forest service, the likelihood of total decline of such spruce-dominated compartments, which also have higher economic value, is relatively low. The response of spruce decline to forest naturalness could thus be partially offset by forest management interventions. Taken together, the mitigating effect of bedrock type was larger compared to that of forest naturalness.

The hypothesized effects of individual quantitative predictors were also only partly supported by our results. Compartments with higher potential evapotranspiration had higher spruce decline on average, which is in line with the results from Hlásny et al. ([18]). This confirms the general expectation that drier and warmer climates are less favourable for Norway spruce. This tree species is more vulnerable to secondary infestations at sites with lower water availability. In comparison to the PET effect, other predictors were either less influential (e.g., heat load index, stoniness/rockiness), or the relationship between response and predictor variables requires careful ecological interpretation (e.g., higher spruce decline in compartments with higher soil clay content). In soils with a high clay content, the water available to plants can be low despite a better water holding capacity because the water is stuck in the fine pores. In addition, soils with a high clay content are prone to waterlogging, which leads to shallow root systems in Norway spruce, making the species susceptible to windthrow or drought effects.

Study limitations

Although we provide a countrywide analysis of the Norway spruce growing stock decline in relation to a broad spectrum of ecological factors, our study has some shortcomings worth to be discussed. Firstly, the low percentage of deviance explained in the GAM-based regression modelling (average was 7.1% for models based on bedrock type and 7.0% for models based on forest naturalness) suggests that some other ecological factors not accounted for in our study likely play a decisive role. The analyses were nevertheless done by maximizing the availability of compartment-level data. On a positive note, the constructed GAMs are fairly simple as they all include only one quantitative predictor, one qualitative predictor and the interaction between the two. When simultaneously including multiple quantitative predictors in the model, the explained deviance reached almost 20% (data not shown).

Secondly, the values of some predictors (e.g., PET and CLAY) were taken from an open-source database with a spatial resolution not fully congruent with the resolution of the forest compartments. In addition, caution is needed with the four classes of forest naturalness as they show quite uneven (unbalanced) ranges, and the reliability of their compartment-level determination might be questionable. Compartments can include various forest stands, but we used a single predictor and response variable value for each compartment, and this might have resulted in the loss of information by upscaling from a smaller (stand) to larger (compartment) spatial scale. A whole analysis based on stand-scale data would likely be more precise.

Thirdly, the investigated period 2010-2018 included a series of different disturbance agents that usually show a differential response to the ecological conditions of forest sites. A diverse combination of disturbances complicates the identification of the proximate causes of spruce decline. We would have gained deeper insights if we had included volumes of salvage logging attributed to different disturbance agents that contributed most to the observed spruce decline (bark beetles, ice storm, windthrows). A possible upgrade of our approach could include temporal dynamics (i.e., annual time series analysis) for individual disturbance agents.

Management implications

Despite the aforementioned limitations, our main findings have important forest management implications. Norway spruce is widely appreciated among foresters and forest owners because of its broad ecological niche, adaptability and valuable timber. However, past and present planting of spruce can create more problems than benefits, particularly in secondary forests on unfavourable sites for spruce growth. These problems seriously question the sustainability of silviculture oriented towards Norway spruce that is broadly rooted in management across Central and SE Europe ([17]). Based on our results, we can conclude that bedrock type is an important mitigating factor of spruce decline. However, it remains unclear whether such a buffering effect would be apparent if the sites on silicate bedrock experienced similar ice storm damage as was the case for compartments on carbonate or mixed geology. Massive insect outbreaks and other intensified disturbances might offset the potential buffering effect of bedrock geology. Large-scale disturbances are likely to overwhelm the meaningful mechanisms related to ecological predictions, as demonstrated by the study of De Groot et al. ([7]). When the level of severity of a disturbance crosses a certain tipping point, the potential of ecological factors to mitigate the damaging effects of natural disturbances and to positively contribute to forest conservation is reduced.

Under epidemic situations, forest sites theoretically less prone to disturbances (e.g., in the case of bark beetle outbreaks: lower temperatures and precipitation deficits, higher clay content with greater soil water storage capacity but not necessarily a greater amount of plant available water, north-facing slopes with a more humid microclimate, uneven-aged structure with mixed tree species) can exhibit substantial tree mortality and stand growing stock declines. Such offsetting effects likely also occurred in Slovenian forests where spruce decline was high on sites with a high degree of forest naturalness. This suggests that appropriate management oriented towards local site-adapted stand composition and structure might lose its ability to increase resistance against disturbances. However, we strongly advocate silvicultural measures that are regularly executed in timely manner, adapted to local site conditions as well as to the landscape-level configuration, and contribute to the greater bio-physical stability of forests. Regular cutting decreases the amount of salvage logging because of bark beetles. Nevertheless, the realization of possible harvest is generally below planned volumes ([9]). Preserved tree species mixture and active forest management reduce the amount of sanitary felling ([10]).

  Conclusions 

Recent evidence of Norway spruce profound disturbance-induce damage and consequent stock depletions appears to be a universal phenomenon in almost all regions in Central and South-eastern Europe. These effects have reinforced justified doubts regarding the sustainability and prospect of this tree species in the future. However, Slovenia is characterized by specificities not common to other Central European countries, such as the 2014 ice storm disturbance that escalated into massive bark beetle outbreaks in the following years. In this study we integrated multiple disturbance agents causing spruce mortality and a growing stock decline in relation to relevant ecological factors.

We found mixed support for ecology-based hypotheses. In agreement with our expectation regarding the influence of bedrock type, we found a lower decline in silicate compartments compared to carbonate sites. The results pointed to unfavourable site conditions for spruce (e.g., carbonate parent material with shallow soils) that likely increase its susceptibility to abiotic and biotic disturbances. In relation to forest naturalness, we did not find support for the hypothesized effects as preserved compartments showed the highest spruce decline. The mismatch between theoretical concepts and observed patterns might be explained by the fact that if intensive large-scale disturbances cross critical thresholds, they have the potential to cancel out the buffering capacity of certain ecological factors. This observation calls into question the ability of forest naturalness to mitigate disturbance effects. Nevertheless, active and local site-adapted forest management still undoubtedly represents the main tool in the fight against threats posed by climate change.

  Acknowledgements 

This research study was supported by Research Program P4-0107 and Grant no. Z4-4543, both funded by the Slovenian Research Agency (ARRS), and by the Public Forest Service (JGS Task 1), funded by the Ministry of Agriculture, Forestry and Food of the Republic of Slovenia. We would like to thank many colleagues from the Slovenia Forest Service (ZGS) for providing various data on the forest stands and forest compartments used in this study. We are grateful to Philip J. Nagel for proofreading and language editing. We would like to thank the reviewers for their thoughtful comments and efforts towards improving the manuscript.

  References

  Supplementary Material

Fig. S1 - Correlation matrix for five numerical predictors of spruce decline.
Fig. S2 - Statistical tests for differences in quantitative predictors between bedrock type categories.
Fig. S3 - Statistical tests for differences in quantitative predictors between classes of forest naturalness.
Fig. S4 - Relative number of forest compartments for bedrock type categories and classes of forest naturalness in subset 1 vs. subset 2.
Fig. S5 - Estimated response of the absolute spruce decline (m³) to selected ecological variables, where bedrock type was included as a categorical predictor in GAMs.
Fig. S6 - Estimated response of the absolute spruce decline (m³) to selected ecological variables, where the degree of forest naturalness was included as a categorical predictor in GAMs.