Modelling Response of Norway Spruce Forest Vegetation to Projected Climate and Environmental Changes in Central Balkans Using Different Sets of Species

Abstract

The structure and function of many forest ecosystems will be modified as a result of air pollution and climate change. Norway spruce (Picea abies L.) forests are among the first terrestrial ecosystems to respond to this change. We analysed how changes in climate and environmental factors will affect vegetation cover in Norway spruce forests and whether it is possible to assemble a list of diagnostically important/sensitive species that would be the first to react to changes in habitats of Norway spruce in Central Balkan. Significant changes in the vegetation cover of Norway spruce forests are mainly influenced by temperature increases (≈4 °C), and precipitation decreases (≈102 mm) by the end of the 21st century. Projections show that vegetation cover changes and future habitat conditions for Norway spruce forests on podzolic brown soils with a low base saturation and soil pH decreases, and temperature growth and precipitation decline, with the worst in the Rodope montane forest ecoregion. In Dinaric Mountain and Balkan mixed forest ecoregions, the range of natural occurrence of Norway spruce forest will shift to higher altitudes, or to the north. One of the cognitions of this paper is that, through available environmental models and their indices, species from the IUCN Red List should be recognised more properly and included in model calculations.

1. Introduction

Mountainous habitats in Europe contain ∼20% of the native flora [1]. They are authentic hotspots of plant diversity, hosting highly specialised vascular plants [2,3] and many endemics [4]. These ecosystems are sensitive to pollution and climate change, as major external environmental factors influencing the forest ecosystems in central Europe, due to sharp transitions in vegetation sequences (ecotones), with often endemic species and a limited species pool [5,6].

Forests are particularly sensitive to climate change because the long lifespan of trees does not allow for rapid adaptation to environmental changes [7]. Common broadleaved species appeared to be less sensitive to climate change, while some of the most important conifer species had already lost a large portion of their climatically suitable habitat [8]. In addition, the Norway spruce (Picea abies L.) migrated from several refugia following the last ice age [9].

The total natural distribution of P. abies covers 31 degrees of latitude from 41°27′ N (the Balkan Peninsula) to 72°15′ N (near the Chatanga River, Siberia). The longitudinal range is from 5°27′ E (the French Alps) to 154° E (the Sea of Okhotsk in Eastern Siberia). The vertical distribution is from 0 m to the altitudes above 2300 m (Italian Alps) [10]. As one of the most widespread types of forest in the higher mountains of Europe, Norway spruce forests prefer cool and moist climatic conditions. In the context of future climate change, this ecologically and economically valuable species may become severely affected under global warming conditions [11]. According to previous studies, the direct effect of air pollution was likely the main and primary stressor causing the weakening of the spruce trees and increased sensitivity to other stressors causing massive Norway spruce forest declines in the past few decades [12,13].

Several factors that are associated with climate change, which can act either independently or in combination, affect a forest ecosystem’s structure and dynamics [7]. They can alter species composition and richness through many pathways: temperature, precipitation, nutrient availability, carbon sequestration, soil chemistry, etc., and consequently lead to homogenisation. By changing gradually, air pollution has a cumulative effect on the biodiversity of natural ecosystems [14]. The atmospheric deposition of reactive nitrogen (N) is considered a global threat to biodiversity [15,16,17,18] and plant N availability over stand development [19].

In the Central Balkans, Norway spruce forests belong to the G3.1E EUNIS type (southern European P. abies forests). They occur in the southernmost part of the species’ range, with clumped distribution determined by the cooler climate of the highest mountains in Serbia (Figure 1). In this part of the European range, spruce mainly grows on shallow mountain soils [20] and on steep slopes, additionally stressed by summer droughts and forest fires, storms, and bark beetles [21].

The evaluation and quantification of the cumulative effects of air pollutants and climate change stressors are very complex and demand dynamic model assessment and different/specific biodiversity metrics systems to be applied [18,22]. According to the results of the MOTIVE project, increasing temperatures and a change in the precipitation regime may lead to a decrease in the area of mesic, cold- and humidity-adapted species such as Norway spruce, which are assumed to move northwards and lose large portions of their growing space. The simulation results for the Norway spruce (for the period of 2051–2080) specifically show that only habitats at higher elevations and latitudes are projected to remain suitable [23].

In this paper, the VSD+PROPS model and Habitat Suitability Index (HSI) were used to predict the influence of future climate changes (A1B scenario) and N deposition on the quality of habitat that might affect species composition. These indicators have been successfully used for policy support of the Convention on Biological Diversity (CBD) and Convention on Long-Range Transboundary Air Pollution (CLRTAP) in Europe and North America, as a link between air pollution and biodiversity. However, there is a lack of clear recommendations given by the CBD and CLRTAP legislation on a representative set of species as an indicator of climate and environmental change for certain habitat types.

One of the cognitions of this paper is that, through modelling processes, UNESCO Red List species are not properly recognised in the calculation of the banded HSI. A separate indicator for biodiversity change assessment should be assigned for Red List species with an additional weighting method.

The aim of this study was to disentangle the following aspects:

• How climate change, according to the projections, would affect vegetation structure in the Norway spruce forests of the Central Balkans;
• How changes in environmental factors and habitat suitability would affect vegetation structure in the named forest habitats;
• Defining a list of plant species as the most suitable for dynamic modelling of HSI metrics and early diagnosis of Norway spruce plant community reaction to multiple stressors.

Projecting the future vegetation structure of Norway spruce forests is of great importance for forest management. Adaptive management plans should acknowledge higher vertical and horizontal structural diversity since it improves resilience and resistance [24]. Another management strategy to be considered is the identification and protection of habitat refugia into which future populations of threatened species could migrate, or use as corridors for dispersal [25].

2. Materials and Methods

2.1. Study Area

This study was conducted in Norway spruce forests at four mountains that belong to three ecoregions in the Central Balkans: Rodope montane forests (Stara Mountain), Dinaric montane forests (Mountain Zlatar), and Balkan mixed forests (Mountain Golija and Kopaonik) (Figure 1). According to The Digital Map of European Ecological Regions- DMEER [26], ecoregions are ecologically distinct areas, with relatively homogeneous ecological conditions on the basis of updated knowledge of climatic, topographic, and geobotanical data within which comparisons and assessments of different expressions of biodiversity are meaningful [27]. Basic information on the studied sites is presented in

Table 1. Basic information on studied sites.

Site	Stara Mt.	Zlatar	Golija	Kopaonik
Ecoregions	Rodope montane forests	Dinaric montane forests	Balkan mixed forests	Balkan mixed forests
Coordinates	43°17′34.33′ N 22°47′09.76′ E	43°24′20.04′ N 19°48′36.19′ E	43°20′6.38′ N 20°17′12.21′ E	43°17′30′ N 20°48′50′ E
Altitude (m)	1450–1700	1220–1400	1300–1720	1520–1800
Precipitation (mm) *	606	826	797	770
Temperature (°C) *	6.84	6.00	5.38	5.86
Soil parent material	Permian red sandstone	Chromic Humic Cambisol/Chert	Phyllite	Granodiorite/Granite
Soil type (WRB)	Podzolic Cambisol	Chromic Humic Cambisol/Dystric Cambisol	Podzolic Cambisol/Dystric Cambisol	Podzolic Cambisol Humic Cambisols
Protection status	Nature Park	Landscape of exceptional natural beauty	Nature Park	National Park
EUNIS habitat type	G3.1E—Southern European P. abies forests
EU HD type	9410 Acidophilous Picea forests of the montane to alpine levels (Vaccinio-Piceetea)
EU HD subtype	42.24 Southern European Norway spruce forests

* Annual average 1961–1990; WRB—World Reference Base for Soil Resources; EUNIS—the European Nature Information System; EU HD—European Union Habitat Directive.

.

2.2. Vegetation Sampling, Data Processing, and Model Description

Vegetation in Norway spruce forests was analysed using the Braun–Blanquet phytocoenological method [28]. A total of 67 relevés were taken from all four sites (400 m² plot size), and 309 species were recorded (Stara Mt.—35 species/10 relevés; Zlatar- 120/22; Golija- 76/27 and Kopaonik- 78/8). The species nomenclature follows the Flora Europaea (Euro+Med) plant base and the plant list.

In order to identify the most useful set of species as an indicator of climate and environmental changes, a list of three groups of species was assembled. Up to 40 of the most frequent species from the relevés per site were extracted as dominant species (dom). Typical species (typ) for the Norway spruce habitat were extracted according to the available literature [29,30,31]. Finally, the EuroVeg Checklist expert system Juice 7.0 [32] was used for extracting the list of diagnostic species (dia) for Vaccinio–Piceetea forests (Table S1, Supplementary Materials). The term diagnostic species refers to species which preferably occur in a single vegetation unit (character species) or in a few vegetation units (differential species), and are particularly useful for the identification of vegetation units [33]

VSD+Studio model, Version 5.5.1, Alterra CEE [34], coupled with the PROPS model (Probability of Occurrence of Plant Species Version 5.5.1. Alterra, CCE) [35], was used for the simulation of ground vegetation cover change due to the changes in climate and soil chemistry dynamics [36]. In conjunction with the VSD+Studio model, two pre-processor models were used. First, MetHyd v1.5.1, CCE [37] was used for calculating evapotranspiration, soil moisture and percolation, and parameters related to nitrification, denitrification, and mineralisation, and secondly, the GrowUp model v1.3.2 Alterra, CCE [37] was used for computing N and base cation uptake from the specified tree-growth and specific forest management. The net nutrient uptake (equal to the removal of harvested biomass) was calculated using the average growth rates measured in the Norway spruce stands with natural rejuvenation forest management in protected areas.

As an output of VSD+PROPS coupled models, the occurrence probability (OP) of the species has been defined as a function of soil chemistry and climatic variables [34,38]. Sørensen’s similarity index (β) [39] was used for the analysis of floristic differences between the OP of the set of species (dom, typ, and dia) from the time of phytocoenological surveying, and modelled OP (for the same sets) in 10-year time slices (1980–2100). Sørensen’s similarity index (β) was used to address the following question: How will vegetation structure change over time, or specifically, what species group will be the first to tell if the overall structure is about to change?

As another output of the same coupled models, the Habitat Suitability Index (HSI) has been defined as the arithmetic mean of the ‘normalised’ probabilities of occurrence of the species of ‘interest’ [40]. In order to discover the connection between different sets of species and environmental changes in Norway spruce communities, three HSIs were calculated: HSI_dom, HSI_typ, and HSI_dia, which were the outputs when we used dom, typ, and dia set of species, respectively. A correlation was performed for the time-series data obtained from the model for the 2020–2100 period in order to determine the significance level of the main abiotic drivers (air temperature, precipitation amount, soil N content, soil pH) and their impact on habitat suitability change. All three HSI indices were included in this analysis (HSI_dom, HSI_typ, and HSI_dia), with OP of Norway spruce, annual temperature, precipitation, N deposition reduction, and soil pH.

2.3. Climatic Data and Modelling

Coupled regional climate model EBU-POM [41,42] was used for the dynamical downscaling of global climate change projections, following moderate IPCC emission scenario A1B [43], in order to obtain a time-series dataset on future climate conditions. The whole model simulation covers the period from 1951 to 2100. The climate variables used in this modelling process were mean monthly temperature, precipitation, insolation, and soil water content. For the surface air temperature and precipitation, model bias was removed using an observed time series of the same variables (on a daily level), taken from the Republic Hydrometeorological Service of Serbia (RHMSS) database, over a normal climate period of 1961–1990, following a statistical bias correction method [44]. For the insolation and soil water content, the same approach was applied using the data from ERA-Interim reanalysis [45] from the European Centre for Medium-Range Weather Forecasts.

Climatic changes were analysed by comparing variables of the 1961–1990 normal period (Republic Hydrometeorological Service of Serbia), with their projected values in the 2071–2100 period (EBU-POM). The time-series data of variables for the 1980–2100 period were used as an addition to the MetHyd model, for the VSD+Studio model run.

2.4. Soil, Deposition Data, and Modelling

Soil data were obtained from the archive of the Faculty of Forestry, University of Belgrade (Serbia). The main physical and chemical soil properties are presented in Table S2, Supplementary Materials. Nitrogen (N) and sulphur (S) deposition data (1980–2013) were obtained from the database of The European Monitoring and Evaluation Programme (EMEP) model results [46], using data from the grid cell corresponding to the location of the relevés (Figure S1, Supplementary Materials). Coupled biogeochemical model VSD + Studio was used for the simulation of future climate, as well as N and S deposition effects on soil chemistry dynamics [34]. It was presumed that the deposition of N and S would stay on the current trend (the average value for the last 10 years); therefore, a corresponding deposition scenario was used, which is in accordance with the A1B climate scenario. Thus, N and S deposition will show a similar decreasing trend, in all studied areas in the 1980-2100 time period.

Canonical Correspondence Analysis (CCA) ordination technique was performed with Past 4.06b software [47], to jointly display the obtained results from the EBU-POM and VSD+ Studio models, for the period 2020–2100, with decade values as sites and OPs as species scores. The analysed environmental variables were temperature (Temp C), precipitation (Precip), soil moisture content (Theta), soil pH value, C:N ratio (CN rat), and plant-available nitrogen.

3. Results

3.1. Pre-Modelling Analyses

All analysed forest communities belong to Piceion excelsae, Pawłowski et al., 1928 alliance, according to the EuroVeg Checklist [48]. In, or in the range of, these forests, there are 15 endemic and Tertiary relict species, and 1 endemic genus (e.g., Bruckenthalia spiculifolia, Pancicia serbica, Aremonia agrimonioides, Asarum europaeum, etc.). The temperature and precipitation for a normal 30-year period (1961–1990) and EBU-POM model projections for the A1B scenario, 2001–2100 and 2071–2100, are presented in Tables S3 and S4 (Supplementary Materials). Very dry and warm summers (covering the whole vegetation period), according to the model simulations, in all studied areas are projected. The temperature will shift by nearly 4 °C until the year 2100 (Table S3, Supplementary Materials), and annual precipitation will decrease in the range of 85.78 mm (Rodope montane forests) to 113.39 mm (Balkan mixed forests—Golija) (Table S4, Supplementary Materials).

3.2. Results of the Norway Spruce Forest Modelling

The analysis of Sørensen’s similarity index (β) showed that, for the dom species, β will be decreasing constantly until 2090, when it should slightly begin to rise (Figure 2a). Considering all studied sites, the β median is 0.54 (min. 0.48 in 2090, max. 0.54 in 2100), meaning that on average 46% of dominant species will change over time. The structure of dominant vegetation will change the slightest in Rodope montane forests (β = 0.93 in 2030), and the strongest in Balkan mixed forests (β = 0.24 in 2040, Golija site).

Changes in typ species are less pronounced. The median β for all sites is 0.62 (min. 0.58 in 2040, max. 0.82 in 1990), and it is pretty much constant throughout the simulation period (specifically from 2010 to 2100). The only ecoregion that stands out is Rodope montane forests, with very high β = 0.95 until 2030, when it will start to decrease to 0.82 in 2100 (Figure 2b).

The changes in vegetation structure are probably the most noticeable when we consider dia species, although the overall median β is the same as with dom species—0.54 (min. 0.53 in 2060, 2070, and max. 0.75 in 1990). Considering all analysed sites, almost half of the diagnostic species were already ‘changed’ in the 1980–2010 period (Figure 2c), and according to the simulation, β will remain almost constant until the end of the century, when it will rise to 0.65. Only Rodope montane forests will not follow such a trend; it appears that β will start to decrease near the end of the 21st century (min. 0.86 in 2090, 2100).

3.3. Habitat Suitability Analysis

HSI was used as a measure of the influence of future climate changes (A1B scenario) and N input on the quality of habitat that might affect the dom, typ, and dia species occurrence. The HSI_dom value increases for the Dinaric and Rodope montane forests and partly in Balkan mixed forests after 2008. Until the end of the simulation period (2100), it remains balanced with only small oscillations (Figure 3a). The HSI_typ increases in Balkan mixed forests after 2008, up until 2083. After 2083, it decreases again to the end of the simulation period (2100). For the Rodope mountain forests, the model shows the earliest (from 2029) and largest decrease (from 0.54 to 0.34) of HSI_typ. The value of the HSI_dia is lower in comparison to the HSI_dom and HSI_typ for all ecoregions, with the exception of the Kopaonik site (Figure 3a–c). It should also be noted that, after the 2080s, the Balkan mixed forest ecoregion (for all three groups of species: dom, typ, and dia) will experience a decline in habitat suitability.

After the correlation analysis, we compared obtained indices across each other (

Table 2. Correlation matrix for Habitat Suitability Indices, P. abies occurrence probability and environmental variables for period 2020–2100.

Sites	Variables	HSIdom	P. abies	T	P	N	pH
Stara Mt.	HSIdom	/	0.648	−0.579	0.492	−0.638	0.748
	HSItyp	0.584	0.944	−0.736	0.617	n.s.	n.s.
	HSIdia	0.474	0.745	−0.571	0.673	n.s.	n.s.
	P. abies	0.648	/	−0.839	0.572	n.s.	0.46
Zlatar	HSIdom	/	0.802	−0.554	0.778	−0.646	0.708
	HSItyp	0.817	0.77	−0.441	0.691	−0.52	0.498
	HSIdia	0.539	0.78	−0.501	0.420	n.s.	n.s.
	P. abies	0.802	/	−0.85	0.592	−0.467	0.608
Golija	HSIdom	/	0.679	n.s.	0.674	−0.801	0.855
	HSItyp	0.966	0.7968	n.s.	0.611	−0.718	0.742
	HSIdia	0.939	0.711	−0.31	0.629	−0.741	0.81
	P. abies	0.679	/	n.s.	0.533	−0.381	0.375
Kopaonik	HSIdom	/	0.588	n.s.	0.66	−0.825	0.862
	HSItyp	0.735	0.799	−0.203	0.606	−0.567	0.452
	HSIdia	0.897	0.571	−0.432	0.55	−0.776	0.886
	P. abies	0.588	/	−0.612	0.703	−0.499	0.484

T—annual temperature; P—annual precipitation; N—soil nitrogen content; pH—soil pH value. All correlations are significant at level p < 0.05; n.s.—not significant.

).

The highest positive correlation was found between HSI_dom and HSI_typ in the Balkan mixed forest (Golija site) ecoregion, and the highest negative correlation (−0.85) was found between the occurrence probability of Norway spruce (Zlatar site) and mean monthly temperature.

There was a strong correlation found between the HSI_dom and habitat suitability of two other groups of species (typ and dia), indicating that all three groups prefer similar habitat conditions, and their species overlap. Additionally, OP of Norway spruce was strongly correlated with all three HSIs, indicating their mutual dependence. Regarding the environmental variables, HSI_dom had a strong positive correlation with annual precipitation and soil pH and a strong negative correlation with annual temperature and N deposition reduction (Table 2). The other two HSIs were similar. The only environmental parameter strongly correlated with all three HSIs is the annual precipitation amount. This result confirms the strong dependence of Norway spruce forests on both aboveground and underground water availability.

Observing the relationship between OP of Norway spruce and environmental variables, the only strong correlation was found for precipitation (Table 2). That led us to the conclusion that Norway spruce, in particular, is the species that reacts the most to a lack of precipitation; thus, the entire habitat becomes less suitable for the entire Norway spruce forest.

3.4. Canonical Correspondence Analysis

The results of CCA gave us a perspective of the species–environment changes over time (Figure 4a–d).

Temperature and C:N ratio are expected to increase, while precipitation, soil moisture content, plant-available N, and soil pH are expected to decrease at all studied sites. An exception is the Dinaric montane forests, where pH values are deviating but remain around 4.7. The C:N ratio will drastically increase by a median of 9.97 (median for 1980 = 11.41, median for 2100 = 22.67) in all analysed ecoregions. Values of plant-available N, as previously stated, will follow the adopted scenario, the CCA confirms. The decrease in soil water content will be very moderate for the analysed ecoregions. It will be less than 2.18% in Dinaric montane forests, to 4% in Rodope montane forests, according to the simulation. Soil pH, although very low (median for 1980 = 4.06, min 3.90, max 4.46), will decrease even more (median for 2100 = 3.92, min 3.71, max 4.11) at all three studied ecoregions. Only a few acidophilus species will prefer these kinds of soil conditions, e.g., Vaccinium myrtillus (V_myr), V. vitis-idaea (V_v-i), Juniperus communis subsp. Nana (J_c_al), while species such as Luzula sylvatica (L_syl), Hieracium murorum (H_mur), and Valeriana tripteris (V_tri) are not tolerable towards low pH and are expected to be near extinction. In Dinaric Mountain forests (Figure 4b), species are concentrated along the axis that represents the pH gradient, indicating that the species distribution is strongly related to soil pH changes. The soil pH appears to have a strong role in the habitat of the Balkans mixed forests as well, where the majority of species are also distributed along this gradient axis (Figure 4d).

4. Discussion

4.1. Norway Spruce Forests’ Response to Climate Changes

According to the simulations undertaken by Mihailović et al. [49] for the A1B scenario, a warmer and drier climate in the region in the study area can be expected in the future. As the overall climate warms, it can be expected that the current Dfwbx’’Köppen zone will shift to higher altitudes in the mountain ranges. As certain types of vegetation are associated with specific climate types (which are shifting, thus warming the surface at higher altitudes), certain native species may no longer grow well in a given altitudes’ climate [49]. In the future, under the A1B scenario, the climate zones of the Central Balkans will shift in altitude relative to the period of 1961–1990 for 100 m (for period 2001–2030) and 400 m (for period 2071–2100). The results of our study completely correspond to those described in the study above, with the Norway spruce forest range expected to shift as well.

High-mountain ecosystems are increasingly threatened by climate change, resulting in loss of biodiversity, habitat degradation, and landscape modifications [50,51]. Although Norway spruce forests are relatively stable ecosystems (in the coenological optimum) and not unduly susceptible to degradation, any long-term negative impact may lead to a reduction in the canopy density, poor rejuvenation, and, ultimately, acidification and further soil physical degradation [31]. As a result, Vaccinium myrtillus and Juniperus sybirica spp. nana will potentially spread throughout these degraded habitats, until the last degradation stadium, Nardetum strictae s.l., is reached.

Boreal forests have a large impact on local temperatures and, of the three main forest types, the largest influence on global mean temperature. Our study shows that the temperature in analysed locations will increase by 3.87–4.00 °C, which is in conjunction with the fourth and fifth IPCC Assessment Reports [52,53]. The precipitation will also decrease, in the range of 12.31% to 14.63%, by 2100. It is expected that a climate warming of 1–2 °C can cause small changes in alpine vegetation, but even greater warming of the climate can cause more pronounced changes [54]. This increment will be met by increased mineralisation, nutrient availability, and uptake by roots, without which the growth response to elevated temperature and CO₂ will stagnate at a lower level [55,56,57,58,59].

Precipitation usually has a strong impact on the dynamics of soil moisture. Its loss by an interception in spruce stands is related to the nature of the rain and forest age. Many studies confirmed that the interception is increasing with the stand age [60,61]. Canopy gap differences were recorded in the dripping zone at the crown periphery and in the central zone of a crown of the stand [61]. Some studies of spruce forests show that the soil water content in the surface layer is higher in older stands, whereas in deeper soil it is higher in the thickets and young spruce stands [60]. Norway spruce needs high relative humidity and high surface moisture [62]. Additionally, even in its seedling stage, Norway spruce is very sensitive to soil moisture change, according to our results (Figure 4). The results from our study are in compliance with the results by Mihailović et al. [63]. The authors argue that for Dystric Cambysol soil type, the decrease in soil moisture will be between 6% and 8% (relative to the period 2021–2050, or 2051–2100).

4.2. Norway Spruce Forests’ Response to Environmental Changes

The results show that sites in the same ecoregion have the highest correlation (Balkan mixed forests), and sites with the lowest similarity belong to those ecoregions (Rodope and Dynaric mountains) with the largest differences in temperature and precipitation values (Tables S2 and S3 in Supplementary Materials). The model output results also show similar trends for the future. The Rodope mountain ecoregion is highlighted as the least suitable for spruce forests (according to the A1B scenario) and the Dynaric mountain ecoregion, as the most suitable. This was additionally proven through the correlation of the HSIs and the main climatic parameters (Table 2). Thus, the delineation of the different habitats’ conditions, according to DMEER, for the comparison and assessment of different expressions of biodiversity are shown to be meaningful for the present and future analyses of plant communities.

A narrow range of tolerance is described as an indicator of a species’ sensitivity to environmental changes. If it is also an endemic or relict species, its response to climate change should be analysed very carefully. Species distribution shifts will also depend on habitat availability. Many already-threatened local endemic species will be affected first [64]. Some of them are range-restricted species and are in the list of characteristic combinations of species or play a pronounced role at the ground level. Some of them have a few generations in the vegetation period and cover the ground level with dense populations. They would be replaced by other species with less demand on the environmental conditions, which would, in turn, signal the first step in the process of the degradation of Norway spruce forests.

4.3. The Effects of Changes in Soil Properties

Norway spruce is typically found growing in acidic soil and is believed to be adapted to low soil pH [65], and it produces slowly decomposing, acidic coniferous litter [66]. According to the projection, by 2100, all studied sites will experience a further soil pH decrease. As previously stated, only a few acidophilus species can overcome such a low soil pH, and they can be considered diagnostic for spruce forests [33]. For all but the Dinaric ecoregion, pH will descend below 4, making those habitats more susceptible to further soil physical degradation.

The C:N ratio of forest soils is one of the best predictors for evaluating the soil functions mainly involved in climate change issues [67]. The C:N ratio is more or less within the range for spruce forests in the studied ecoregions and soil types [68]. In spruce forest soil, the C:N is affected mostly by humus type, soil type, ecoregion, and parent material. In our study, according to the A1B scenario and corresponding temperature increase, the C:N ratio will increase by a median of 9.097, and it is expected that it will affect the occurrence of some plant species with a narrow range of tolerance to a C:N ratio change. Furthermore, the species will react to N deposition as well.

4.4. The Effects of Deposition Changes

Thus far, the deposition of N and S has decreased considerably as a result of political initiatives in Europe—Convention on Long-range Transboundary Air Pollution (CLRTAP)—with S reduction from 1980 ranging up to as high as 90%. There has been a small decline in N deposition in some regions in recent years but with occasional rises after 2006. The deposition of N remains high in many areas and critical loads are still exceeded in many parts of Europe. Many studies have observed high N inputs as a factor in soil biogeochemistry changes, community composition change, and, finally, loss of species diversity [69,70,71,72,73,74,75,76]. For all sites we have examined in this study, the total N deposition is within or near the lower limit of the recommended range (10–15 kg ha⁻¹ yr⁻¹) for spruce forests [77]. At the beginning of the research period, in 1980, N deposition was near the lower limit of the recommended range, and for the period 2020–2100, it is constant and below the lower empirical range for a spruce habitat. N deposition below 10 kgha⁻¹yr⁻¹ could still have a negative impact on most plant species and their occurrence within the spruce habitats (Figure S1, Supplementary Materials).

5. Conclusions

In conclusion, we can emphasise that according to the projections, the vegetation structure of Norway spruce forests will change. The outlook for Norway spruce forests is the worst on podzolic brown soils, with a low content of base saturation, temperature growth, precipitation decline, and soil pH decrease in the Rodope mountain forest ecoregion. Changes in environmental factors will result in large changes in vegetation cover and may further lead to physical soil degradation. No tree species would be able to adapt to these conditions, except maybe the woody species from the Juniperus and Vaccinium genus. In Dinaric Mountain and Balkan mixed forest ecoregions, the Norway spruce forest limit will change, and their range of natural occurrence will shift to higher altitudes, or to the north. The composition of species and the relationship between dominant and diagnostic species and the environment can be used as predictors for the determination of habitats suitable for Norway spruce.

Species of importance for biodiversity maintenance and conservation (e.g., endemic and relict plant species or species from the IUCN Red List), should be considered in habitat suitability modelling. Uncertainties in the habitat suitability modelling for these species are reflected in the following: not all species are represented in the used plant database, and these species are not usually particularly abundant, or they are not typical for the specific habitat. Consequently, they are not properly recognised in the calculation of the banded HSI. A separate indicator for biodiversity change assessment should be assigned for Red List species with an additional weighting method. Further research and monitoring of Norway Spruce habitats are needed in order to understand these relationships and to prove or disprove the simulation results.

However, there are more synecological and autecological processes in the very communities that might be initiated by changes in ecosystem composition, phenology, reproduction rate, tree growth, habitat disturbance, range, genetic diversity, adaptations, competition, and many others. Clumped distribution of Norway spruce forests in the Central Balkans signals the additional sensitivity of these habitats and those spatial changes in these fragmented parts would be useful to monitor in the context of climate change.

Forest ecosystem modelling studies should be based on reliable and long-term time-series data. All results should be validated with on-site observations in order to improve the level of model certainty, as support for adaptive forest management/conservation and decision making under climate change.