Spatial Shifts in Species Richness in Response to Climate and Environmental Change: An Adaption of the EUROMOVE Model in the Czech Republic

: Climate change has greatly altered plant habitats, resulting in greater biodiversity loss at different scales. Therefore, it is important to quantify such changes for better monitoring and conservation. In this study, we adapt the EUROMOVE model and its mean stable area indicator (MSAi) to the conditions in the Czech Republic. Our objective was to predict change in species richness from a representative pool of 687 species from 1990 to 2100 under the RCP 8.5 climate scenario, focusing on the current period (2018). Another objective was to assess the effectiveness of the MSAi as a tool for quantifying landscape vulnerability. Our result shows that species habitat expanded between 1990 and 2018, although about 2 per cent of species were lost. The average MSAi of the most favourable highland habitats may decrease from 0.85 to 0.65 by 2100 as >20% of baseline species may be lost. Indicator species of Alnus (alder) and Festuca (fescue), typical of lowland habitats, are among the most vulnerable, already showing a net loss of their current habitat extent. The MSAi can be applied as a comprehensive tool to quantify the impact of climate change on landscape vulnerability as more survey data becomes available.


Introduction
Climate change continues to accelerate species diversity loss, shifting species and their habitats northward [1][2][3][4][5][6][7]. Such changes have been captured at different spatial and temporal scales with the help of species distribution models SDM [8][9][10][11][12]. However, this evidence of change does not adequately quantify the state of landscapes for monitoring biodiversity loss and related ecosystem function [13].
The statistical correlation between species occurrence and the current environmental conditions is often the basis for projecting future distribution patterns from SDMs [2,4,7]. Subsequently, the results from SDMs can be used to choose conservation sites [7,14] and assess the risk of extinction [2,4,5]. In the case of multiple plant species, individual binary species maps can be additively aggregated into stacked species distribution models (S-SDMs) [15][16][17]. The approach is based on the fact that each species in a community will adapt differently to changing environmental conditions [17,18]. The importance of how modelled changes may affect ecosystems can be drawn from experimental studies showing a positive correlation between species richness, stability and functioning [19][20][21]. More importantly, there is evidence that above a 20% threshold, the impact of species loss on some ecosystem functions is comparable in magnitude to that of climate change and other major drivers of change [22]. Thus, detailed changes in species richness and habitat extent from S-SDMs is needed to comprehensively quantify and facilitate the application of S-SDMs in monitoring and conserving biodiversity [13,23].
Several biodiversity models and indicators have been proposed to quantify biodiversity loss at the global and regional scales. These metrics, including the mean species abun-dance (MSA) from the GLOBIO model [24,25], the biodiversity intactness index (BII) [26] and the living planet index (LPI) [27], have been tested at different scales. They are all variants of the natural capital index (NCI) indicator and measure population trends among species, ecosystems' state, and their degradation extent relative to a predefined historical species abundance [28][29][30]. A less commonly tested indicator, the mean stable area index (MSAi) from the EUROMOVE model, has been used to quantify the impact of climate change on plant species diversity within the European region [2,5]. EUROMOVE is an S-SDM based on multiple logistic regression between representative plant species and climate data for different climate scenarios in the European region. The MSAi indicator was used to quantify the stability of the entire European landscape on a scale of 0.0 for highly impacted landscapes to 1.0 for unimpacted landscapes [2,5]. Like most global and regional biodiversity indicators, MSAi was calculated from coarse climate data, ignoring local topography, hydrology, and geological variations. Hence, it is likely that the estimated MSAi values for most European countries may not reflect the actual situation. A detailed estimation at a national scale is necessary to improve conservation management.
The current locations of vascular plant species diversity in the Czech Republic, like in most parts of the European region, are restricted to national parks and protected areas in mountains and hilly areas [31]. In the Czech Republic, the distribution of high plant species has been classified into nine forest vegetation zones (FVZ) to reflect changes in microclimatic conditions and the transition from one type of species to another with changes in altitude [1,[32][33][34]. Oaks-beech species dominate in lowlands; fir-beech are common at moderate altitudes. Spruce-beech dominates in high altitudes, while dwarf pines are restricted to the most elevated areas [1,33,34]. Under a partially mitigated climate scenario, the FVZs have been projected to contract or shift to higher altitudes. Although their response under the baseline climate scenario is not well known, related studies have shown that high altitudes are equally at risk of losing their climate space by 2100 [35,36].
This study uses a representative sample of plant species to predict and quantify species richness patterns in the Czech Republic from 1990 to 2100 under the baseline climate scenario, focusing on the current period (2018). The main goal is to adapt the EUROMOVE model to the environmental conditions in the Czech Republic. We seek to assess the vulnerability of landscapes in the Czech Republic to species loss. Achieving this aim should improve our understanding of landscape vulnerability to biodiversity loss and its implications for policies and conservation management in the Czech Republic.
We present the results of a detailed multispecies distribution model capturing simultaneous changes in species richness and potential habitat extent with time. We also include the distribution of species indicative of distinct vegetation zones to understand the impact of climate and environmental factors. Our results confirmed that the most vulnerable areas in the Czech Republic are lowlands. Highland habitats are the most stable but are highly at risk of losing their climate space. Their current MSAi values range from 0.8 to 1.0 and may decrease from 0.5 to 0.8 by 2100 without efforts to mitigate climate change. About 2% of the baseline (1990) species pool was lost by 2018. However, the average MSAi for the entire landscape has remained the same. While MSAi is a broad and comprehensive tool for quantifying the impact of climate change on landscape vulnerability, there is a need to quantify vulnerability at the ecosystem level in detail. This paper is structured as follows: Section 2 describes data sources and the modelling approach used to achieve the research aim; Section 3 describes the results in four subheadings; Section 4 discusses the results, limitations and suggestions for future research and ends with concluding remarks in Section 5.

Materials and Methods
We processed and selected relevant species and environmental data in R (R Development Core Team 2021), QGIS 3.x and ArcMap 10.x software. Next, we modelled species distribution in R using Maxent [37,38]. Lastly, we used the model output to calculate the Diversity 2022, 14, 235 3 of 16 mean stable area index for the Czech Republic. A description of the data and modelling approach is presented as follows.

Species Data
We obtained georeferenced survey records (presence only) of higher vascular plant species from 1960 to 1990 (baseline species data) from the Agency for Nature Conservation and Landscape Protection of the Czech Republic. The data which excluded alien species consisted of 3044 species and~1.9 million records. We created 315,960 spatial sampling grids with a dimension of 500 m × 500 m for the entire Czech Republic to process species data. We joined each species' data to the spatial grids by location, counting the total number of each species in each grid (species abundance). We then reduced species abundance to presence-only records by assigning a value of 1 to each grid having species abundance ≥1 to estimate species richness. We removed all species with less than 50 observations for the entire country because we thought these could not be accurately modelled. Hence, this resulted in a total of 1117 baseline species.

Climate and Environmental Data
We obtained the following category of environmental data: climate, soil, topography, geology, hydrology and drainage and prepared them at a resolution of 500 × 500 m. The sources and scale of each environmental data are summarised in Table 1. The climate variables selected to reflect varying conditions of moisture and energy included: the mean annual temperature (antemp), mean annual sum of rainfall (anrain), mean annual temperature of the coldest month (tcold), mean temperature of the growing season above 5 • C (tempgs), the effective sum of temperature above 5 • C (efstemp) and the length of the vegetation season (lenvegt). We obtained the monthly averages of these variables from 1961 to 2099 as shapefiles or as raster layers. From this time span, we considered four modelling periods, namely, 1961-1990 (baseline climate), 1990-2018 (current climate), 2041-2060 (nearfuture climate) and 2081-2100 (distant future climate). Future climate data sets were based on the unmitigated representative concentration pathway -RCP 8.5 [39], which was the only projected data set of the Czech Republic with the finest resolution (500 × 500 m) available at the time of this study. The projected datasets were dynamically downscaled from the HadGEM 2-ES global model, chosen because it was the most accurate model, capturing change in precipitation patterns in the Czech Republic. These datasets excluded the unavailable monthly averages from 1961 to 1979. In ArcMap 10.x Geostatistical analyst, we interpolated the historical climate data using ordinary kriging and fitted a spherical model to each climate variable. The Pearson correlation coefficient between the observed and predicted historical climate variables ranged from 0.80 to 0.88.
We obtained non-climatic data in vector or raster format, which included: vegetation cover, soil texture, soil depth, geology, distance to water bodies, soil drainage, slope and aspect. Soil depth and texture were scaled from 1 to 2, where 1 represented shallow and coarse soils, respectively, while 2 represented deep and fine soils. Distance to water body varied from 15 to 3539 m. Soil drainage varied from 1 for poorly drained to 5 for welldrained soil. Geology considered different rock fragments and varied from 1 for coarser and heterogeneous rocks to 1.9 for finer and homogeneous rocks. Slope angle range was from 0 degrees for flat areas to~35 degrees for mountainous areas. All vector layers were converted to raster and resampled to match the resolution of the climate data. Before modelling, we removed redundant variables (annual temperature and the effective sum of temperature) having a variance inflation factor (VIF) > 10. We also eliminated vegetation cover data because it was missing in most of the chosen species locations.

Spatial Models
We adapted the EUROMOVE modelling approach to predict species richness from 1990 to 2100 and then calculated the mean stable area index (MSAi). The EUROMOVE model, as mentioned earlier, is a European regional model for higher vascular plants. The model was used to capture changes in species composition from the additive overlay of binary species maps for different modelling periods. It calculated the mean area covered by species for each modelling period and then aggregated these outputs into the mean stable area indicator. Thus, on the regional scale, the model predicted both appearing and disappearing species in each European country based on the regional species pool. In the case of the Czech Republic, species distribution was estimated on grids of about 500 m 2 . We replaced logistic regression with the maximum entropy algorithm-Maxent [37,38], accepting all default settings. Maxent optimises prediction by comparing the probability density of environmental conditions where a species is observed to the probability density of background environmental conditions in an area based on minimum distance [40]. Thus, we sampled 10,000 background points and their environmental attributes for each species and implemented Maxent in R from the dismo package [41]. Before modelling, we randomly partitioned (without replacement) 70% of each of the 1171 species datasets for model training and 30% for model evaluation. We used the trained model to predict species richness for subsequent periods. We also assessed selected species' response to climate-environmental gradients using Random Forest [42,43] and logistic regression [44,45]. Random Forest is based on an ensemble of decision trees built on different bootstrapped samples and predictors. The algorithm is robust to nonlinear variation and depends on the most votes from decision trees to classify unique cases in an observation. In contrast, logistic regression is a flexible extension of linear models with the possibility to handle nonlinear variation and different families of probability distribution in the data [44]. Thus, we studied response curves from the best of the three modelling approaches and further assessed landscape vulnerability as the net loss or gain of the simulated current and future climate space. The workflow is summarised in Figure 1. [44,45]. Random Forest is based on an ensemble of decision trees built on different bootstrapped samples and predictors. The algorithm is robust to nonlinear variation and depends on the most votes from decision trees to classify unique cases in an observation. In contrast, logistic regression is a flexible extension of linear models with the possibility to handle nonlinear variation and different families of probability distribution in the data [44]. Thus, we studied response curves from the best of the three modelling approaches and further assessed landscape vulnerability as the net loss or gain of the simulated current and future climate space. The workflow is summarised in Figure 1.

Model Evaluation
Model evaluation was performed to assess performance and classification ability. We calculated the area under the receiver operator curve (AUC) [46,47] and the true skill statistics (TSS) [48]. The AUC gives the overall performance of a mode. It is a probability measure of a model's ability to distinguish between presence and absence locations. The

Model Evaluation
Model evaluation was performed to assess performance and classification ability. We calculated the area under the receiver operator curve (AUC) [46,47] and the true skill statistics (TSS) [48]. The AUC gives the overall performance of a mode. It is a probability measure of a model's ability to distinguish between presence and absence locations. The AUC ranges from 0.5 for poor models without discrimination ability to 1 for perfect models [46,47]. The TSS is a robust threshold-dependent matrix used to determine a model classification ability not happing by chance [48]. It measures the level of similarity between observed and modelled distributions. The TSS index ranges from −1 to 1, where ≤0 indicates random agreement and values approaching 1 indicate a strong agreement between predicted and observed distribution [48]. To further reduce the baseline species pool to 687 species, we only considered species models with Maxent TSS ≥ 0.4. We chose TSS over AUC because the latter largely depends on the prevalence of sample locations and tends to be biased [47]. Subsequently, we used these 687 representative baseline species models to predict their distribution under current and future environmental conditions. From the 687 species, we chose eight common species in the Czech Republic [1,[32][33][34]49,50] to compare and assess their performance using Random Forest and GLM Logistic regression. Where necessary, we included the subspecies of the selected species among the 687 species to avoid ambiguity and ensure sufficient observations for each species in both monitoring periods (1960-1990 and 1991-2018). Monotypic species included spruce (Picea abie) and beech (Fagus syvestica L.). In contrast, those with two to 11 subspecies included alder (henceforth as Alnus sp.), fescue (henceforth as Festuca sp.), lipnice (henceforth as Poa sp.), oak (Quercus sp.), blackberry (Rubus sp.) and willow as (Salix sp.).

Mean Stable Area Index (MSAi)
The MSAi quantifies changes in habitat extent and indirectly measures conservation status. It was computed in two steps: first, we reclassified each species distribution map into binary maps based on each species-specific threshold for each modelling period. Second, for each modelling period, we compared (ratios) each binary species map with its corresponding baseline map in terms of habitat extent and calculated the mean of the comparison. We used the simulated species distribution from the baseline species data and not current empirical species data to ensure consistency and comparability. Hence, MSAi was computed considering changes in species composition and habitat extent as Equation (1).
where A i1,y1 is the area of species i for the baseline period and A i1 , y2 is the area of species i for a later modelling period. N is the total number of species that should be the same for the two modelling periods, irrespective of whether some species have disappeared in the future.

Climate Indices
The range of key climate indices for the different modelling periods is summarised in Table 2. The table shows a general increase in the average values of each climate index by 2100. Although the minimum and maximum annual rainfall (rain) may vary irregularly, the average value is projected to increase by~154 mm. Likewise, the average minimum temperature may increase by 5.2 • C by the end of the century. The average temperature of the growing season (tempgs) and the length of the vegetation periods are also projected to increase by 2.3 • C and 53 days, respectively. Compared to the current situation (2018), the table shows that the average minimum temperature has increased by 0.15 • C, the temperature of the growing season has increased by 0.85 • C while the average length of the growing season and annual rainfall have increased by~3 days and 132 mm respectively.

Model Performance
A summary of the model evaluation statistic for the initial 1177 baseline species modelled is presented in Table 3. The AUC varies from 0.43 to 0.99 and shows that fewer species models were randomly classified, while the majority was moderate to very good. The model TSS varies from 0.026 to 8.3 and shows that the correlation between observed and predicted species distribution for 50% of the species modelled was greater than 0.43. The performance measures of the eight selected species are presented in Table 4. The table shows an overall improvement in AUC based on the Random Forest classifier compared to Maxent and the GLM logistic regression. The Random Forest model for all selected species also estimated a moderate TSS comparable with the Maxent model. However, Maxent consistently predicted a higher TSS and showed a better match at the observed species locations than Random Forest or GLM models. Notwithstanding this difference, model performance was also consistent among modelling methods for some species. For example, the AUCs for the species beech, spruce and willow were almost consistent among the classification methods. While the AUC based on the GLM model was generally comparable to Maxent, it was the least accurate as it consistently underestimated the TSS value for all species.

Species Richness and Distribution
The estimated species richness for the baseline period is presented in Figure 2. Species richness is skewed to the right as expected and varies from 0 to~580. The mean species richness of the initial 11,117 species selected is~26. The modelled species richness from the pool of 687 baseline species is presented in Figure 3. The baseline species richness varies from 1 to 576, with about 80% of the landscape having 1 to 200 species (Figure 3a). Areas with high species diversity are mainly in the highlands of the northeast and some isolated plains. Compared to the current situation (2018), species richness varies from 39 to 531. About 70% of species richness is below the mean value (Figure 3b). The figure shows a spatial shift in species richness toward the highlands at borders. Figure 3c shows a drastic reduction in habitat extent and species richness by 2060, which will continue until 2100 (Figure 3c,d). Unlike the current and baseline situation with a few suitable areas on plains, the most favourable areas in the future will be restricted to highlands. The mean species richness is expected to drop from~200 species to~84 species by 2100.

Species Richness and Distribution
The estimated species richness for the baseline period is presented in Figure 2. Species richness is skewed to the right as expected and varies from 0 to ~580. The mean species richness of the initial 11,117 species selected is ~26. The modelled species richness from the pool of 687 baseline species is presented in Figure 3. The baseline species richness varies from 1 to 576, with about 80% of the landscape having 1 to 200 species (Figure 3a). Areas with high species diversity are mainly in the highlands of the northeast and some isolated plains. Compared to the current situation (2018), species richness varies from 39 to 531. About 70% of species richness is below the mean value (Figure 3b). The figure shows a spatial shift in species richness toward the highlands at borders. Figure 3c shows a drastic reduction in habitat extent and species richness by 2060, which will continue until 2100 (Figure 3c,d). Unlike the current and baseline situation with a few suitable areas on plains, the most favourable areas in the future will be restricted to highlands. The mean species richness is expected to drop from ~200 species to ~84 species by 2100. A general estimate of the simultaneous change in habitat extent (area) and species composition, summarised as MSAi for each modelling period for the entire landscape, is presented in Table 5. The estimated MSAi calculated as the product of the change in species richness and mean habitat extent between the baseline and each modelling period is projected to decrease by the end of the century. However, the current MSAi is the same as in the baseline period. The table also shows that ~11 out of the 686 (about 2%) of the baseline species may have been lost by 2018. About 25 species (~4.1%) of the baseline species may be lost by 2060. About 140 species (~20%) of the baseline species will become   A general estimate of the simultaneous change in habitat extent (area) and species composition, summarised as MSAi for each modelling period for the entire landscape, is presented in Table 5. The estimated MSAi calculated as the product of the change in species richness and mean habitat extent between the baseline and each modelling period is projected to decrease by the end of the century. However, the current MSAi is the same as in the baseline period. The table also shows that~11 out of the 686 (about 2%) of the baseline species may have been lost by 2018. About 25 species (~4.1%) of the baseline species may be lost by 2060. About 140 species (~20%) of the baseline species will become extinct by 2100. While the mean habitat extent expanded by~7.0% in 2018, it is projected to decrease to~49 and 46% by 2060 and 2100, respectively. The variability in the mean stable area (MSAi) values for the different modelling periods is shown in Figure 4. In Figure 4a, near-natural to natural areas with MSAi values close to 1 are currently confined to highlands around the borders and a few patches of flat areas inland. However, the extent of these areas are projected to rapidly decrease in the future (Figure 4b,c). extinct by 2100. While the mean habitat extent expanded by ~7.0% in 2018, it is projected to decrease to ~49 and 46% by 2060 and 2100, respectively.
The variability in the mean stable area (MSAi) values for the different modelling periods is shown in Figure 4. In Figure 4a, near-natural to natural areas with MSAi values close to 1 are currently confined to highlands around the borders and a few patches of flat areas inland. However, the extent of these areas are projected to rapidly decrease in the future (Figure 4b,c).

Species Respond to Climate and Environmental Conditions
The species selected to assess the role of climate, topography and environmental conditions showed diverse and nonlinear response patterns as expected. The four main variables affecting the distribution of all species included annual rainfall (anrain), minimum temperature (tcold), the temperature of the growing season (tempgs) and slope. The length of the vegetation period (lenvegt), the type, drainage, and to a lesser extent, the type of geologic material were equally important. The response patterns for the top four most important predictors, explaining at least 70% variability in distribution, are shown in Figures 5 and 6.
The species selected to assess the role of climate, topography and environmental conditions showed diverse and nonlinear response patterns as expected. The four main variables affecting the distribution of all species included annual rainfall (anrain), minimum temperature (tcold), the temperature of the growing season (tempgs) and slope. The length of the vegetation period (lenvegt), the type, drainage, and to a lesser extent, the type of geologic material were equally important. The response patterns for the top four most important predictors, explaining at least 70% variability in distribution, are shown in Figures 5 and 6. Figure 5. Species response to climate and environmental gradient: panel 1: alder; panel 2: beech; panel 3: fescue; panel 4: spruce; "anrain" = annual rainfall; "tcold" = average minimum temperature; "tempgs" + temperature of the growing season; "lenvgt" = length of the vegetation period and "geomat" = geological material. (a-d) annual rainfall, (e,o) average minimum temperature, (f-h) variation with slope angle, (i) average temperature of the growing season, (j-l) length of the vegetation season in days, (m) geological material scaled from 0 for very coarse to 1.9 for very fine material, (n,p) soil drainage scaled from 1 for poorly drained to 5 for well-drained soils Figure 5. Species response to climate and environmental gradient: panel 1: alder; panel 2: beech; panel 3: fescue; panel 4: spruce; "anrain" = annual rainfall; "tcold" = average minimum temperature; "tempgs" = temperature of the growing season; "lenvgt" = length of the vegetation period and "geomat" = geological material. (a-d) annual rainfall, (e,o) average minimum temperature, (f-h) variation with slope angle, (i) average temperature of the growing season, (j-l) length of the vegetation season in days, (m) geological material scaled from 0 for very coarse to 1.9 for very fine material, (n,p) soil drainage scaled from 1 for poorly drained to 5 for well-drained soils.
Alder (Alnus sp.) responded positively to precipitation. Their distribution was optimal in areas with annual total precipitation of 900 mm. Their distribution was optimal on gentle slopes (5 degrees) and generally facing the north or east. They also showed a near-stable response to minimum temperature (tcold) up to −5 • C, followed by a strong negative response with rising temperatures (Figure 5b,e,i).
Beech (Fagus syvestica L.) response was similar to alder. Their distribution was optimal on gentle slopes (10 to 25 degrees) and facing the north. They also preferred well-drained soils because their distribution expanded with distance from water bodies. Beech distribution increased with the length of the vegetation period (Figure 5b,f,j,n).
Fescue (Festuca sp.) The distribution of fescue generally decreased with precipitation and slope angle. However, their distribution increased as minimum temperatures rose above −4.5 • C. They tend to be common on gentle and south-facing slopes and showed a strong negative response with distance to water bodies. Their distribution was optimal in areas where the duration of the vegetation period was up to 220 days (Figure 5c,g,k,o).
Spruce (Picea abie) showed a positive response to precipitation and distance to water bodies similar to alder and beech and a similar response pattern to minimum temperature and slope variation as beech. They preferred well-drained soils typically on moderate slopes facing north to northeast and will grow optimally at 12 • C like beech species (Figure 5d,h,l,p). "anrain" = annual rainfall; "tcold" = average minimum tempera-ture; "tempgs" + temperature of the growing season; "lenvgt" = length of the vegetation period; and "geomat" = geological material. (a-d) annual rainfall, (e,o) average minimum temper-ature, (e) geological material scaled from 0 for very coarse to 1.9 for very fine material, (f-h) slope angle, (i,j,l) average minimum temperature, (k,n) soil drainage scaled from 1 for poorly drained to 5 for well-drained soils, (m,o,p) average temperature of the growing season.
Alder (Alnus sp.) responded positively to precipitation. Their distribution was optimal in areas with annual total precipitation of 900 mm. Their distribution was optimal on gentle slopes (5 degrees) and generally facing the north or east. They also showed a nearstable response to minimum temperature (tcold) up to −5 °C, followed by a strong negative response with rising temperatures (Figure 5b,e,i).
Beech (Fagus syvestica L.) response was similar to alder. Their distribution was optimal on gentle slopes (10 to 25 degrees) and facing the north. They also preferred welldrained soils because their distribution expanded with distance from water bodies. Beech distribution increased with the length of the vegetation period (Figure 5b,f,j,n) Fescue (Festuca sp.) The distribution of fescue generally decreased with precipitation and slope angle. However, their distribution increased as minimum temperatures rose above −4.5 °C. They tend to be common on gentle and south-facing slopes and showed a strong negative response with distance to water bodies. Their distribution was optimal in areas where the duration of the vegetation period was up to 220 days (Figure 5c,g,k,o) Spruce (Picea abie) showed a positive response to precipitation and distance to water bodies similar to alder and beech and a similar response pattern to minimum temperature and slope variation as beech. They preferred well-drained soils typically on moderate slopes facing north to northeast and will grow optimally at 12 °C like beech species (Figure  5d,h,l,p).
Lipnice (Poa sp.). The distribution of lipnice gradually increased with precipitation Lipnice (Poa sp.). The distribution of lipnice gradually increased with precipitation above 750 mm. They showed a stable response to minimum temperature variations but an irregularly positive response above −4.5 • C. They became common above 12 • C and were optimal on south-facing slopes. They also preferred well-drained soils (Figure 6a,e,i,m).
Oak (Quercus sp.) responded negatively to annual rainfall. However, their distribution increased with slope values and distance to water bodies, as minimum temperatures warmed above −4 • C. They also show a bimodal positive response to changes in temperature of the growing season (Figure 6b,f,j,n).
Blackberry (Rubus sp.) distribution increased with annual rainfall and was optimal at 1000 mm. Their distribution also increased with the temperature of the growing season and optimal on gentle slopes having well-drained soils (Figure 6c,g,k,o).
Willow (Salix sp.). The impact of annual variations in precipitation on the distribution of willows was like that of oaks and generally decreased with precipitation. Their distribution also decreased with slope values and annual rainfall. They may be optimally distributed on southeast-facing slopes where annual rainfall and slope values are 650 mm and 10 • , respectively. (Figure 6d,h,l,p).
The selected species' habitat extents for the current period are summarised in Table 6. In general, the average climate space of these species expanded by~12% between 1991 and 2018. Spruce and willow are the most adaptable, accounting for about 72% of net habitat expansion. The distribution of fescue and alder is already declining, while Lipnice and oak are equally at risk of losing their climate space.

Discussion
The impact of climate and environmental change on individual species distribution is very diverse but varies with the local topography. Rising temperatures have generally shifted species towards highlands (Figure 3b). Species richness has slightly declined under the current climate as more than 97 per cent of the representative baseline species are currently preserved in most areas ( Table 5). The change is due to the near stable climate between the two modelling periods ( Table 2). Our results show little change in climatic conditions. The average minimum temperature was nearly the same between these two periods. The mean temperature of the growing season increased by 0.85 • C, while the mean length of the vegetation period increased by three days. Although species richness is nearly the same, species habitat expanded remarkably between the two modelling periods as growth conditions have become more favourable for most species. Although these conditions have extended highland habitats where low temperature is a limiting factor for growth [51], further increases in temperature will be devastating to these habitats. Our results show that species composition and habitat extent will decline as the average minimum temperature and the growing season temperature rise by +5 • C and +3 • C, respectively ( Table 2). More than 20% of the baseline species may be at risk of becoming extinct at the end of the 21st century (Table 5). These results are comparable to [1,33], who predicted frequent heat spells in lowland habitats based on a moderately mitigated climate scenario. As growth conditions under the baseline climate scenario may become too extreme for most species, these results should be interpreted with caution because they are only a simulation of what may be possible [39,52,53]. They do not consider global or country-specific efforts such as substituting coal with clean and renewable energies to reduce greenhouse gas emissions [53]. Hence an objective assessment based on a mild or moderate climate scenario will be necessary for future studies.
The variability in MSAi values has reaffirmed that the most stable areas of the Czech Republic are currently in protected and mountainous areas (Figure 5f,g,h). Their current MSAi values range from 0.7 to 0.94 and may drop from 0.5 to 0.8 by 2100 without intervention or mitigation efforts. Lowlands with the least species variety are the least stable and the most vulnerable. Our results show more variability in the MSAi ratio for the Czech Republic than the regional EUROMOVE model for Europe [2,5]. A possible reason for the difference could be that we modelled change based on 686 species for the Czech Republic compared to 430 species for the entire Czech Republic, Slovakia, and Hungary in the regional EUROMOVE model [2]. The extra details also highlight the benefits of using high-resolution climate and environmental data to account for local variations [6,54].
The advantage of quantifying change in terms of MSAi is that additional information about the state of the landscape, which is more related to ecosystem functions than species richness alone, is known. Experimental studies have generally associated a decline in species richness with a decline in biomass production, leading to a 20 per cent loss in species as a proposed threshold for stable ecosystems [22]. The application of such speciesbased thresholds in nature has been questioned due to inconsistencies in the underlying processes that affect species richness [55]. Our results also show that species loss may not be proportionate to the loss of potential habitat extent. (Table 5). Second, the loss of a few dominant species may drastically shrink or expand habitats, impacting selected ecosystem function and services. Thus, integrating both parameters to obtain information about the state of landscapes, we expect vulnerability thresholds established from MSAi to be more reliable and applicable than those based solely on species richness. While assessing ecosystem function is not our focus, our result also shows that MSAi may be used as a validation tool or dataset to supplement such studies because changes in stable areas are based on surveyed records. Stable areas can be compared to favourable areas or classes of land use and land cover data preserved or appearing over time as criteria for assessing ecosystem function and services [56,57]. Therefore, the detailed spatial variation in MSAi has highlighted highly vulnerable areas where a decline in species richness relative to habitat extent should be accompanied by a loss of key ecosystem functions and services.
The response of the eight selected species justifies grouping species with nearly the same ecological requirements into distinctive FVZs as an effective management option for biodiversity in the Czech Republic. Generally, species most tolerant of high precipitation, including alder, beech and spruce, become more adaptable as minimum temperature decreases [1,33]. These species will typically prefer well-drained soils and are expected to thrive on moderate to steep slopes. In contrast, species tolerant of low and moderate precipitation, including fescue, lipnice, oak, blackberry, and willow, prefer gentle slopes where soil moisture is high and soil drainage is low to moderate. Because these major species impact the distribution of understory species, current efforts should focus on their preservation. Appropriate species mixing in the different FVZs could be adapted as a long-term strategy to buffer the most vulnerable species and minimise further species loss [58]. The selected species' response also reflects how habitats have shifted within the main forest vegetation zones of the Czech Republic from 1990 to 2018. The results in Table 6 show that habitat contraction has occurred primarily in the first, second and third FVZs, typical of lowlands and dominated by alder, oak, and fescue species. Habitat contraction has been accompanied by a wider expansion in the higher (sixth to eighth) FVZ where spruce and willow species dominate. These trends are consistent with the works of [1,32], who attributed the shift to rising temperatures and decreasing precipitation in the lower FVZs.
The main limitations of this study are that while MSAi broadly captures the impact of climate change on species richness and landscape vulnerability, it does not adequately quantify ecosystem vulnerability. The impact of climate change on key biomes or individual FVZs of the Czech Republic has not been assessed to understand habitat vulnerability and biodiversity loss. Because resistance to climate change varies across ecosystems, we hope to address this limitation in a future study by adapting the GLOBIO modelling framework to understand how different drivers of change related to land use (agriculture and infrastructure) will affect biodiversity [25,26]. Second, we did not evaluate change under a mild or moderate climate scenario because such data were not available at the appropriate resolution when this study was carried out. Future studies should be based on recently downscaled climate models within the coupled model intercomparison project (CMIP6) framework, which are more likely to show variability in climate indices on a smaller scale than CMIP5 models on which our results are based [59]. Third, EUROMOVE and its MSAi indicator heavily depend on species data quality. They may produce inconsistent and incomparable results for species records acquired over different periods. We have observed (results not presented) a rather weak correlation between the species richness map for 2018 projected from 1990 SDM models compared with that directly modelled from the 2018 species records. Although the discrepancy may have been due to slight differences in species composition between the two modelling periods, it also suggests possible uncertainties in the MSAi. Uncertainties may also be introduced because the interaction we assumed to be present among species both in the present and in the future may not be the case. In addition, species richness was probably overestimated because we did not model rare species and the level or type of association among common species [16,60]. Lastly, although our results reflect the response of key and representative species to climate and environmental change, they are based on species richness rather than species abundance or both. As such, we acknowledge the fact that information about population trends between species and ecosystem stability has been lost or not adequately captured [30].

Conclusions
The richness of plant species in the Czech Republic is slowly declining in response to climate change. Although the current species habitat has expanded compared to the baseline, species typical of lowland habitats are the most vulnerable. The decline in species richness may be drastic in the future without any mitigation efforts. Our assessment of the entire landscape stability using the mean stable indicator (MSAi), adapted from the EUROMOVE model, has reaffirmed that the most stable areas are currently restricted to highlands. However, their stability will rapidly decline if local temperatures rise to extreme levels. Thus, there is a need to objectively quantify change based on a mild or moderate climate scenario, given that such extreme conditions are unlikely with the current efforts to curb greenhouse gases. The MSAi does not specifically quantify the stability of key biomes and forest vegetation zones (FVZs) of the Czech Republic. However, it is a comprehensive tool that broadly summarises the impact of climate change on habitat extent and species composition into an indicator of area stability that can also be considered an indicator of landscape vulnerability. It may be used to validate habitat-based models in assessing ecosystem function and services as survey data become available.
Funding: This research was supported by the Internal Grant Agency of Palacky University, Olomouc (IGA_PrF_2022_027) as part of the project "Advanced application of geospatial technologies for spatial analysis, modelling, and visualisation of the phenomena of the real world".
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Original datasets are not available to the public but can be obtained by contacting the corresponding institution listed in Table 1.