Extinction risks of a Mediterranean neo-endemism complex of mountain vipers triggered by climate change

Climate change is among the most important drivers of biodiversity decline through shift or shrinkage in suitable habitat of species. Mountain vipers of the genus Montivipera are under extreme risk from climate changes given their evolutionary history and geographic distribution. In this study, we divided all Montivipera species into three phylogenetic-geographic Montivipera clades (PGMC; Bornmuelleri, Raddei and Xanthina) and applied an ensemble ecological niche modelling (ENM) approach under different climatic scenarios to assess changes in projected suitable habitats of these species. Based on the predicted range losses, we assessed the projected extinction risk of the species relative to IUCN Red List Criteria. Our result revealed a strong decline in suitable habitats for all PGMCs (63.8%, 79.3% and 96.8% for Xanthina, Raddei and Bornmuelleri, respectively, by 2070 and under 8.5 RCP scenario) with patterns of altitudinal range shifts in response to projected climate change. We found that the mountains close to the Mediterranean Sea are exposed to the highest threats in the future (84.6 ± 9.1 percent range loss). We also revealed that disjunct populations of Montivipera will be additionally highly isolated and fragmented in the future. We argue that leveraging climate niche projections into the risk assessment provides the opportunity to implement IUCN criteria and better assess forthcoming extinction risks of species.

In recent decades, risks of population declines and local extinctions in many species of reptiles and amphibians have been reported from across the globe 1 . A multitude of factors has been identified as potential main drivers behind these declines, including climate warming and climate change as a widely approved factor 2,3 . However, determining the role of climate change on biodiversity decline is challenging as numerous factors are involved and climate changes vary depending on spatial scales 4 . Conservationists have assessed the biological risks from projected impacts of climate changes using empirical species distribution models 2,5-7 , demographic simulations [8][9][10] or phylogenetic proxies [11][12][13][14] .
One of the most obvious effects of climate change on biodiversity is the projected change in species geographical ranges 15,16 . The geographic distribution of any species, with some degrees of disequilibrium, is a spatial representation of its 'ecological niche' . Extrinsic factors such as abiotic conditions and species interactions, and intrinsic ones such as dispersal capability and genetic plasticity determine form, dynamic, and evolution of the niche 17 . The interactions between these factors through past evolutionary time-scales determine the current distribution range of the species 18 . In the Anthropocene, the rapid rate of climate change has caused signatures of range loss and range shift in many species, with an observed time lag of population dynamic responses relative to climatic trends 8 . Over the past century, increasing evidence of climate change impacts on biodiversity patterns has been documenting a globally consistent rate of poleward or upward shifts in many species' geographic distributions 6,7,12,15 , although some species show deviating responses 19,20 .
The consequences of climate change on conservation status and extinction risks have been assessed for a wide array of taxonomic groups. For example, Thomas, et al. 21 estimated that 18-35% of the analyzed species are at risk

Results
The ecological niche models revealed that suitable habitats of mountain vipers, particularly those of the Bornmuelleri and Raddei clades are patchily spread in mountain ranges under current climates. For the Bornmuelleri clade these areas mainly stretched along the Anatolian diagonal and along two mountainous landscape parallel to the Mediterranean Sea, Mount Lebanon and Anti-Lebanon Mountain. Suitable habitats for the Raddei clade were found in the highlands of the Zagros and Alborz mountains within Iran, and in Ararat and Lesser Caucasus in Turkey and Armenia. For the Xanthina clade the predicted suitable habitats were stretched close to coastal lowlands of the Mediterranean Sea in Turkey. The area under the ROC curve (AUC) and true skill statistics (TSS) values of the 10-fold replicated split-sample tests showed excellent prediction accuracy, as all models revealed an average AUC > 0.8 (except the SRE model for Bornmuelleri and Xanthina) and an average TSS > 0.6 among replicates (Fig. 2). Among the four algorithms, GBM presented the highest and SRE the lowest model accuracies, both for AUC and TSS. Highest model accuracies were generally obtained for the Raddei clade. The AUC of the consensus models of the Bornmuelleri, Raddei, and Xanthina clades were 0.925, 0.939 and 0.911, respectively, and the according TSS values were 0.786, 0.825 and 0.733. The consensus distributions of all three Montivipera clades indicated a continuous pattern of range contraction (reduction in projected AOO) under changing climates during the twenty-first century (from the current to 2050 and to 2070) and for RCP scenarios (from RCP 4.5 to RCP 8.5) (see Figs 3,4,5 for the three species and Table 1). The extent of range shifts was quite  www.nature.com/scientificreports www.nature.com/scientificreports/ consistent and similar under the two dispersal scenarios, i.e. no and unlimited dispersal, hence we focused on the no dispersal range changes hereafter. However, the magnitude and direction of range changes vary by clade and emission scenario, with the Bornmuelleri clade showing the highest and the Xanthina group representing the lowest loss in projected AOO (Fig. 6). As we expected, range contractions increase over time and under the more severe emissions scenarios. Accordingly, the extent of projected AOO would get the least reduction by 2050 and using RCP 4.5 (34.6%, 52.1% and 77.6% for the Xanthina, Raddei and Bornmuelleri clades, respectively), while the highest reduction was found under the severe scenario, RCP 8.5 by 2080 (63.8%, 78.3% and 96.8% for Xanthina, Raddei and Bornmuelleri, respectively; Table 1).
The assessment of extinction risks revealed that the Bornmuelleri clade is most at risk (i.e. Threat 2) as it is projected to lose 97% of its AOO by 2070 under RCP 8.5. For the Raddei and Xanthina clades the highest risk was projected to be categorized as Threat 3 (78% and 64% range loss, respectively; Fig. 6) by 2070 under RCP 8.5. Regarding range gain statistics, we also found that the Bornmuelleri clade exhibited a complete absence of gaining newly suitable habitats (percent gain = 0) under any future climate scenario, whilst the Xanthina group was projected to gain the highest proportion of newly suitable habitats (10-11%; Fig. 6).
The results from the landscape analysis showed a decreasing trend in the extent and number of patches projected to remain climatically suitable for all Montivipera clades. While, the decreasing trend in the number of patches for the Xanthina clade was low, the decreasing trend of the average patch size was striking for this group (Fig. 7a). Also, both the patch isolation (i.e. aggregation index) and the patch fragmentation (i.e. splitting index) showed increasing trends from current to projected future climates (Fig. 7a). This was particularly the case for the mountain-dwelling clades Bornmuelleri and Raddei, which in turn were reasonably consistent with their projected altitudinal upward shifts ( Table 2).  www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
Most previous studies have found that the impact of climate change on biodiversity range dynamics emerges through projected range shifts and/or range shrinkage 6,12,22,47 . While range shift can be attributed to poleward latitudinal range displacements, as is a case for species inhabiting extensive geographic areas, the range shrinkage often occurs in range-restricted species, often inhabiting mountains, and results from upward altitudinal range contractions 16,48 . The result of projections under GCMs and emissions scenarios in our study revealed a consistent response of mountain vipers in terms of range contractions, originating from an altitudinal shift to higher elevations (Tables 1 and 2; Fig. 2-6). However, the three clades showed varying degrees of range loss under climate change. The Bornmuelleri clade was projected to lose almost all suitable habitats under the RCP 8.5 by 2070. From a geographic and geomorphologic point of view, the remaining suitable habitats in this clade, which are projected to shift towards a warm Mediterranean climate regime, consist of mountainous landscapes close to flat coastal areas.
With respect to the significant climate shifts projected by 2100 49,50 , the Mediterranean Basin will become an important hotspot for biodiversity conservation 51 . Our results are thus consistent with recent updates on projected species habitat losses due to global warming in the Mediterranean and Middle East 7,52,53 . Nevertheless, the effects of climate change on species range shifts in the eastern Mediterranean has not properly been investigated in the literature, and the few existing studies have only focused on very local patterns 7,54 . A key finding of this study is that there is a strong difference in the climate change impacts among the three species clades inhabiting three different eastern Mediterranean sub-regions. Accordingly, our results highlight a significant sub-regional variation in the impacts of climate change on biodiversity, specifically in wide-ranged species such as mountain vipers. Previous studies on climate change projection over the period 1979-2100 have revealed that among the Middle East -North Africa domain the alpine ecosystems of Turkey and Iran are likely to face the largest increase in temperatures, particularly during March-May 49 . This time of the year coincides with the post-hibernation mating activities for Montivipera species. Accordingly, apart from projected reductions in the mountain vipers' suitable habitats and increasing habitat fragmentation, climate changes might also interfere with the species reproductive cycles, and thus their demographic rates.
While there is a proliferation in studies analyzing the projected impact of climate change on biodiversity range shifts 16,22 , bioclimatic niche projections mostly neglect important drivers of ecological integrity 55 , and often fail to address habitat connectivity, despite being an important element for successful conservation of biodiversity under a changing climate 56,57 . Here, by using metrics of patch configuration, we incorporated the concept of landscape integrity to assess the impact of climate change on the capacity for dispersal and long-term resistance of the species. Furthermore, we eliminated areas that are unsuitable because of human land transformation as extracted from a human footprint dataset. This allowed us to incorporate human-induced barriers into our landscape integrity analysis and increased the realism in projections to the landscape scale. Our approach demonstrated that significant proportions of climatically suitable habitats of mountain vipers are both inaccessible and/or highly www.nature.com/scientificreports www.nature.com/scientificreports/ fragmented due to the dominance of humanized developments. The decline in both the extent and number of suitable habitat patches was consistent among all species, with the Raddei clade showing a lower decrease in the number of patches compared to the Xanthina and Bornmuelleri clades, although its suitable habitat range was highly reduced (Fig. 7a). We thus conclude that the currently suitable habitat patches for the latter two groups will likely become too small for supporting viable populations in the future, as they will lose many of these small marginal patches. For the former, the currently suitable habitat represents core patches that will still remain suitable in the future, albeit losing extent.
The configuration-isolation indices (i.e. AI and SI; Fig. 7a,b) revealed that mountain vipers will be facing an extreme rate of degradation in habitat quality due to increased isolation and fragmentation under projected future climates. The decrease in the aggregation index is caused by increased distance among subpopulations which is the primary factor of reduced connectivity 58 . The impact of climate change on connectivity is more deteriorating for range-restricted species, such as Montivipera species, as their population dynamics is suffering from their strongly-limited dispersal behaviors 55,59 . In addition to isolation, mountain vipers' suitable ranges are increasingly subdivided into smaller patches under projected future warming, which in turn will increase their vulnerability to climate change. Generally, species' vulnerability to environmental stress is thought to depend on three main elements: exposure, sensitivity and adaptive capacity 60 . In the case of mountain vipers, the high magnitude of  www.nature.com/scientificreports www.nature.com/scientificreports/ range loss due to climate change (i.e. exposure 21,61 ) and human-induced barriers to connectivity impose high levels of vulnerability. Mountain vipers are ecological specialists distributed among isolated populations, thus exhibiting a patchy distribution. Their reduced population sizes and life histories exhibiting long generations limit their adaptive capacity to changing climates (i.e. sensitivity 10,62 ), which additionally intensify their vulnerability to changing climates. Finally, ongoing trailing-edge fragmentation, limited dispersal capability, lack of gene flow between disjunct populations and low effective population size due to declining population trends all lead to likely low genetic diversity within populations (i.e. adaptive capacity 22,63,64 ) which adds again to the vulnerability of mountain vipers to climate change. Given the reduction in the projected AOO of the PGMCs we then assigned them to Threat levels (see Methods for more details on Threat levels). Of the three PGMCs, Bornmuelleri is at highest risk to go extinct and would be classified as IUCN Threat 2, but very close to critically threatened (Threat 1) supported by projections of all GCM models. Overall, the uncertainty in range loss values due to the different climate change models that we considered was very low (SD = 3.8-9.2% for Raddei, SD = 1.1-9.7% for Bornmuelleri, SD = 3.2-12.26% for Xanthina) leading to all three groups being assigned into only one threat level. There are additional drivers that could intensify the risk of species extinctions. For example, catastrophic  www.nature.com/scientificreports www.nature.com/scientificreports/ events such as climate extremes, pathogens, large fires or spread of invasive organisms, although occurring at a low or unpredictable frequency, may impose severe consequences and irreversible effects 65 . Increasing human pressure and land use change, which have been intensifying particularly in the eastern Mediterranean Basin and the Middle East 53 are important threats to the connectivity of the Montivipera's remnant populations. In addition to habitat destruction, land use degradation hampers the effective efficiency of mountain vipers to track climate change and hence increases the vulnerability of these species. The species have to disperse to areas that are becoming newly suitable under climate change. If populations cannot disperse at the velocity of environment shifts, the species may only persist in remaining refugia 66,67 . The capability to disperse across landscapes may be further reduced by human developments that dissect suitable habitats into isolated patches, thus reducing habitat connectivity and increasing barriers to dispersal 68 . Altogether, the likely fact that most species will not be able to balance their distributional dynamics with the velocity of environmental changes (e.g. future climate change), highlights the necessity of considering 'assisted migration' as an adaptive conservation strategy 67 . Moreover, regarding the recent decline in natural populations of some of the Montivipera species (e.g. M. latifii and M. wagneri), ex situ conservation, for example captive breeding and introduction-reintroduction programs, could be considered as a strategy to assist the long-term persistence of the species in areas that are highly sensitive to climate changes.
Although EOO and AOO are globally accepted as surrogates of species' geographic range to assess extinction risk, their application to assign species to risk categories have remained controversial 65 due to issues of spatial scale and availability of data 33,36 . The finer the scale at which the distribution of species is mapped, the larger will be the area within which the taxon is recorded to be absent, and the more likely the assessed extinction risks will be higher. Conversely, coarse-resolution mapping of species distributions fails to map unoccupied areas, resulting in lower assessed extinction risks. Moreover, mapping species at coarse resolution increases the likelihood of neglecting environmental heterogeneity and hidden local microclimatic refugia, particularly in mountain areas 5,69 . The assessment of climate change impacts and extinction risks is thus scale-dependent, with finer spatial scales usually being more representative of a species' biology than coarse resolution assessments. In our study, the ensemble models were developed at a relatively fine grid size (ca. 1 km 2 ) through which we included local topographic heterogeneity effects and cryptic climate refugia into the range estimations and risk assessments.
In this study, we performed ENM analyses at a clade level above the species rank rather than at the species level due to the sparse availability of data on species distributions. This procedure assumes that niche conservatism holds for recently diversified species in mountain vipers, which is likely the case in most species that have evolved mostly from allopatric speciation. This speciation mode tends to retain ancestral ecological traits (e.g. ecological niche) 70 . In fact, niche conservatism in climatic tolerance might have limited the geographic range expansion of mountain vipers and may have resulted in very similar species-specific ecological niches in this genus. However, we admit that our approach of phylogenetically lumping occurrence points may have introduced some bias. First, our approach expands geographic range size of the species by assigning them to clades. Given the positive correlation between species' range size and niche breadth 71 , and given that aggregations may smooth over smaller lower-ranked taxonomic units 72 this method may inflate niche breadth and overestimate the range of suitable habitats for species-specific local populations. The resulting habitat suitability maps may thus be rather optimistic. In reality, there will likely be fewer suitable patches available to any species than what is modelled at the clade level, particularly in the absence of unlimited dispersal. Another key assumption of ENMs is that a species is in equilibrium with the climate 67,73 , meaning that the observed distribution represents the total of suitable habitats of the species. Regarding the evolutionary history of Montivipera, high degrees of new-endemism, isolated distribution to mountain and limited dispersal abilities 45,74 might have eradicated the species from some areas that once were suitable. This could mean that Montivipera species might be isolated from otherwise suitable habitats, thus violating to some degree the equilibrium assumption. If true, this would underestimate the true width of suitable habitats, and thus, counteract the overestimation imposed by aggregating species to clades. We believe that both effects are considerably small given the very clear and specific niche identified by our models.
Combining results of population demographic models (e.g. population viability analysis) with distribution models allows for a more accurate prediction of a species extinction risk due to a more detailed representation of the demographic processes 75 . Due to the lack of demographic data, we were not able to construct more dynamic models (e.g. mechanistic models), and we focused on projections of AOO as a proxy of effective habitat loss. Our results thus need to be considered with caution. However, our model-derived estimates of AOO could be useful to design conservation areas 76 , spatial priority-setting implications and conservation planning of the highly endangered mountain vipers.

Methods study area and Montivipera occurrences. The study area encompasses parts of the eastern
Mediterranean Basin, extends to the Near and Middle East and contains two biodiversity hotspots, the Mediterranean Basin and the Caucasus, which are both rich in endemics from a faunal and also from a floristic viewpoint. From a geomorphologic viewpoint, the study area is predominantly mountainous with a small coastal fringe in the west to an increasing alpine system stretching towards the east, starting from the Taurus Mountain in Turkey and ending with the Lesser Caucasus in Armenia and with the Alborz and Zagros Mountains in Iran. Although this region mainly receives its humidity and rainfall from the Mediterranean Sea, long dry seasons have caused the formation of a heterogeneous mixture of sclerophyllous cover types spanning from alpine steppe grasslands to mixed and deciduous wood-forest belts dominated by oaks, junipers, and pines.
Species occurrences used as input to ENM were obtained from field sampling between 2011 and 2016 (particularly for species of the Raddei complex, n = 84) and from other herpetologists' field studies (n = 62), extractions from existing databases (i.e. Herp-Net and GBIF, n = 16), and from scientific publications (n = 10). In order to avoid effects from spatial autocorrelation of the presence points due to different protocols of field sampling 77 Scientific RepoRts | (2019) 9:6332 | https://doi.org/10.1038/s41598-019-42792-9 www.nature.com/scientificreports www.nature.com/scientificreports/ we randomly sub-selected one point per 10 km radius buffer using the sp package 78 in R v. 3.4.3 79 . In summary, we retained 152 records ranging from 7 for M. b. albizona to 37 for M. xanthina.
Mountain vipers have a highly isolated distribution with incompletely sampled data. Hence, the insufficient data of the species occurrence is too sparse for performing ENMs at the species level. By lumping occurrence points of the closely related sister-clades we developed niche models above species level 80 . To do so, using a recently constructed phylogenetic tree of the genus 45 (see Fig. 1), which subdivides the genus into geographically well-separated clades, we divided the current Mediterranean distribution of the genus into three phylogeographic Montivipera clades (PGMC). These PGMCs represent three different bioregions including coastal habitats (western Turkey), nearshore mountain habitats (southern Turkey, Lebanon, and Syria) and inland mountainous habitats (northern and western Iran and southern Caucasus), which encompass the core area of the distribution of the three phylogenetic clades Xanthina, Bornmuelleri, and Raddei, respectively (Fig. 1a,b). The ENM analyses were performed such that all occurrences were pooled by PGMCs for further modelling in order to have sufficient observations to build models, obtaining 37, 80 and 35 occurrence points for Bornmuelleri, Raddei and Xanthina clades, respectively. However, we removed the Taurus occurrences of the Xanthina sub-clade because of its more mountainous distribution compared to the rest of the Xanthina sub-clade (see Fig. 1c). ecological niche modelling approach and climate change. To track the impact of climate change on the distribution of Montivipera clades we conducted ENM based on climatic variables at a ~1 km resolution derived from WorldClim 81 . We only used climate variables rather than other vegetation-derived predictors in the ENM analysis in order to assess changes of the climatic suitability and to avoid circular reasoning, as vegetation properties themselves depend on climate. In order to sample the whole study area for background contrast, we allocated 10,000 pseudo-absences spatially at random. We then extracted the values of all 19 bioclim variables for presence and pseudo-absence points and calculated the degree of multicollinearity among them based on the variance inflation factor (VIF). To do so, we used the 'usdm' package 82 and set a VIF value of 6 and a correlation threshold of 0.75, as recommended by Guisan, et al. 83 . As a result, we retained 6 bioclim variables for modelling, including annual mean temperature (bio1), temperature seasonality (bio4), maximum temperature in the warmest month (bio5), annual precipitation (bio12), precipitation of driest month (bio14) and precipitation seasonality (bio15).
To assess the habitat suitability of the species, we fitted four ENM methods, including generalized linear models (GLM), generalized boosting models (GBM), surface range envelop (SRE, or commonly called as BIOCLIM), and maximum entropy (MaxEnt) and combined them into a final ensemble model using biomod2 package 84 . Using the four methods allowed us to combine simple regression (GLM) and envelop (BIOCLIM) models and complex machine learning (MaxEnt, GBM) models 85 and thus to assess uncertainties arising from the choice of methods. For GLM we used simple and quadratic terms and stepwise selection procedure based on Akaike Information Criteria (AIC) to fit the best model. GBM was set to allow up to 2500 trees, and we set the learning rate to 0.001 and the bag fraction to 0.5. To conduct SRE the species envelop was set up with a quantile 0.05% to avoid the most extreme values. For MaxEnt all feature types (linear, quadratic, product, threshold and hinge) were allowed and a maximum iteration of 200 was used. For all models we weighted presences and pseudo-absences inverse-proportional so as to give them equal weights and the models' predictive performance was evaluated by means of a repeated (10×) split-sample test (75% training, 25% evaluation) for which we calculated the AUC and the TSS. The final ensemble model was obtained by weighted-averaging the individual models proportionally to their AUC scores 84 .
The future habitat suitability of the clades was estimated by projecting fitted models to four different global circulation models (GCM), including CCSM4 86 , CNRM-CM5 87 , HadGEM2-ES 88 and MIROC5 89 for two time periods of 2050 (average for 2041-2060) and 2070 (average for 2061-2080). Using four different GCMs allowed us to assess the uncertainty from selecting GCMs as these differ clearly among regions 90 . We included a moderate and an extreme greenhouse gas emission scenario in our analyses 91 . For each GCM, we considered two representative concentration pathways (RCPs), namely RCP 4.5 and RCP 8.5. To obtain a final ENM projection, we averaged the ensemble rasters arising from projecting the ENMs to the four GCMs for both scenarios and for the two time periods. Data for bioclim variables based on the four GCMs scenarios and time periods were downloaded from WorldClim at a ~1 km resolution. For each PGMC, we therefore implemented 640 projections (4 ENM algorithms × 10 replications × 4 GCMs × 2 RCPs × 2 time-periods).
To track changes in suitable habitats we calculated the percentage of cells that gained or lost climatic suitability for the four GCM models, for the two RCPs, and for 2050 and 2070, compared to the current area of suitable habitat for each PGMC. To identify suitable cells from unsuitable ones, we first converted continuous ensemble models into binary presence/absence maps based on the minimum habitat suitability score at occurrence points per PGMC as a threshold. Using this approach, we calculated the Minimal Predicted Area (MPA), a minimum area where the species is projected to exist 92 , for establishing the potential AOO for all clades. For this step, we assumed that all current occurrence points represent suitable habitat condition for the species and we did not remove marginal occurrences as "unsuitable". Then, based on the number of pixels gained or lost under future climatic scenarios we calculated the respective changes in projected AOO (i.e. range shifts) of the groups. We analysed shifts in suitable habitats by employing two different dispersal scenarios, i.e. unlimited dispersal and no dispersal as two extremes (for more details on dispersal scenarios see reference number 93 ).
To obtain an accurate prediction of suitable habitats and to avoid inaccessible areas being identified as suitable for mountain vipers, such as developed areas, we masked these human transformed areas from climatically suitable habitats. Accordingly, we also incorporated the human landscape context into suitability projections. We identified developed areas based on the human footprint model developed and implemented by Sanderson et al. (2002) 94 . This layer incorporates data on population density and the presence of infrastructures including road networks, land transformation and human access with pixel values ranging from zero (undeveloped) to 100 (cities and highly developed areas). For our study, we assigned areas with values greater than 50 as developed areas and (2019) 9:6332 | https://doi.org/10.1038/s41598-019-42792-9 www.nature.com/scientificreports www.nature.com/scientificreports/ masked them from the predicted habitat suitability maps. We admit that these humanized areas will expand in the future. Yet, due to the lack of available data for future conditions, the mask was applied unchanged to the future climate scenarios. However, urban growth and change in developed areas might not severely alter the results given the large spatial scale of our analyses.
Next, we assessed the severity of future extinction risks of the clades by categorizing them into four threat levels using IUCN Red List Criteria following Maiorano, et al. 53 . Threat levels are assigned based on the projected AOO losses, with 1: 100%, 2: <100% and > = 80%, 3: <80% and > = 50%, and 4: <50% and > = 30%. Finally, we assessed the spatial configuration of the suitable habitats of the clades by borrowing landscape ecology metrics that express facilitation or suppression of dispersal and long-term persistence of the clades. To do so, we calculated the number of patches (NP), the proportion of the landscape (PL), the aggregation index (AI), and the splitting index (SI) of suitable habitat patches under current and projected future climates. AI is a metric of spatial adjacency and demonstrates the frequency at which different pairs of patches appear side-by-side in a landscape 95 . The more clumped and compact the patches, the higher are the values of AI. SI measures the degree of fragmentation and can be interpreted as the number of patches with a constant size when the landscape is subdivided into S patches, where S is the value of the splitting index 95 . SI increases as the landscape is increasingly subdivided into smaller patches. All metrics were calculated using the 'SDMTools' R package. A more detailed description of the metrics is available in McGarigal, et al. 95 .

Data Availability
The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.