Topography of the Dolomites modulates range dynamics of narrow endemic plants under climate change

Climate change is expected to threaten endemic plants in the Alps. In this context, the factors that may modulate species responses are rarely investigated at a local scale. We analyzed eight alpine narrow endemics of the Dolomites (southeastern Alps) under different predicted climate change scenarios at fine spatial resolutions. We tested possible differences in elevation, topographic heterogeneity and velocity of climate change among areas of gained, lost, or stable climatic habitat. The negative impact of climate change ranged from moderate to severe, depending on scenario and species. Generally, range loss occurred at the lowest elevations, while gained and stable areas were located at highest elevations. For six of the species, climate change velocity had higher values in stable and gained areas than in lost ones. Our findings support the role of topographic heterogeneity in maintaining climatic microrefugia, however, the peculiar topography of the Dolomites, characterized by high altitude plateaus, resulted in high climate change velocity in areas of projected future climatic suitability. Our study supports the usefulness of multiple predictors of spatio-temporal range dynamics for regional climate-adapted management and eventual assisted colonization planning to not overlook or overestimate the potential impact of climate change locally.

www.nature.com/scientificreports/ dynamics 19 . Consequently, a complex topography may favor the presence of microrefugia due to increased habitat diversity and niche space availability 20,21 . Microrefugia have increased the likelihood of species persistence during periods of environmental alterations (e.g., during past climatic fluctuations 22,23 ). For these reasons, topographically complex alpine areas are considered a much 'safer' place to live under climate change than flat areas, which offer no short-distance escapes from new climatic conditions 24 .
Even if mountain topography buffered the effects of climate change on plants 25 , the poor ability to disperse of endemic plants 26 can impede them from tracing the geographical shift in climatically suitable environments. For this reason, areas with a low rate of climate change may have favored 27 , and likely will favor, the persistence of endemic species 28 .
The Dolomites are one of the most relevant areas for plant biodiversity in the Alps 29 , being recognized as one of the richest centers of endemism and one of the main areas of glacial refugia 30,31 . Moreover, these mountains have unique geological and resulting topographical structures 32,33 and are one of the most protected regions in the Alps, with ten natural parks and nine systems according to UNESCO within the World Heritage Sites 34 . In these areas, vegetation changes (i.e. increasing species numbers at high elevations and 'thermophilization') are currently taking place in the alpine zones 35,36 . However, as the Dolomites likely acted as refugia for cold adapted species during the last glaciations ensuring climatically stable areas, they might similarly enable species to survive future climate change in local refugia 28 .
Species distribution models (SDMs 37 ) allow to predict the shift in suitable climatic conditions of species under environmental change 38 . However, because coarse resolution climate modeling approaches do not account for the importance of microclimatic conditions, they may underestimate the habitat suitability and overestimate the species' rate of extinction 39 . For such reason, fine-scale modeling based on climatic downscaling techniques can be more effective in considering such microhabitat conditions 40 .
Using species distribution models, we analyzed the potential effects of climate change on eight plants endemic to the Dolomites under different climate change scenarios. Moreover, we tested whether areas where species are projected to lose, gain or maintain their suitable habitat differ in elevation, topographic heterogeneity and in velocity of climate change. First, we expected an extensive reduction of the distributional range of the endemics in the Dolomites, little balanced by gain of newly suitable areas. Secondly, we expected that endemics lose climatically suitable habitats at low elevation, in areas with low topographic heterogeneity and high climate change velocity, and that they maintain or gain suitable habitat at high elevation, in areas with high topographic heterogeneity and a low climate change velocity.

Methods
Study area and selected species. The study area includes the Dolomites and surrounding areas with a total area of about 8325 km 2 . In the last centuries, the field surveys carried out by several botanists provided an extensive documentation on the distribution of endemic plants 29,[41][42][43] . In this study, we selected eight out of 25 endemic species with a restricted range present in the Dolomites 43 (Table 1,  The species were selected based on their sub-alpine/alpine distribution and on the availability of more than 30 points of occurrence, which has been shown to provide reliable results in modeling studies 44 . Moreover, these species represent the taxonomic and ecological diversity of endemics in the alpine zone of the Dolomites belonging to eight families and growing in different habitats (i.e. screes, rocks, ridges and alpine grasslands).
All eight taxa are of conservation concern 45-48 due to their narrow distribution range, rarity and high habitat specificity (Table 1). We obtained occurrences from public regional floristic databases and expert observations www.nature.com/scientificreports/ of Barbet-Massin et al. 56 (Supplementary Table S2). For each pseudo-absence set, a split-sample cross-validation was repeated 10 times, using a random subset (30%) of the initial data set. We assessed model performance by using the area under the receiver operating characteristic curve (AUC) 57 and the true skill statistic (TSS) 58 . As the choice of threshold may affect projection bias, we used three different thresholds performing equally or better than others 59,60 to convert the continuous suitability maps into binary projections of species presence and absence: the threshold selection method based on (i) equal training sensitivity and specificity, (ii) maximizing training sensitivity and specificity, and (iii) minimizing the distance between the curve and the upper left corner of ROC plot. Finally, we generated a total of 165 binarized projections for each species: 15 for baseline climate (3 threshold * 5 algorithms) and 75 for each future climate rcp scenario (3 threshold * 5 algorithm * 5 GCMs). Then, we applied the majority consensus rule among the binarized projections: we considered a species as occurring in a cell if at least 50% of all the models predicted its occurrence there. Further, as the study species grow only on limestones and dolomitic rocks 42 , we used a geological layer obtained rasterizing at 50 m the geolithological map of Italy downloaded from the Geoportale Nazionale (http:// www. pcn. minam biente. it/ mattm/) to mask SDMs outputs (Supplementary Figure S1). All analyses were carried out on the Vienna Scientific Cluster (VSC) by using R (version R/3.6.2) 61 . The full-dispersal scenario is likely to be unrealistic because it assumes that a species can colonize all locations without physiological, environmental or geographical limitations, while the no-dispersal scenario is likely to be more appropriate for poor dispersers as endemic species frequently are 26 (Table 1).

Range analysis and maps.
Climatic habitat-suitability relationship with elevation, topographic heterogeneity and climate change velocity. We explored the relationship among change in habitat suitability (loss or stability in the future) and elevation, topographic heterogeneity and climate change velocity (CCV). Elevation was assessed from the DEM at 50 m resolution (see above). The topographic heterogeneity was calculated with two different indices: (i) the Terrain Ruggedness Index (TRI) 62 with "tri()" function of "spatialEco" package in R that estimates for the among-cells grid complexity and is calculated as the mean of the absolute differences between the value of a cell and the value of its eight surrounding cells, (ii) the intrinsic topographic complexity index (TCI) to account for within-cells grid complexity (following Irl et al. 63 ), calculating the ratio between the 3D and 2D surface area, including the Tinitaly DEM at 10 m resolution (from http:// tinit aly. pi. ingv. it/) 64 , with the following equation: The multivariate climate change velocity surface was calculated using the protocol described in Hamann et al. 65 . In short, climate change velocity was calculated as the minimum geographic distance to the nearest analogous climate, divided by the number of years between the baseline climate period and the future projection for both the intermediate (rcp 4.5) and the realistic (rcp 8.5) climate scenario. The CCV value is calculated as log 10 of velocity (m*yr -1 ) multiplied by 100. We used the first two PCA axes calculated from the 19 bioclimatic variables and 120 bins of unique climates to define climate matches for each of the 5 GCMs, then we averaged them to obtain two scenarios of CCV (rcp 4.5 and rcp 8.5).
We assessed the correlation among the five variables (elevation, TCI, TRI, CCV rcp4.5, CCV rcp8.5) over the whole study area using the Pearson correlation coefficient implemented in R.
We created a matrix of each loss, gain and stable spatio-temporal category of difference of climatic habitat suitability between the baseline and the future predicted conditions for each species. We analyzed the relationship among habitat suitability and elevation, topographic heterogeneity and climate change velocity for each species using the Kruskal-Wallis test, with 0.05 significance. Subsequently, we tested the difference of elevation, topographic heterogeneity and climate change velocity among areas of climatic loss, stability or gain by using the post-hoc Mann-Whitney test (pairwaise.wilcox.test() R function) with 0.05 significance and "Bonferroni" correction 66 for each species. To test for the robustness of these results and for pseudo replication issues, we run the analysis randomly sampling 10% and 1% of cells of each range change category. We repeated the randomization procedure 1000 times, then we counted how many times the Kruskal-Wallis test was significant (p value < = 0.05). Second, we performed the previous procedure setting the maximum number of cells of each range change category to 1000.

Results
Overall patterns of change in habitat suitability. Under baseline climatic conditions, model evaluation showed a good model performance for the majority of modelling techniques in all analysed species (Supplementary Table S3). Considering each species' response, range loss was different between the two future sce- www.nature.com/scientificreports/ narios. In particular, the percentage of loss ranged from 51 to 79% under the rcp 4.5 scenario, while it was higher under the rcp 8.5 scenario ranging from 72 to 92% (Table 2, Table S4). The habitat loss mainly occurred at lower elevations and at the geographical periphery of distributional range of each species (Supplementary Figures S3-S17). The percentage of range gain was low (< 10%) for C. The correlation among changes in climatic habitat suitability and elevation, topographic heterogeneity and climate change velocity. Elevation was correlated by 18% and 19% respectively with TCI and TRI and by 26% with the intermediate CCV, while by 49% with the realistic CCV. TCI and TRI showed a correlation of 65%. A low correlation was found for the topographic heterogeneity parameters with CCV, with a little increase in the realistic scenario (Table 3).
For each species, each category of habitat suitability (loss, gain and stability) resulted always significantly different among the different indexes (Tables S5-S6). Similar results were obtained by using randomization, in fact the number of times that p values were significant was very high with a minimum of 926 times. Nevertheless, the pairwise test detected the following exceptions that were not significant (Table S7): stable areas compared to and lost areas for TCI in F. austrodolomitica for rcp 8.5; stable areas compared to lost and gained areas for both TCI and TRI in N. buschmanniae; gained areas compared to lost and stable areas for TCI and TRI for rcp 4.5 in R. alpina and S. facchinii. In the majority of the species (except for S. facchinii), climatically stable areas were expected to occur at high elevation. Also, together with C. morettiana and R. alpina, S. facchinii showed lower values of elevation for the areas of gain in the rcp 4.5 scenario, due to the low extension of gain areas (Table S4). In the rcp 8.5 scenario, the elevational trend was comparable but magnified for all species, compared to the rcp 4.5 scenario for at least the loss/stable comparison, with losses and stability shifted at higher elevations (Figs. 2, 3, Table S7).
Under the rcp 4.5 scenario, TCI had only few times high values in stable areas for six species, while N. buschmanniae and G. brenate had only few times higher values in lost areas (Fig. 2). The gained areas had low topographic complexity in all the species (except S. dolomiticum that had no gain areas). A similar response was detected under the rcp 8.5 scenario but differences were not significant in F. austrodolomitica and N. buschmanniae. On the contrary, G. brentae showed significantly lower TCI values in stable areas (Fig. 3,  www.nature.com/scientificreports/   Table S6). The trend of terrain ruggedness (TRI) was similar to that of TCI for four species (C. morettiana, F. austrodolomitica, N. buschmanniae, P. tyrolensis), while for G. brentae the TRI was only few higher in areas of loss compared to both stable and gained areas (Fig. 3, Supplementary Table S7). Moreover, R. alpina did not show significantly higher values in areas of gain. A similar pattern was detected under the rcp 8.5 scenario (Fig. 3, Supplementary Table S7). Under the rcp 4.5 scenario, six out of eight species (C. morettiana, F. austrodolomitica, G. brentae, N. buschmanniae, P. tyrolensis and S. dolomiticum) had lower values of CCV in lost areas than in stable and gained areas. Differently, in R. alpina and S. facchinii the lowest values of CCV were recorded in gained areas. In particular, S. facchinii was the only species having the highest CCV in the lost areas. Similarly, under the rcp 8.5 scenario, significantly lower values of CCV were recorded in areas expected to be lost by S. facchinii. In N. buschmanniae, areas projected to be lost and stable had low but not significantly different values of CCV. Moreover, among the species projected to gain range, F. austrodolomitica, G. brentae and N. buschmanniae showed highest CCV in the areas of gain, while P. tyrolensis had lower values of CCV in areas projected to be gained.

Discussion
This work is the first study that used high resolution SDMs to predict the impact of future climate change in a center of endemism in the Alps (the Dolomites, Italy). Our fine resolution modeling approach is relevant to project the risks of biodiversity loss in mountain systems, where marked microclimatic variation can allow species to persist locally. Further, it is well-suited for the assessment of species-specific climate change threats in particular regions. In line with previous findings in similar parts of the Alps and in other continents 8,9,14,[67][68][69] , our results suggest a likely extensive species-specific reduction of the distributional range of the endemics, little balanced by gain of newly suitable areas. However, the peculiar and complex topography of the Dolomites, where the slope is sharpest at middle elevation, while the declivity decreases forming highlands near the top of the mountains 32,33 , resulted in higher climate change velocity at high elevations and low climate change velocity in middle altitude slopes, where species are projected to lose their suitable habitats.
We projected a reduction in the ranges of high-altitude endemics of the Dolomites, despite the range contraction strongly depended on the climatic rcp scenario. Despite the fact that all the species were predicted to be negatively affected by climate change, for the subalpine species occurring in the southernmost periphery of the Dolomites (i.e. P. tyrolensis), the range loss might be partially buffered by a certain amount of potential gain towards north and towards high altitude. Similarly, high range gain values were detected in the alpine species G. brentae and N. buschmanniae. Most of the alpine endemics analyzed here occupied only a small part of climatically suitable areas suggesting that post-glacial range-shift limitations may have resulted in a lack of range filling, due to barrier effects and dispersal limitations, as previously detected for other alpine plants 70 . Nevertheless, other factors such as competition, soil conditions and nutrients' availability could have prevented the range filling 5,71 . However, for these species, a large part of the study area was climatically suitable, and consequently the areas of gain were quantitatively small. Therefore, for these narrow endemic species, the few currently unsuitable areas that became suitable under future climate, were mainly proximal to the currently suitable areas. In contrast, areas currently suitable for G. brentae and N. buschmanniae occured principally in areas that were occupied extensively by the species. Therefore, considering the poor dispersal capability of endemic species, only P. tyrolensis may likely reach the new suitable areas. Differently, the other species that showed a range gain (i.e. G. brentae and N. buschmanniae), occuring almost only in the Brenta group 41,42 , could hardly keep pace with climate shifts due to the high distance between the current occurrences and the new suitable areas. In line with the hypothesis of a taxon-specific response to climate change, previously detected in Mediterranean mountains 72 , this speciesspecific response suggested that endemics may respond differently to future climate change, as detected for their response to past climate changes in the Alps 17,73 .
Furthermore, the current global change is predicted to be more extensive than the past changes 4 . The latitudinal and altitudinal shift, which we projected in some endemics, is in line with previous findings on plants both locally and globally 74,75 . Nonetheless, latitudinal and altitudinal shifts to suitable areas may differently affect species survival 76,77 . Indeed, the geographical distance between different climatic zones is shorter along altitudinal than along latitudinal gradients 78 . Consequently, endemics, despite their poor dispersal abilities 79,80 , may more likely keep pace with the upward shift than with the latitudinal one, since lower distances are required to follow the climate shift upward. The eight studied narrow endemic species of the Dolomites have a low dispersal ability (boleochory and blastochory) limited by both the absence of specialized diaspores for dispersal and by the short stem height (Table 1) 26 . Thus, within the time interval of the predicted future climate projections (75 years), all the studied species (except F. austrodolomitica) could probably move only a hundred meters away from their current presence (i.e. dispersal distance of about one meter per year), which is very close to a no-dispersal scenario.
According to the general expectation of an upward range shift 35,81 and the local high-elevation persistence hypothesis 16 , the majority of endemic species of the Dolomites were projected to lose their range at lowest elevations and to gain (few) or maintain range at highest elevations under both scenarios. Nevertheless, high-altitude endemics like S. facchinii, which currently grows above 2600 m a.s.l., and N. buschmanniae, with a current range between 2100 to 2400 m a.s.l. only in the Brenta group, could undergo a "no-where to go" scenario only under the rcp 8.5 scenario 82 (Figures S9-S15).
In the Dolomites, geomorphologically heterogeneous environments were only slightly more complex in areas climatically stable for endemics than in areas where climatic suitability will be lost. Previous studies 18,21,83 suggested that local geomorphological heterogeneity may increase the availability of microclimatic refugia where alpine species may move or persist 15 . Topographic heterogeneity may cause strong temperatures differences over short distances 15 and may decouple the local climate from the surrounding landscape buffering the effect of extreme temperatures and favoring survival of relict and endemic species 84 . In endemic centers such as the www.nature.com/scientificreports/ Dolomites, the rugged topography and resulting microclimatic diversity are thought to have locally buffered the endemics from extinction during past climate change 27,85 . The lower topographic heterogeneity in stable areas detected for G. brentae and N. buschmanniae is likely explained by the fact that these species were expected to lose all climatically suitable areas within their distributional range and to maintain habitat suitability in the future in northernmost high-altitude flat areas, where topographic complexity was low. The difference in topographic complexity was significant but low between areas of range loss and stability, suggesting that the availability of microclimatic refugia is similar between these areas. Therefore, differences in microrefugia will play a secondary role in assuring endemics survival in the Dolomites. In general, areas with a slow rate of climate change velocity were expected to provide important refugia for species survival under climate change 28,86 . Contrary to this general expectation, however, we detected that the velocity of climate change was lower in areas where species were projected to lose their suitable habitat than in areas where species were projected to retain or gain suitable habitat. Climate change velocity was inversely related to slope, in fact, the spatial gradient of climatic change was greatest on steep slopes where consequently modest shifts in space are required to meet similar climatic conditions 25 . Although counter-intuitive, the low climate change velocity at low altitudes may be explained by the peculiar topography of the Dolomites, where the slope is sharpest at middle elevation, while the declivity decreases forming highlands near the top of the mountains 87 . The populations of alpine endemics growing at the lowest elevation limit of the species were expected to occur near the species' warm climatic boundaries 67 . Therefore, they were more prone to extinction risk due to global warming, regardless the low climate change velocity. By contrast, the highlands, where climate change velocity was high, will probably harbor climatic conditions that currently prevail at middle elevations, becoming suitable for endemics in future. In fact, the high elevation areas in the northern and central Dolomites showed the highest values of climate change velocity, making them a difficult terrain for the survival of existing species but possibly also new territory for species projected to occur in these areas under future climate.

Study limitations.
Our results should be considered with particular attention to the main limitations of the approach we adopted. Apart from the use of a different set of GCMs that could provide different results 51 and the inclusion of missing occurrence records that could affect SDMs 88 , we identified three further main limitations. First, range loss and gain were calculated on binary models' predictions. A binary cutoff may reduce the predictive ability of models, therefore, such results should not be used to draw inferences on single species but they are useful to infer trends at the level of geographic areas 88 . Second, we modeled the species without any biotic interaction, which has been demonstrated to drive key ecological and evolutionary processes, mediating ecosystem responses to climate change 89 . Third, even if our spatial resolution (50 m) would account for microclimatic variability 90 , small topographical structures under 10 m of resolution like small depressions and crevices could be hardly captured, but they could affect microclimates.

Conclusive remarks.
Our results suggest that alpine narrow endemic plants will be at high risk under climate change in the Alps 9,14 . Topographic and microclimatic conditions may still provide extended areas of refuge only if interventions aimed at reducing and/or offsetting emissions would keep climate change within the rcp 4.5 scenario. Furthermore, the study of the relationship among multiple spatio-temporal predictors (e.g. CCV, TCI etc.) and the species range dynamics might be useful for regional climate-adapted management to avoid overlooking or overestimating the potential impact of climate change locally. The stable areas highlighted in this study are potential "refugia", areas that should be conservation priorities with the goal of maintaining current patterns of biodiversity. In case of predicted high range loss, the identification of future suitable areas, together with the evaluation of the likelihood of dispersal, are important tools for future conservation planning 91,92 , to maximize the conservation benefit in terms of range loss compensation for rare species at risk of local extinctions.

Data availability
We report detailed references of occurrence data in the 'Supporting information' file (Table S1). Data on species occurrences as centroid for the 50 m × 50 m cell grid are available at the Figshare data repository https:// doi. org/ 10. 6084/ m9. figsh are. 18401 966 = data will be made publicly available once the manuscript is accepted for publication. www.nature.com/scientificreports/