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 Reference Trenberth2011, Brown & Caldeira Reference Brown and Caldeira2017). 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. Reference Chen, Hill, Ohlemüller, Roy and Thomas2011, Botts et al. Reference Botts, Erasmus and Alexander2015, Tiberti et al. Reference Tiberti, Mangiacotti and Bennati2021). 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. Reference Foden, Mace, Vié, Angulo, Butchart, DeVantier, Vié, Hilton-Taylor and Stuart2008). These increase the risk of extinction for these species (Thomas et al. Reference Thomas, Cameron, Green, Bakkenes, Beaumont and Collingham2004, Foden et al. Reference Foden, Mace, Vié, Angulo, Butchart, DeVantier, Vié, Hilton-Taylor and Stuart2008).
Mountains have been important climatic refugia for many species in the past due to their heterogeneous topography and elevational range (Bennett et al. Reference Bennett, Tzedakis and Willis1991, Picker & Samways Reference Picker and Samways1995). 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. Reference Keppel, Van Niel, Wardell-Johnson, Yates, Byrne and Mucina2012) 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. Reference Keppel, Van Niel, Wardell-Johnson, Yates, Byrne and Mucina2012) 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. Reference Hugall, Moritz, Moussalli and Stanisic2002, Keppel et al. Reference Keppel, Van Niel, Wardell-Johnson, Yates, Byrne and Mucina2012). 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 Reference Heller and Zavaleta2009). 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. Reference Nori, Leynaud, Volante, Abdala, Scrocchi and Rodriguez-Soto2018, Zhu et al. Reference Zhu, Papeş, Giam, Cho and Armsworth2021). ENMs are especially useful when modelling restricted and understudied species with limited species-specific data (Pearson & Dawson Reference Pearson and Dawson2003, Fordham et al. Reference Fordham, Bertelsmeier, Brook, Early, Neto and Brown2018). 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. Reference Penman, Pike, Webb and Shine2010, Croak et al. Reference Croak, Pike, Webb and Shine2012). Since rocky habitat is often fragmentary and does not change in response to climatic shifts (Penman et al. Reference Penman, Pike, Webb and Shine2010), climate tracking in rupicolous species may be dependent on serial jump dispersal events rather than diffusion or secular migration (sensu Mason Reference Mason1954). 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. Reference Tolley, Weeber, Maritz, Verburgt, Bates and Conradie2019), with more than 100 species occurring on the mountain range (Alexander Reference Alexander2009). Many of these species have restricted ranges, and 16 are strictly rupicolous, of which 5 are endemic to the mountains (Petford et al. Reference Petford, van Huyssteen and Alexander2019). 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. Reference Slatyer, Hirst and Sexton2013). 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. Reference Botts, Erasmus and Alexander2015).
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 Reference Hahn2006, Alexander Reference Alexander2009).
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 Reference Hahn2006, Reference Hahn2011). 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 Reference Hahn2006). 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 Reference Hahn2006).
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. Reference Aiello-Lammens, Boria, Radosavljevic, Vilela and Anderson2015) 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. Reference Pintor, Schwarzkopf and Krockenberger2016, Bradie & Leung Reference Bradie and Leung2017). 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. (Reference McSweeney, Jones, Lee and Rowell2015) 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 CO2 concentrations (IPCC Reference Stocker, Qin, Plattner, Tignor, Allen and Boschung2013). 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 Reference Stocker, Qin, Plattner, Tignor, Allen and Boschung2013). 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 Reference Stocker, Qin, Plattner, Tignor, Allen and Boschung2013). 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 Reference Kujala, Moilanen, Araújo and Cabeza2013).
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. Reference Phillips, Anderson and Schapire2006). Maxent was chosen due to its credibility in modelling species environmental niches with presence-only data and small sample sizes (Elith et al. Reference Elith, Graham, Anderson, Dudík, Ferrier and Guisan2006, Pearson et al. Reference Pearson, Raxworthy, Nakamura and Peterson2007). Kaky et al. (Reference Kaky, Nolan, Alatawi and Gilbert2020) 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. Reference Hao, Elith, Lahoz-Monfort and Guillera-Arroita2020). Thus, the ENMeval package (Muscarella et al. Reference Muscarella, Galante, Soley-Guardia, Boria, Kass, Uriarte and Anderson2014) 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. Reference Phillips, Dudík, Elith, Graham, Lehmann, Leathwick and Ferrier2009) 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. (Reference Muscarella, Galante, Soley-Guardia, Boria, Kass, Uriarte and Anderson2014). 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. Reference Muscarella, Galante, Soley-Guardia, Boria, Kass, Uriarte and Anderson2014).
Optimal settings for each species were selected using four different criteria following Muscarella et al. (Reference Muscarella, Galante, Soley-Guardia, Boria, Kass, Uriarte and Anderson2014). Firstly, the model with the lowest Akaike information criterion metric corrected for small samples sizes (AICc) was inspected, as these are considered to be best in terms of goodness of fit and complexity (Burnham & Anderson Reference Burnham and Anderson2004, Warren & Seifert Reference Warren and Seifert2011). Following this, the threshold-independent metric (area under the curve; AUC), the difference between the test and training AUC (AUCdiff) and the minimum training presence omission rate (ORmtp) were assessed to ensure models could discriminate and prevent overfitting (Anderson & Gonzalez Reference Anderson and Gonzalez2011). Response curves were checked for ecological realism (Guevara et al. Reference Guevara, Gerstner, Kass and Anderson2017), and the most appropriate models were projected onto 2070 climate scenarios.
For comparison between current and future distributions, simple mean ensemble approaches (Palm & Zellner Reference Palm and Zellner1992) 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. Reference Liu, Berry, Dawson and Pearson2005). 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. Reference Hijmans, Phillips, Leathwick, Elith and Hijmans2017) 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. Reference Moilanen, Kujala and Leathwick2009) 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 Reference Moilanen and Kujala2008).
Species weightings are user-defined (Lehtomäki et al. Reference Lehtomäki, Moilanen, Toivonen and Leathwick2016). 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. Reference Kujala, Moilanen, Araújo and Cabeza2013).
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 Reference Moilanen and Kujala2008). 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. Reference Kujala, Moilanen, Araújo and Cabeza2013). Estimating species-specific dispersal metrics requires information on dispersal ability (e.g., Salas et al. Reference Salas, Seamster, Boykin, Harings and Dixon2017), 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. (Reference Kujala, Moilanen, Araújo and Cabeza2013), 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 (AICc), threshold-independent metric (AUCtest), difference between test and training AUC (AUCdiff) and the minimum training presence omission rate (ORmtp).
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 km2.
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 Reference Pearson2006). 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 Reference Davis and Shaw2001, Brooker et al. Reference Brooker, Travis, Clark and Dytham2007). 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 Reference Petford and Alexander2020), and due to anthropogenic habitat modification (Honnay et al. Reference Honnay, Verheyen, Butaye, Jacquemyn, Bossuyt and Hermy2002, Opdam & Wascher Reference Opdam and Wascher2004).
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. Reference Seebacher, White and Franklin2015). 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. Reference Urban, Richardson and Freidenfelds2014, Buckley et al. Reference Buckley, Ehrenberger and Angilletta2015, Wiens et al. Reference Wiens, Litvinenko, Harris and Jezkova2019). The distributions of the species examined are also likely to be impacted by biotic factors that can also affect range edges (Wiens et al. Reference Wiens, Stralberga, Jongsomjita, Howella and Snyderb2009). 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. Reference Carroll, Dunk and Moilanen2010). 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.
Supplementary material
For supplementary material accompanying this paper visit https://doi.org/10.1017/S0376892921000370.
Acknowledgements
We thank Ryan van Huyssteen for his assistance with fieldwork and comments on the manuscript and Ruan Stander for allowing the use of his image of Platysaurus intermedius. We are grateful to all landowners who allowed access onto their property and for their hospitality, particularly Jabu and Bibi Linden from Lajuma Research Centre and Marietjie and Hannes Underhay, Oldrich and Judy van Schalkwyk and the Endangered Wildlife Trust from Medike Nature Reserve. We thank three anonymous reviewers for their comments that improved the manuscript.
Financial support
This research did not receive any specific grant from funding agencies in the public, commercial or not-for-profit sectors.
Conflict of interest
None.
Ethical standards
Ethical approval to record rupicolous reptiles in the Soutpansberg was granted by the Animal Ethics Screening Committee of the University of the Witwatersrand – Clearance Certificate Number 2017/08/53/B.
Open access
