Shifting aspect or elevation? The climate change response of ectotherms in a complex mountain topography

Climate change is expected to cause mountain species to shift their ranges to higher elevations. Due to the decreasing amounts of habitats with increasing elevation, such shifts are likely to increase their extinction risk. Heterogeneous mountain topography, however, may reduce this risk by providing microclimatic conditions that can buffer macroclimatic warming or provide nearby refugia. As aspect strongly influences the local microclimate, we here assess whether shifts from warm south‐exposed aspects to cool north‐exposed aspects in response to climate change can compensate for an upward shift into cooler elevations.

However, mountains exhibit a complex topography and factors like aspect affect the local temperature and water balance and lead to heterogeneous microclimates at small spatial scales, which often differ starkly from the regional climate (Austin & van Niel, 2011 ;Dobrowski, 2011 ;Scherrer & Körner, 2011 ;Winkler et al., 2016 ).
These microclimates can locally buffer against the effects of regional warming, offer stepping stones for migration or provide suitable alternative habitats within short distances and hence weaken the impact of climate change and may even prevent an upslope shift (Kulonen et al., 2017 ;Meineri & Hylander, 2017 ;Opedal, Armbruster, & Graae, 2015 ;Scherrer & Körner, 2011 ;Suggitt et al., 2018 ). Thus, as an alternative scenario, mountain species may move to nearby locations with suitable microclimates (e.g., cool north-exposed locations) instead of migrating upward to avoid, for example rising temperatures.
In the Northern Hemisphere, south-exposed locations are associated with higher maximum temperatures, shorter duration of snow cover, longer warm seasons and more evaporation than north-exposed ones (Dobrowski, 2011 ). As the local preference of a certain aspect of a species can change with elevation (i.e. different regional climatic regimes), distributional shifts towards cooler north-exposed locations might be more likely at lower elevations (warmer temperature regime), while at higher elevations (cooler temperature regime) warmer south-exposed locations might be preferred (Dobrowski, 2011 ). Despite the fact that fine-scale microclimatic features may influence the distribution of species, analyses of climate change effects often have a coarse spatial scale. For instance, the spatial resolution of input data of most species distribution models (SDM) (generally 1 km or 10 km for larger extents) is orders of magnitude larger than the size of animals under study (Potter, Arthur Woods, & Pincebourde, 2013 ). Understanding how climate change impacts species at a fine scale-the scale at which microclimatic variability may buffer against the effects of climate change-is a critical step to take for climate change biology, as data availability is still limiting high-resolution analyses (Kearney et al., 2014 ;Lembrechts, Nijs, & Lenoir, 2019 ;Potter et al., 2013 ;Suggitt et al., 2018 ).
Here, we assessed the fine-scale impact of climate change on the distributions of two ectothermic vertebrate species. Using data from the Swiss Alps, we built SDM based on fine-scale environmental data and projected them into multiple future climate scenarios to examine whether the study species respond to climate change through upslope shifts or through shifts in aspect or both. Climate change can have particularly strong impacts on ectotherms, in which body temperature and physiology are heavily related to external temperature, especially if they are viviparous (Buckley, Hurlbert, & Jetz, 2012 ;Deutsch et al., 2008 ;Dillon, Wang, & Huey, 2010 ;Hoffmann, Chown, & Clusella-Trullas, 2013 ;Sinervo et al., 2010 ).
We focused on the Alpine salamander ( Salamandra atra ) and the Common lizard ( Zootoca vivipara ). The former is endemic to the Alps while the latter has a wide geographic range; both often co-occur in syntopy in the study area. As both species are viviparous, they are independent of special requirements for reproduction, such as the fine-scale distribution of water bodies in case of the Alpine salamander or specific egg deposition sites in the case of the Common lizard (Dely & Böhme, 1984 ;Guex & Grossenbacher, 2004 ). Such a dependency of specific breeding sites might limit the ability of the species to respond to climate change. The Alpine salamander is considered to be vulnerable to climate warming due to its restricted distribution in a climate change sensitive area, its preference of cool humid habitats and low dispersal ability (Rabitsch et al., 2010 ;Schlumprecht et al., 2010 ). Rising temperatures may also have negative impacts on populations of the Common lizard due to more restricted activity periods, declined juvenile dispersal or decreased adult survival (Bestion, Teyssier, Richard, Clobert, & Cote, 2015 ;Massot, Clobert, & Ferrière, 2008 ;Wang, Ma, Shao, & Ji, 2017  between 372 and 2,849 m a.s.l. (mean 1,330 m a.s.l.) which we used for analyses (Figure 1 ).

| Environmental data
We used fine-scale climate data for Switzerland with a spatial resolution of 100 m and a yearly temporal resolution. Climate maps were spatially interpolated from national "MeteoSchweiz" climate station data using Daymet (Thornton, Running, & White, 1997 ). It covered current conditions (1977-2006, hereafter "2000" period).
In addition, two future periods (2041-2070, hereafter "2050" period; 2071-2100, hereafter "2080" period) for three different climate models and two representative greenhouse gas concentration pathways (rcp 4.5 and 8.5) were calculated, resulting in 12 future scenarios (Appendix S2 : Table A1 ). Future projections were generated using model data from the EURO-CORDEX set of climate projections (Jacob et al., 2014 ) and using the delta change method (Anandhi et al., 2011 ). To this end, we downloaded EURO-11 data that is available at a ca. 12.5° spatial resolution, calculated anomalies at this coarse resolution between the current period and the two future periods, downscaled the anomalies to a 100-m spatial resolution (cf. Guisan, Thuiller, & Zimmermann, 2017 ), and added the anomalies to the current period climate data. The three regional climate models (RCMs) from which we used EURO-11 model runs from CORDEX were as follows: CCLM (Rockel, Will, & Hense, 2008 ), Hirham (Christensen, Christensen, Machenhauer, & Botzet, 1998 ) and RACMO (van Meijgaard et al., 2012 ). Of the yearly means of temperature and precipitation, we calculated 30-year averages and standard deviations (SD) for all periods, which are suggested as relevant predictors of distributions for amphibians and reptiles (Buckley & Jetz, 2007 ;Clusella-Trullas, Blackburn, & Chown, 2011 ).
More proximal variables relevant for our study like solar radiation, soil temperature (also temperatures under the snow cover in the winter) were unfortunately not available in our analysis. However, air temperature is a proxy for these variables; for example, higher temperatures in southern expositions come along with more solar radiation in these locations and vice versa. In our study area future minimum, median and maximum values of average temperatures increased in all scenarios, with maximum average temperatures increasing between 2 and 5 °C. For average precipitation, the future scenarios mostly showed lower maximum values and stable or slightly increasing minimum and median values. Temperature SD revealed spatial patterns of a lower-resolution grid for the future climate models. Therefore, we dropped this variable and used only average temperature, average precipitation and precipitation SD as variables (pairwise Pearson r < .8; Appendix S2 : Table A2 ). We additionally calculated MESS (multivariate environmental similarity surface) maps for the future predictors to assess the extent of model extrapolation, which can be problematic (Elith, Kearney, & Phillips, 2010 ). For the analysis of elevation patterns and aspect, we used a 100-m resolution digital elevation model (Bundesamt für Landestopografie swisstopo, 2018 ). Aspect was calculated using ArcMap 10.2.1. As aspect is a circular variable and as we emphasize to contrast northern to southern exposition, we transformed it into "northness" (cos(aspect)), which is scaled from −1 (south-exposed) to 1 (north-exposed) (Guisan, Weiss, & Weiss, 1999 ;Lassueur, Joost, & Randin, 2006 ).

F I G U R E 1
Presence locations and elevation density plots (presence vs. background locations) of the Alpine salamander ( Salamandra atra, n = 2,051) and the Common lizard ( Zootoca vivipara, n = 4,561)

| Species distribution modelling
We used an ensemble approach with three SDM algorithms to account for variation among modelling algorithms (Araújo & New, 2007 ;Watling et al., 2015 ). We used boosted regression trees (BRT) (Elith, Leathwick, & Hastie, 2008 ), generalized additive models (GAM) (Guisan, Edwards, & Hastie, 2002 ) and Maxent (Phillips, Anderson, Dudík, Schapire, & Blair, 2017 ;Phillips, Anderson, & Schapire, 2006 ). We ran each model ten times, which was enough to capture the environmental background variability (Appendix S2 : Table A3 ). In each run, we randomly used 75% of the presences for model training and the remaining 25% as test presences. To account for a possible sample selection bias, we choose the target-group background approach for the selection of pseudo-absence and (in Maxent) background points (Phillips et al., 2009 ). That is, we randomly chose 10,000 background points for model training and 2,500 as test points out of the ca. 21,000 sampled locations (which share the possible sample bias of presences) (Barbet-Massin, Jiguet, Albert, & Thuiller, 2012 ). In each replicate run, we then projected the models into the different future scenarios. Finally, for current and future predictions we calculated the mean of all replicates within and between algorithms to get our final ensemble suitability maps (from now on we refer to these ensemble maps when talking about suitability maps) (Marmion, Parviainen, Luoto, Heikkinen, & Thuiller, 2009 ).
For model evaluation, we present test AUC (area under the receiver operating characteristic curve) as a threshold-independent accuracy metric. Although it is rather a measure of how restricted a distribution is than a measure of performance, AUC is useful to compare model outputs for the same species in the same study region (Lobo, Jiménez-Valverde, & Real, 2008 ;Merow, Smith, & Silander, 2013 ).

| Threshold selection
As our aim was to compare elevation and northness of current and future locations of the study species, we first had to define when a current location became unsuitable and where the closest and therefore most likely future destination was located. To do the former, we transformed suitability maps into binary presence-absence maps.
This can be a critical step, as these transformed maps can heavily depend on the choice of threshold (Fielding & Bell, 1997 ;Merow et al., 2013 ;Nenzén & Araújo, 2011 ). We thus used four threshold rules in our analysis to study whether the results were consistent between them. The threshold dependency issue is often neglected (but see Nenzén & Araújo, 2011 ;Steen, Sofaer, Skagen, Ray, & Noon, 2017 ); therefore, we present our approach here in detail. We used two "standard" threshold rules, which are applicable to the entire study area, and a more biologically justified location-dependent approach in two implementations. As standard rules, we chose the average suitability at presence locations ("avg_suit") ( S. atra : 0.407, Z. vivipara : 0.467) and the maximizing the sum of test sensitivity and specificity ("maxSSS") threshold rules, which both perform well even for presence-only data (Liu, Berry, Dawson, & Pearson, 2005 ;Liu et al., 2013 ;Nenzén & Araújo, 2011 ;Steen et al., 2017 ). To calculate the maxSSS threshold, as in model evaluation, we randomly chose 25% of the presences and 2,500 pseudo-absences (background points) and extracted the suitability scores at these locations in 100 replicates. The R package "presenceabsence" (v1. In these approaches, one threshold value is applied to every grid cell of the study area to define whether the species is present or absent therein. These threshold rules inevitably lead to false-negative current predictions of actual presence locations with suitability scores below this threshold value, which biases the downstream analyses (cf. distribution of current suitability at presence locations in Appendix S1 : Figures A2 and A3 ). All studies using SDM to calculate range changes in species have these caveats, and this approach is generally accepted as it is the only way to calculate binary presence-absence maps. To avoid these conceptual caveats, we propose an additional approach that focuses on a location-dependent loss of relative suitability. We defined a location as unsuitable if the future suitability was lower than a certain percentage of its current suitability. Due to the probabilistic nature of predictions, this might be biologically more reasonable than calculating one universal threshold value for all presence locations. If there are erroneous presences in the data set (e.g. misidentifications, migrating individuals in unsuitable locations), this approach might give them more weight.
However, as our study species are very unique, have small home ranges and as we used a very large number of observations, erroneous presences are very unlikely to occur in our data. We considered two thresholds to define unsuitable conditions: a decrease of more than 20% ("0.8_suit") and a more conservative decrease of more than 40% ("0.6_suit") of the current suitability of each presence location. As this approach calculates a different threshold value for every presence location (and only them), it cannot be used to calculate binary presence-absence maps.

| Data analysis
Our final analysis included the following steps and was done for each threshold and climate scenario, respectively. We first compared current and future suitability at each presence location to identify which locations were predicted to become unsuitable. For these locations, we looked for the nearest predicted suitable location, which we assumed to be the most likely destination. This approach might include highly improbable situations where species have to cross inaccessible terrain or barriers (valleys or mountain tops) to reach that nearest suitable location. However, due to our large number of presences and the resulting distances, these cases are not likely to affect our main conclusions, which describe overall trends. We then extracted elevation and northness of the new location and calculated the horizontal distance to the original location. This 1st-Nearest Neighbour analysis was done using the "knnLookup" function of the R package "searchtrees" (v0.5.2, Becker, 2012 ). For locations predicted to remain suitable, neither elevation nor aspect shift was recorded.
We used a Wilcoxon signed rank-test to examine the potential shift in the distributions of current and future elevation and northness for the entire data set (paired data). For further analysis of northness patterns, we additionally divided the data set into lower and higher elevation locations to account for a possible relationship between aspect preferences and elevation (Dobrowski, 2011 ).
We ran a linear model on the current data to test for a dependency between northness and elevation. We then used an elevation of 1,700 m to separate lower and higher elevations, which represents the transition zone of the montane and alpine zone of the central Alps (which is not constant throughout its north-south profile) (Veit, 2002 ). When comparing northness between high-elevation and low-elevation locations within climate models or between current and future periods, we used a Wilcoxon rank sum-test (unpaired data). We additionally ran Anderson-Darling (AD) tests to compare the shapes of these distributions. We chose this test, as it is more sensitive towards differences at the tails of distributions, which are of special interest in comparing northness distributions compared to, for example the Kolmogorov-Smirnoff test (Engmann & Cousineau, 2011 ). To assess the impact of climate change on the geographical distributions of our study species, we ran some further analyses. We calculated Schoener ' s D ("dismo" package) (Rödder & Engler, 2011 ;Warren, Glor, & Turelli, 2008 ), a measure of niche overlap, to compare the similarity of current and future SDM predictions in geographical space. Finally, future range changes in the binary predictions ((occupied future grid cells/occupied current grid cells)-1) were calculated for the maxSSS and avg_suit threshold rules (the other threshold rules cannot be used to build binary maps)).

| Distribution modelling and future changes
All SDM performed well for both species (test AUC > 0.77, Appendix S2 : Table A3 ). The modelled response curves corresponded to biological expectations (see Section 4 ) and were consistent across models (Appendix S2 : Figures A1-A7 ). While both average precipitation and average temperature contributed similarly to the salamander models, lizard models were mainly influenced by average temperature (Appendix S2 : Figures A1-A8 ). SDM predictions for all periods and scenarios are presented in Figures A9-A14 (Appendix S2 ). MESS maps revealed that only small areas were outside training conditions for most scenarios (Appendix S2 : Figures A15-A26 ). Suitability at the presence locations decreased for both study species in almost every future scenario (Table 1 ). When applying the thresholds, the proportion of locations with unsuitable future conditions were strongly dependent on the threshold rule. The 0.6_suit threshold was the most conservative, while the avg_suit rule leads to the highest number of unsuitable locations. Differences between thresholds were smaller and less variable for the lizard models. For locations that became unsuitable under future conditions, in most cases suitable future locations were only a couple of kilometres away (Appendix S1 : Tables A1 and A2 ). For instance, the median third quartile of these distances over all future scenarios and thresholds was 1.4 km (± 1.0 km SD) for the Alpine salamander and 5.0 km (± 5.7 km SD) for the Common lizard. When using Schoener ' s D to assess the similarity between current and future predictions, we found strong heterogeneity, with some predictions showing high similarity, and others showing moderate dissimilarity ( Table 2 ). The conversion of the continuous maps into binary maps revealed different range change patterns among scenarios and study species, with partly large differences between thresholds (Table 2 ). For the Alpine salamander, in half of the scenarios the suitable area increased (by up to 194%), while it decreased in the other half (by up to 85%). In contrast, binary maps showed an overall decrease in suitable area for all scenarios (by 9% -40%) for the lizard.

| Shifts in elevation and aspect
We found significant future shifts in elevation (for the paired data) and changes in its distribution between current and future locations in almost every scenario for both study species (Appendix S1 : Tables A4 and A5 ). The effect of these shifts on the overall elevation distribution heavily depended on the climate scenario, the period and especially the threshold rule (Figures 2 and 3 , Appendix S3 : Figures A2-A13 ). While for the most conservative scenarios and thresholds, there was no or only a slight upward shift at the lower range limit, in the extreme cases lower ranges shifted up to 750 m upwards for S. atra , even up to 1,100 m for Z. vivipara (Appendix S1 : Table A3 ). We also observed a shift of northness under all the future scenarios and thresholds for both species (Appendix S1 : Tables A4 and A5 ). Overall, northness distributions differed also for all scenarios for both species, except a few scenarios for S. atra . The elevation-dependent analysis of northness revealed some more detailed patterns (Appendix S3 : Figures A2-A13 ). The linear regressions between northness and elevation were highly significant for both species, but they explained almost no variance of the current data (adj. r 2 around .03 for both species). When comparing the distributions of northness of locations at lower (< 1,700 m) and higher (> 1,700 m) elevations, they differed within almost every scenario and for every threshold (Appendix S1 : Tables A6-A7 ). At current lower elevation locations, salamanders showed no aspect preference, while south-exposed locations were preferred at higher elevations (Figure 2 , Appendix S3 : Figure A1 ). The lizards preferred southern expositions at both elevations, though north-exposed locations were relatively more frequent at lower elevations (Figure 3 , Appendix S3 : Figure A1 ). When accounting for elevation in the analysis between periods, the differences at lower elevations caused the northness differences observed for lizards (Appendix S1 : Table A9 ). Under future conditions at lower elevations, lizards were no longer associated with south-exposed locations, but northern expositions became as or even more frequent (cf. Figure 3 , Appendix S3 : Figures A8-A13 ).
For the salamander models, these differences between elevations were less pronounced and more scenario and threshold dependent.
For some scenarios, the northness distributions of future locations differed also from current locations at higher elevations, while for some scenarios no differences could be seen even at lower elevations (Appendix S1 : Table A8 ). In contrast to current conditions, for most scenarios at low elevations salamanders were more frequent at northern expositions in future periods (cf. Figure 2 , Appendix S3 : Figures A2-A7 ). Contrary to most lizard models, at higher elevations the salamanders in some scenarios also occupied northern expositions as often or even preferred them to south-exposed locations.

| Distribution modelling and future changes
Our high-resolution SDM performed well and reflected the biology of our study species. Temperature was the strongest predictor of the distribution of the Common lizard, which is a diurnal reptile requiring sun for thermoregulation. On the other hand, precipitation was more or equally important for the Alpine salamander, which is a mostly nocturnal amphibian strongly dependent on humid conditions (Werner, Lötters, Schmidt, Engler, & Rödder, 2013 ). Precipitation SD did not have a strong impact on the distributions of our study species. This might be an effect of the independency from small water bodies or their ability to withdraw into crack systems and crevices to shield negative effects of, for example drought, a possible negative manifestation of a high precipitation variability in humid locations. Under future conditions, the suitability generally decreased for both species, indicating a possible negative impact of climate change (Table 1 ). Currently, global emissions are tracking the trajectory of the rcp85 greenhouse gas concentration scenario (Sanford, Frumhoff, Luers, & Gulledge, 2014 ). In these scenarios, the lower elevations were predicted to become unsuitable at the end of the century for both species at the northern range edges, and the centres of their distributions were predicted to shift more or less southwards. However, for the majority of the current locations suitable alternatives were only a few kilometres away and could be colonizable by the study species, justifying the use of these future locations in our analysis. When assessing the binary future range changes in our models, the results were rather different between our study species (Table 2 ). While for the Common lizard, a minor range loss was consistently predicted for all scenarios, the predicted Alpine salamander ranges varied drastically between large gains (+194%) and losses (−85%), suggesting a strong uncertainty on the potential effects of climate change on this species (but see the next paragraph). This variation is unexpected given that climate variables were highly correlated between future scenarios (Pearson r > .98 for between-scenario correlations pairs of average temperature and precipitation) (but see Braunisch et al., 2013 ). The differences between the study species are explainable by the different species-specific variable contributions to the models, which represents a different model complexity. Models for the Alpine salamander were mainly based on two more or less equally contributing variables, resulting in more complex predictions, while mean temperature was a single strong predictor for the Common lizard. The diverse predictions for the Alpine salamander highlight the strong impact of the complex microclimatic conditions in the mountains on distribution models and the need of high-resolution data for climate change risk assessments. For instance, when information on these fine-scale conditions is missing, future SDM projections using coarse resolution data (~20 km) suggested a consistent range loss for the Alpine salamander (Beierkuhnlein, Jentsch, Reineking, Schlumprecht, & Ellwanger, 2014 ). We did not explicitly build a coarse resolution SDM with our data for comparison, as a spatial resolution of 1 km or larger just is not plausible in an environment, where you can have over 1,000 m elevation differences in such a grid size.

| Thresholding caveats
The conversion of continuous maps into binary ones using threshold values is extremely challenging. When comparing the future range change between binary maps, the differences between the climate scenarios and threshold rules varied strongly and inconsistently compared to the similarity of the original continuous model outputs (measured by Schoener ' s D ) ( We therefore share the view of Merow et al. ( 2013 ) to avoid thresholding whenever possible. It masks the probabilistic nature of predictions and the various uncertainties already associated with future predictions (Braunisch et al., 2013 ;Peterson, Cobos, & Jiménez-García, 2018 ;Steen et al., 2017 ;Wenger et al., 2013 ).
When thresholding is essential for a certain research question, it can be useful applying multiple threshold rules, which should be carefully chosen based on the existing literature (Jiménez-Valverde & Lobo, 2007 ;Liu et al., 2005Liu et al., , 2016Nenzén & Araújo, 2011 ;Steen et al., 2017 ). Furthermore, we strongly advise to compare the similarity of current and projected binary suitability maps with the F I G U R E 2 Violin plots of elevation and northness (for lower and higher elevations) of S. atra locations for current ("2000") and future periods and thresholds (a-d), exemplary for the rcp85_EC_DMI scenario. a: 0.6_suit, b: maxSSS, c: 0.8_suit, d: avg_suit. Maximum width of violin plots is proportional to the sample size. Elevation is shown in m a.s.l.. Northness is scaled from 1 (north-exposed) to −1 (south-exposed). Horizontal lines represent the first, second and third quartiles similarity of the corresponding continuous maps (as in Table 2 ) to assess the ecological plausibility of the binary range differences.

| Shifts in elevation and aspect
Fine-scale projections of species response to climate change showed both upward and northward shifts in alpine landscapes.
When comparing current and future locations, these two strategies were not mutually exclusive, as there was both a northward shift in aspect and an upward shift in elevation. Regression analysis indicated a linear relationship between northness and elevation for both species (likely due to the high number of presences), but it did not explain any of the variance in the data. However, comparing distributions of low and high elevations showed that the association of certain aspects indeed varies with elevation. At current conditions, both the Alpine salamander and the Common lizard preferred warmer microclimates on southern expositions at high elevations. At lower elevations, the association with south-exposed locations was weaker. This pattern suggests that the occupation of cool, north-facing microclimates is a strategy for dealing with warm conditions, so we expect an even more pronounced pattern under future warming conditions. In fact, in most future scenarios, a shift in the distribution of both species towards cooler north-exposed locations was predicted, especially at low elevations. This clearly indicates that many southern slopes at lower elevations may become too warm and/or too dry in the future. At higher elevations, however, this pattern is weaker, with some predicted northward F I G U R E 3 Violin plots of elevation and northness (for lower and higher elevations) of Z. vivipara locations for current ("2000") and future periods and thresholds (a-d), exemplary for the rcp85_EC_DMI scenario. a: 0.6_suit, b: maxSSS, c: 0.8_suit, d: avg_suit. Maximum width of violin plots is proportional to the sample size. Elevation is shown in m a.s.l.. Northness is scaled from 1 (north-exposed) to −1 (south-exposed). Horizontal lines represent the first, second and third quartiles shifts of the salamander and the lizard even showing no response to future warming. These results show that climate change response of biota in complex mountain topography is also complex and that high-resolution data are needed to identify terrain-based microrefugia that can mitigate the effect of climate change (Bennie, Wilson, Maclean, & Suggitt, 2014 ;Dobrowski, 2011 ;Meineri & Hylander, 2017 ;Potter et al., 2013 ). So far, it has been shown that microclimates in mountains buffer regional warming and provide nearby suitable habitats, especially for plants (and some insects) (Kulonen et al., 2017 ;Meineri & Hylander, 2017 ;Opedal et al., 2015 ;Scherrer & Körner, 2011 ;Suggitt et al., 2018 ). Our new study highlights the need to consider microclimate variation in SDM to more realistically assess the effect of climate change on ectothermic vertebrates which are usually less mobile and highly dependent on the microclimate and the existence of nearby refugia compared to endothermic species with often larger home ranges (Bennie et al., 2014 ;Ficetola et al., 2018 ;Potter et al., 2013 ).

| CONCLUSIONS
The limited potential of cold-adapted high mountain species, of which many are even mountain endemics, to respond to climate change is of major concern to conservationists. It is essentially based on the long-standing view that elevational upward shift is their main response to warming (Marris, 2007 ;Moritz et al., 2008 ;Parmesan, 2006 ;Wilson et al., 2005 ). Using high-resolution SDM data and an exceptional large numbers of precise locations, we show for the first time that the responses to climate change in two alpine vertebrate species may also be much more complex than previously thought (Kulonen et al., 2017 ;Opedal et al., 2015 ;Scherrer & Körner, 2011 ;Suggitt et al., 2018 ). Both upward and northward shifts may occur, with aspect preferences being highly dependent on elevation and vice versa, as a shift from a warmer southern to a cooler northern aspect is not decoupled from elevational shifts.
For most current locations of our two study species, there exist suitable and nearby alternative locations. Therefore, the negative effects of climate change might be mitigated, at least partially, by the complex topography of the mountains, so an upward movement may not necessarily be the only option (Suggitt et al., 2018 ). Our study also implies that the use of high-resolution data is critical to identify suitable microrefugia and thereby improving impact assessment of climate change. This will make SDM projections to future climate change scenarios more realistic and may be of particular importance to ectothermic species, such as all invertebrates, whose home ranges are smaller than the extent of microclimatic variation.

ACK N OWLED G E M E NT S
We thank Jonas Lembrechts and an unknown reviewer for valuable comments on the manuscript, which helped to improve the paper.
We also thank all the herpetologists and volunteers who submit-