Potential range shifts and climatic refugia of rupicolous reptiles in a biodiversity hotspot of South Africa

Summary

Climate change is causing the geographical ranges of some species to track suitable conditions. Habitat specialists, range-restricted species and species with limited dispersal abilities may be unable to track changing conditions, increasing their extinction risk. In response to changing conditions and species movement patterns, there is a need to account for the effects of climate change when designing protected areas and identifying potential climate refugia. We used ecological niche models projected into future climates to identify potential impacts of climate change on the distribution of 11 rupicolous reptile species in the Soutpansberg Mountains, South Africa. Lygodactylus incognitus, Lygodactylus soutpansbergensis, Platysaurus relictus and Vhembelacerta rupicola were identified as being vulnerable to climate change due to substantial reductions in suitable habitat and low spatial overlap between current and future niche envelopes. We identified areas of high conservation importance for the persistence of these species under present-day and projected future conditions. The western Soutpansberg was identified as an area of high conservation priority as it is a potential refuge under future projections. Projecting distributions of vulnerable species into future climate predictions can guide future research and identify potential refugia that will best conserve species with restricted ranges in a world with climate change.

Introduction

Climate change is one of the main threats to biodiversity, with changing temperatures and alterations in rainfall patterns predicted on a global scale (Trenberth 2011, Brown & Caldeira 2017). In response to changing climates, species’ geographical ranges are shifting to track suitable conditions, with movements usually aligning along latitudinal and elevational gradients (Chen et al. 2011, Botts et al. 2015, Tiberti et al. 2021). However, physical barriers and limited dispersal ability hinder the expansion of ranges into suitable areas for some species, resulting in an overall contraction in distribution, making them vulnerable to the negative effects of climate change (Foden et al. 2008). These increase the risk of extinction for these species (Thomas et al. 2004, Foden et al. 2008).

Mountains have been important climatic refugia for many species in the past due to their heterogeneous topography and elevational range (Bennett et al. 1991, Picker & Samways 1995). The presence of refugia has also supported the maintenance of biodiversity and speciation through evolutionary history, particularly during times of climatic extremes, because they act as safe havens for many species under these conditions (Keppel et al. 2012) and promote the allopatry that is conducive to speciation. However, areas that have functioned as refugia in the past may not necessarily do so for future climate change (Keppel et al. 2012) since change is complex, multidimensional and may not be simply repetitive. Areas where current and future distributions will likely overlap have the greatest potential to act as future climatic refugia (Hugall et al. 2002, Keppel et al. 2012). The identification of these areas of overlap is crucial for long-term conservation planning.

The likely effects of climate change on species distributions are being seen as increasingly important when designing protected areas (PAs) (Heller & Zavaleta 2009). Incorporating estimates of how species distributions are predicted to move, habitat connectivity between these areas and potential climatic refugia into PA designation will improve the effectiveness of PAs under a changing climate. Ecological niche models (ENMs) are increasingly used as a tool to estimate how climate change will likely affect species distributions in the future and to identify priority areas for conservation planning (Nori et al. 2018, Zhu et al. 2021). ENMs are especially useful when modelling restricted and understudied species with limited species-specific data (Pearson & Dawson 2003, Fordham et al. 2018). Thus, on condition that modelling is based on realistic predictions of future climate, this approach can provide information for long-term conservation strategies and planning.

Identifying priority conservation areas is particularly important for range-restricted species that are susceptible to climate change and have limited dispersal abilities. In this context, rupicolous species are especially at risk since their distribution is dependent on the presence of suitable rock substrates, often resulting in patchy occurrence (Penman et al. 2010, Croak et al. 2012). Since rocky habitat is often fragmentary and does not change in response to climatic shifts (Penman et al. 2010), climate tracking in rupicolous species may be dependent on serial jump dispersal events rather than diffusion or secular migration (sensu Mason 1954). This results in rupicolous species being less able to climate track, especially when climate change is rapid, and suitable rock substrates may not exist in the new niche envelope space. Thus, it is important to identify potential climatic refugia and suitable areas of connectivity between current and future niche envelopes for these species.

The Soutpansberg Mountains in Limpopo Province, South Africa, have extremely high reptile species richness, being the most species-rich area in South Africa (Tolley et al. 2019), with more than 100 species occurring on the mountain range (Alexander 2009). Many of these species have restricted ranges, and 16 are strictly rupicolous, of which 5 are endemic to the mountains (Petford et al. 2019). We evaluated the likely impacts of climate change on these rupicolous species using ENMs. We focused on rupicolous reptiles as they are likely to be highly sensitive to changing climatic conditions due to their limited dispersal ability, particularly for endemics, as range-restricted species likely have narrower climatic niches (Slatyer et al. 2013). We hypothesized that the endemic species in our study are the most likely to be negatively impacted by change, with substantial predicted decreases in future range size compared to those with wider distributions (sensu Botts et al. 2015).

Materials and methods

Study area

The Soutpansberg Mountains are located in the north-eastern corner of Limpopo Province, South Africa (Fig. 1). Extending over 210 km from Pafuri in the east to Vivo in the west and ranging from 250 to 1748 m above sea level, the mountains have high habitat, topographical and climatic heterogeneity (Hahn 2006, Alexander 2009).

Fig. 1. The Soutpansberg Mountains study site in relation to South Africa.

The mountain range has three predominant homoclinal ridges on an east–west axis, each pushed up in the south, forming steep southern cliffs with more gentle slopes to the north (Hahn 2006, 2011). The Sand River Valley is the only corridor bisecting the range in the western Soutpansberg. The northern side of the mountain is more arid than the south. Waterpoort (northern slope) receives 367 mm of rainfall per annum, while Entabeni (southern slope) receives 1874 mm of rainfall per annum (Hahn 2006). The Soutpansberg similarly becomes drier and hotter towards the east, with tropical elements entering the eastern extremities. The axis of the mountains also results in moisture-rich air from the Indian Ocean transporting mist to the southern and central slopes at high elevations; this effectively doubles annual rainfall measurements in these areas (Hahn 2006).

Species occurrence records

Occurrence records were gathered for 11 rupicolous reptiles occurring on the Soutpansberg Mountains: Afroedura pienaari (endemic), Chondrodactylus turneri, Cordylus vittifer, Lygodactylus incognitus (endemic), Lygodactylus soutpansbergensis (endemic), Platysaurus intermedius, Platysaurus relictus (endemic), Smaug depressus, Trachylepis margaritifer, Trachylepis punctatissima and Vhembelacerta rupicola (endemic). See Supplementary Appendix S1 (available online) for details. Records were gathered during fieldwork from 2014 to 2018 and from the Global Biodiversity Information Facility (GBIF 2018, 2019). The oldest locality records postdate any significant transformation of the land. The dataset was inspected for spatial errors and cleaned. Combined datasets resulted in 3875 locality records (7.8% GBIF data and 92.2% own data). To reduce spatial autocorrelation, spThin (Aiello-Lammens et al. 2015) was run in R 3.5.1 (R Core Team 2018) to spatially filter occurrence points using a 1-km radius. This reduced the number of records to 366 (n = 12–79 per species; Table S1).

Environmental variables

Nineteen bioclimatic variables for present and future conditions were downloaded from www.worldclim.org at a 30 arc-second (≍900 m at the equator) resolution. Elevation (GTOPO30 (30 arc-second resolution); Earth Resources Observation and Science Center et al. 1997), geology (Council for Geoscience 2018 (scale = 1:1 000 000)) and vegetation type (bgis.sanbi.org (scale = 1:1 000 000)) were also downloaded. To reduce collinearity effects, a Pearson’s correlation was performed on all environmental variables in R 3.5.1. Variable pairs with r ≥ 0.75 were inspected and the variable considered to be the least important for the distribution of rupicolous reptiles was removed. This was assessed by considering the biology of our study species, the biological system of our study area and primary limiting factors for ectotherm distribution, specifically prioritizing variables relating to extremes and variability (Pintor et al. 2016, Bradie & Leung 2017). The remaining variables used to perform ENMs were: mean diurnal temperature range (BIO2); isothermality (BIO3); temperature seasonality (BIO4); mean temperature of the coldest 3 months (BIO11); annual precipitation (BIO12); precipitation seasonality (BIO15); and geology and vegetation type (categorical).

Due to uncertainties surrounding future climatic conditions, three global circulation models (GCMs) were used for 2070 projections: Hadley Global Environment Model 2 – Carbon Cycle (HadGEM2-CC); Model for Interdisciplinary Research on Climate 5 (MIRCP5) and Goddard Institute for Space Studies E2-R (GISS-E2-R). These GCMs were selected using McSweeney et al. (2015) as a guide for which models work well over Africa. Two Representative Concentration Pathways (RCPs) of 4.5 and 8.5 were selected. These are considered pathways with moderate and high predicted CO₂ concentrations (IPCC 2013). Models with an RCP of 4.5 predict that carbon emissions will peak in 2040 and mean global temperature will increase by 1.8°C by 2100 (IPCC 2013). This RCP was selected over 2.6 as it is a more moderate estimate, with 2.6 RCP predicting decreasing carbon emissions throughout the century (IPCC 2013). An RCP of 8.5 is the worst-case scenario that predicts carbon emissions will continue to rise until 2100 with a temperature increase of 3.7°C by 2100 (IPCC 2013).

Ecological niche models

Species distributions were modelled using a maximum entropy (Maxent) approach, with presence-only species occurrence records and randomly generated background points (Phillips et al. 2006). Maxent was chosen due to its credibility in modelling species environmental niches with presence-only data and small sample sizes (Elith et al. 2006, Pearson et al. 2007). Kaky et al. (2020) demonstrated that Maxent performs equally to an ensemble approach, but with less computational time, and a recent study found no benefit to using an ensemble approach over an appropriately tuned model (Hao et al. 2020). Thus, the ENMeval package (Muscarella et al. 2014) in R 3.5.1 was used to construct Maxent models with different parameter settings and perform model evaluations.

Models were built with different combinations of linear (L), quadratic (Q), hinge (H), product (P) and threshold (T) feature classes (LQHPT; LQHP; LQH; L; LQ; H) and levels of regularization (0.5–4.5 with increments of 0.5), resulting in 48 different combinations of model runs for each species. To reduce spatial sampling bias, 10 000 background points were selected from a bias layer (Phillips et al. 2009) containing 8111 records of all reptile data points from fieldwork and all rupicolous records accessed from GBIF. Data were partitioned into testing and training bins using the ‘block’ method, as recommended when projecting data across time by Muscarella et al. (2014). The block method partitions occurrence and background data into four bins based on longitude and latitude, three of which are used for training while the fourth is retained for testing; this is repeated four times, with a different bin used for testing each time (Muscarella et al. 2014).

Optimal settings for each species were selected using four different criteria following Muscarella et al. (2014). Firstly, the model with the lowest Akaike information criterion metric corrected for small samples sizes (AIC_c) was inspected, as these are considered to be best in terms of goodness of fit and complexity (Burnham & Anderson 2004, Warren & Seifert 2011). Following this, the threshold-independent metric (area under the curve; AUC), the difference between the test and training AUC (AUC_diff) and the minimum training presence omission rate (OR_mtp) were assessed to ensure models could discriminate and prevent overfitting (Anderson & Gonzalez 2011). Response curves were checked for ecological realism (Guevara et al. 2017), and the most appropriate models were projected onto 2070 climate scenarios.

For comparison between current and future distributions, simple mean ensemble approaches (Palm & Zellner 1992) were used to create combined estimates from the three GCMs for each species at 4.5 RCP and 8.5 RCP. The resulting maps were modified into binary maps (1 = predicted suitable habitat; 0 = predicted unsuitable habitat). The binary maps were created using individual threshold values for current distributions using the 10-percentile training presence value, a conservative method allowing for potential errors in the GBIF datasets (Liu et al. 2005). Predicted suitable habitat was calculated using the GRASS r.report function in QGIS 3.4.2, and the average elevation of the suitable areas were derived using the SAGA zonal raster statistic. Schoener’s D statistic for niche overlap was calculated between current and future predictions using the nicheOverlap function of the dismo package (Hijmans et al. 2017) in R 3.5.1. A Kruskal–Wallis H test was used to identify significant differences in change of suitable habitat between endemic and widespread species using IBM SPSS Statistics for Windows Version 26 (IBM Corp. 2019). A non-parametric test was used as the variables did not meet the assumptions of normality or heteroscedasticity.

Priority areas for conservation

The program Zonation 4.0 (Moilanen et al. 2009) was used to highlight areas of conservation priority for the rupicolous species of the Soutpansberg whilst taking into account climate change projections and connectivity between present and potential future distributions. Zonation selects areas of high conservation priority by using a hierarchal prioritization based on biodiversity features such as species distributions. The Core-area Zonation model was used, which prioritizes areas of high species richness while accounting for species with high priority ratings; therefore, areas with low species richness can be highlighted as important if they contain species with higher weightings (Moilanen & Kujala 2008).

Species weightings are user-defined (Lehtomäki et al. 2016). Species not endemic to the study area were weighted with a value of 1 and endemic species were weighted with a value of 2. Vhembelacerta rupicola was assigned an additive weighting of +1 as it is monotypic. Due to the uncertainty associated with projections into the future, weights given to species for future projections were half of the values described above (Table S2) (Kujala et al. 2013).

Transformed areas with a 500-m buffer were assigned negative weights to ensure that these areas were not highlighted as high priority by the Zonation program. Transformed areas include human settlements and agricultural and silvicultural areas. Unsuitable areas were selected using a land-use map – DEA/CARDNO SCPF002: Implementation of Land Use Maps for South Africa (GEOTERRAIMAGE – 2014).

Connectivity between current and future projections was accounted for by implementing the ecological interactions connectivity function (Moilanen & Kujala 2008). This function ensures that conservation prioritization is assigned to areas where current and future distributions of species overlap and to areas providing connectivity between the two with a dispersal kernel. Reptiles are not considered vagal and are thus not expected to easily track climate (Kujala et al. 2013). Estimating species-specific dispersal metrics requires information on dispersal ability (e.g., Salas et al. 2017), which we do not have for our study species; therefore, the dispersal kernel was set to 0.2 km/year, a moderate estimate used in Kujala et al. (2013), who set the parameter based on the available literature information.

Zonation output creates a map with cells ranging in value from 0 to 1, with 1 being cells considered the highest conservation priority. For this study, we allocated cells with a score of >0.85 as high conservation priority. The post-processing function Landscape Comparison (LSC) was used to identify differences between high-priority areas identified in current and future climate maps. The high-priority areas were compared to the South African Protected Area Database (SAPAD 2020) to measure the percentage of current PAs included in the Zonation-derived high-priority areas using QGIS 3.4.2.

Results

Maxent ENMs

The results from the ENM evaluation revealed that all models performed from fair to excellent with AUC values ranging between 0.70 and 0.92 (Table 1). Optimum models in each species had different parameter settings (Table S3). The training omission rate was seemingly high for A. pienaari, L. incognitus and C. vittifer, suggesting that these models may be overfitting. In general, climatic variables had a larger permutation importance than geology or vegetation for all species (Table S4). The distributions of L. incognitus, L. soutpansbergensis and P. intermedius were largely influenced by BIO11, while the distribution of T. punctatissima was mostly influenced by BIO15. All other species had large permutation importance values across two or more variables. No species were strongly influenced by BIO3 or BIO4.

Table 1. Sample sizes and evaluation metrics of the optimum maximum entropy (Maxent) models constructed using the package ENMeval for 11 rupicolous species present in the Soutpansberg Mountains, South Africa. Metrics shown are Akaike information criterion corrected for small sample sizes (AIC_c), threshold-independent metric (AUC_test), difference between test and training AUC (AUC_diff) and the minimum training presence omission rate (OR_mtp).

AUC = area under the curve.

Future predictions

Of the 11 distributions projected into future climate, 7 were predicted to undergo decreases and 3 were predicted to experience increases in suitable habitat by 2070 under both 4.5 and 8.5 RCP (Fig. 2 & Table 2). One species, S. depressus, showed a marginal predicted increase for 4.5 RCP (11%) and a marginal predicted decrease for 8.5 RCP by 2070 (1%). Of the seven species with predicted decreases, L incognitus, L. soutpansbergensis, P. relictus and V. rupicola are considered particularly vulnerable, with reductions in suitable habitat of more than 85% for 8.5 RCP, with increases in average elevation and relatively low niche overlaps (Table 2). All species predicted to show large decreases in suitable habitat shifted westwards (Fig. 2).

Fig. 2. Binary suitable/unsuitable maps showing the predicted distributions of Afroedura pienaari, Chondrodactylus turneri, Cordylus vittifer, Lygodactylus incognitus, Lygodactylus soutpansbergensis, Platysaurus intermedius, Platysaurus relictus, Smaug depressus, Trachylepis margaritifer, Trachylepis punctatissima and Vhembelacerta rupicola in current and ensemble future 2070 scenarios. Black = suitable; grey = unsuitable. RCP = Representative Concentration Pathway.

Table 2. Predicted changes in suitable habitat for 11 rupicolous species from the current climate to the 2070 climate for two different carbon emission scenarios based on maximum entropy (Maxent) model outputs. Predicted habitat loss (%), predicted habitat gain (%) and predicted habitat to remain suitable (%) to the nearest 1% were calculated by comparing areas predicted as suitable and unsuitable for the current climate with those predicted for future climates. Net predicted suitable habitat change (%) indicates the overall loss or gain of habitat. Overlap of niche envelopes was calculated by comparing the current distribution with each future prediction. Predicted area of suitable habitat is rounded to the nearest 10 km².

RCP = Representative Concentration Pathway.

The Kruskal–Wallis H test identified a statistically significant difference in net predicted suitable habitat change between endemic and widespread species in both 4.5 RCP (H(1) = 6.82, p = 0.006) and 8.5 RCP (H(1) = 5.77, p = 0.010), with endemic species typically having larger reductions in suitable habitat.

Priority areas for conservation

Zonation outputs were colour coded to identify areas with a conservation priority score above 85%, highlighting the high-priority areas in the Soutpansberg. The current prediction and all future connectivity scenarios identified the western part of the Soutpansberg Mountains as most important for the conservation of the rupicolous reptiles, with other important areas on ridges in the east (Fig. 3). Notably, there are few differences between current and future connectivity scenarios, with only slight deviations occurring in the far eastern Soutpansberg (Fig. 3). Few of the currently existing PAs in the Soutpansberg overlap with the Zonation priority areas (Fig. 4), with a 19.8% overlap.

Fig. 3. Zonation outputs showing the high-priority areas for each scenario and Landscape Comparison (LSC) maps to identify differences between high-priority areas in current and future scenarios. The current scenario includes degraded land and weighted current distributions only, whilst all future scenarios consider current distribution, degraded land and weighted species with connectivity to future predictions. Arrows are used to highlight small areas in the LSC maps. RCP = Representative Concentration Pathway.

Fig. 4. Zonation outputs highlighting the high-priority areas for conservation using current scenarios with existing protected areas (South Africa Protected Areas Database (SAPAD_OR_2020_Q1)).

Discussion

Of the 11 species modelled, suitable habitat was predicted to increase for 3 and to decrease for 7 across all climate scenarios. One species showed a marginal increase for 4.5 RCP and a marginal decrease for 8.5 RCP. Of the seven species with predicted decreases, four were highlighted as being particularly vulnerable due to decreases in suitable habitat of more than 85% in at least one RCP; these were all species endemic to the Soutpansberg Mountains: L. incognitus, L. soutpansbergensis, P. relictus and V. rupicola. The suitable habitats of the vulnerable species were predicted to track upslope and shifted westward, indicating that species may also move along a longitudinal gradient. The Zonation output recommended that the western Soutpansberg is an area of high conservation priority, both at present and when accounting for potential impacts of climate change on the species distributions.

Most of the species modelled do not appear to be especially vulnerable to the anticipated climate changes, with suitable habitat either predicted to increase (C. turneri, P. intermedius, S. depressus (4.5 RCP), T. punctatissima) or predicted to show minor decreases (T. margaritifer, S. depressus (8.5 RCP)). These species have wide-ranging distributions spread across much of southern Africa. Since wide-ranging species tend to have broader niches, they have improved capacity to track and adapt to changing climatic conditions. However, this is not the case for C. vittifer. Although the range of this species extends much further south, with the Soutpansberg making up the northern extremes of the range, it is predicted to experience a decline in suitable habitat in the Soutpansberg of just over 40% by 2070 under 8.5 RCP. The effect that climate change will have on the persistence of C. vittifer throughout its range goes beyond the scope of this study, yet it highlights the fact that this species will likely experience range contractions in parts of its distribution in the future.

We hypothesized that Soutpansberg endemics with restricted ranges would be the most vulnerable to climate change, and the four species predicted to be most at risk are indeed all endemic to the mountains. Lygodactylus incognitus, L. soutpansbergensis, P. relictus and V. rupicola are predicted to undergo substantial range reductions of over 85% by 2070 under 8.5 RCP, and range reductions of over 75% for L. incognitus and L. soutpansbergensis and of over 65% for P. relictus are predicted under 4.5 RCP. Conversely, in comparison, the Soutpansberg endemic A. pienaari is predicted to show a smaller decrease in suitable habitat under both RCPs (–27% (4.5 RCP); –28% (8.5 RCP)). The smaller reduction in A. pienaari range compared to the other endemics may be because it has a wider distribution across the warmer northern slopes of the mountains. Therefore, its tolerance of warm areas may result in A. pienaari not being as vulnerable to the warming effects of climate change. The Soutpansberg endemics at greater risk are those with the smallest distributions and thus being most likely to occur over a smaller climatic range.

Risk of extinction due to climate change is strongly linked to the ability to track favourable conditions (Pearson 2006). Not only do L. incognitus, L. soutpansbergensis, P. relictus and V. rupicola have low capacity to track climate change due to their restricted ranges, their dispersal ability is further limited by their rupicolous habits. Our models show that the overlap of niche envelopes between current and future scenarios for these species is low, particularly for the highest RCP (Table 2). This low overlap indicates that these species may not be able to track to suitable areas due to large expanses of unsuitable habitat acting as barriers, particularly if climatic conditions change rapidly (Davis & Shaw 2001, Brooker et al. 2007). Rapidly changing conditions and large dispersal distances are particularly problematic when the current suitable habitat is already fragmented due to the organism’s habitat requirements, as is predicted for L. incognitus (Petford & Alexander 2020), and due to anthropogenic habitat modification (Honnay et al. 2002, Opdam & Wascher 2004).

Complete loss of suitable habitat, as defined by our modelling, does not guarantee extinction, as correlative ENMs do not take into account adaptive capabilities or phenotypic plasticity, which may facilitate adaptation under changing conditions (Seebacher et al. 2015). However, due to a lack of data, there is much uncertainty about the impact of phenotypic plasticity in native ranges of species and whether it will alleviate the risk of local extinctions (Urban et al. 2014, Buckley et al. 2015, Wiens et al. 2019). The distributions of the species examined are also likely to be impacted by biotic factors that can also affect range edges (Wiens et al. 2009). Nonetheless, our predictions of range changes provide the most informed assessment on the impact of climate change for these species, highlighting areas that are likely to act as refugia in the face of climate change. Thus, the four species predicted to be vulnerable to climate change require conservation measures to protect areas where their distributions are likely to occur in the future in order to minimize extinction risk.

Our results indicated that the western Soutpansberg is likely to act as a refuge and is thus an area of high conservation priority presently and under all tested future carbon emission scenarios. All of the modelled endemic species and the majority of the widespread species occur in this portion of the western Soutpansberg. Furthermore, all modelled species currently occurring there are predicted to persist there under all climate change scenarios tested, and, importantly, this includes the four species identified as most vulnerable. It is therefore likely that this region will act as a refuge not only for the rupicolous species studied here, but also for many species that occur in the area.

Currently, only 19.8% of existing PAs in the Soutpansberg occur in areas considered of high conservation priority. Current land use in the western Soutpansberg is predominantly game farming with a few silviculture and agriculture areas. The low levels of land transformation mean that there is a high potential for large portions of this region to be declared as formal PAs. In addition, the Endangered Wildlife Trust (EWT), a non-governmental organization for animal conservation, is now active in the area and intends to establish a large PA (EWT 2020). This proposed PA is in the centre of the area that Zonation predicted as a priority.

Despite the pitfalls associated with modelling future distributions, models can provide predictions to strengthen conservation plans and to promote the need to consider future distributions of restricted species in the light of a changing climate (Carroll et al. 2010). Additionally, they can highlight future research needs and identify which species are especially vulnerable. Our study revealed that four Soutpansberg endemics were vulnerable to the impacts of climate change, but that there are areas that could act as refugia to aid their persistence. Future studies in the region should focus on identifying shifts in species distributions to observe whether they follow those predicted here. More in-depth species-specific studies on physiological capabilities and tolerances would also be beneficial, providing a basis to observe climatic resilience and to construct mechanistic niche models.

In this study, we identified potential climatic refugia for rupicolous reptile species in the Soutpansberg Mountains and demonstrated their importance for conservation planning in response to impending climate change, as well as showing that endemic species with restricted ranges will likely be more vulnerable. We contend that our method applies to conservation planning for other mountainous areas where range-restricted species may be threatened by change. ENMs provide a robust exploratory tool for assessing likely responses to change, guiding the placement of PAs, identifying areas that are important for habitat connectivity between present-day and future niche envelopes and identifying potential climatic refugia. By taking future distributions of restricted species into account when designing PA networks, these reserves will be more effective at ensuring the persistence of species that are vulnerable to climate change, a necessity in a changing world where anthropogenic activities have placed a significant number of species at risk.