Assessing future distribution, suitability of corridors and efficiency of protected areas to conserve vulnerable ungulates under climate change

Central part of Iran harbours populations of wild ungulates that are threatened or extinct over large parts of the region, and are likely to be impacted by climate change. In this study, we predicted the impact of climate change on the distribution of three vulnerable ungulates in central Iran. We then evaluated future suitability of corridors connecting the protected areas for movement of the ungulates in response to climate change.


| INTRODUC TI ON
Climate change is now recognized as one of the main threats to biodiversity (IPCC, 2013). It has been estimated that as a result of climate change, nearly a quarter of all species may face increased extinction risk, among which are many large mammals (Thomas et al., 2004).
Large mammals (adult body weight > 5 kg; Bourliere, 1975) are already at high risk of extinction as their intrinsic traits such as low population density, low reproduction rate, low life history and large body size (Fisher & Owens, 2004) make them highly vulnerable to anthropogenic factors including habitat loss, poaching and anthropogenic land uses. Synergistically acting with these factors, climate change could also put large mammals at greater risk of extinction in the future (Cardillo et al., 2005). These impacts may be more challenging when climate change forces large mammals to shift their ranges. Although their large body size enables large mammals to readily reach suitable habitats (Schloss, Nuñez, & Lawler, 2012), habitat fragmentation and reduced dispersal ability as a result of anthropogenic land uses could substantially hamper their range shifts in response to climate change (Robillard, Coristine, Soares, & Kerr, 2015). Therefore, to effectively conserve large mammals under climate change, high priority should be given to conservation plans aiming at incorporating their future distributions and enhancing landscape connectivity for their populations (Groves et al., 2012;Hannah et al., 2007).
Evaluating the impact of climate change on distribution of species is mainly accomplished by applying species distribution models (SDMs; Araújo & Peterson, 2008;Pearson & Dawson, 2003).
Application of SDMs to predict the future distribution of species is subjected to several limitations as they do not include other factors limiting their distribution such as competition, dispersal and predation (Pearson & Dawson, 2004); however, they still remain the most widely used tools for this purpose. SDMs find the statistical relationship between a set of environmental variables (i.e. climate variables) and species geographic distribution which is then applied to different scenarios of greenhouse gas emission to predict climatically suitable habitats in the future (Pacifici et al., 2015). From future predictions, it is then possible to evaluate changes in the extent of species distribution by estimating proportion of habitats predicted to remain suitable and become unsuitable in the future. SDMs also enable us to evaluate the impact of climate change on protected areas and identify those resilient to climate change (Araújo, Alagador, Cabeza, Nogués-Bravo, & Thuiller, 2011;Carroll, Dunk, & Moilanen, 2010).
The feasibility of species range shifts in response to climate change largely depends on the availability and accessibility of suitable habitats in the landscape (Warren et al., 2001). That is why increasing landscape connectivity-the degree to which it facilitates or hampers species movements-is one of the most important adaptive strategies to conserve species under climate change (Watson, Iwamura, & Butt, 2013). When highly suitable corridors are identified and maintained, species will have the chance to escape from unsuitable conditions resulting in increasing their resiliency to climate change (Hannah, 2011;Morecroft, Crick, Duffield, & Macgregor, 2012). Evaluating landscape connectivity and identifying corridors could be done by employing different approaches, one of the most common of which is using SDMs (Littlefield, Krosby, Michalak, & Lawler, 2019). When applied in the context of climate change, using SDMs offers three main advantages including: (a) identifying opportunities to maintain landscape connectivity based on current distribution of suitable habitats (i.e. stable habitats); (b) considering species-specific needs for connectivity under climate change with respect to the changes in distribution and availability of suitable habitats in the future (Littlefield et al., 2019); and (c) incorporating anthropogenic impacts into future predictions of species distribution, and identifying pathways assisting their future movements through the least disturbed parts of the landscape.
Across western Asia, Iran is regarded as one of the most important countries for future conservation of threatened wild ungulates.
At global scale, Iran includes considerably large extent of suitable habitats for populations of goitered gazelle (Gazella subgutturosa).
Until 2003, goitered gazelle was classified as a "near-threatened" species in IUCN red list; however, large reductions in population size caused this species to move into the "vulnerable" category (IUCN SSC Antelope Specialist Group, 2017). More topographically rich areas in the country harbour large populations of two vulnerable mountain ungulates (Valdez, 2008;Weinberg et al., 2008), the wild sheep (Ovis spp) and wild goat (Capra aegagrus) for which Iran encompasses major portions of their global populations. These three ungulates used to be largely distributed across the country; however, habitat degradation and poaching caused their populations to suffer from substantial reductions in distribution and population size (Bashari & Hemami, 2013;Esfandabad, Karami, Hemami, Riazi, & Sadough, 2010;Hemami & Groves, 2001). Currently, the majority K E Y W O R D S CIRCUITSCAPE, corridors, ensemble modelling, landscape connectivity, latitudinal shifts, protected areas, vulnerable ungulates of these ungulate populations live within protected areas, specifically those in central part of the country (Ansari, 2016;Esfandabad et al., 2010;Hemami & Groves, 2001). Central Iran has a large network of protected areas (PAs) and no-hunting areas (NHAs) with high potential to protect habitats and populations of the ungulates.
However, the efficiency of this protected network and likely changes in suitability of protected habitats for these mammals has not been evaluated yet under climate change. The level of protection and restriction on anthropogenic activities differ within the central Iranian network of protected areas. Therefore, the efficiency of some of these sites under climate change may largely depend on the future management strategies implemented. This particularly holds true for protected areas with high future conservation capability, but low protection level, which will need more management considerations. Furthermore, considering the key role of the protected areas in survival of the ungulates in Iran, predicting future distribution of these mammals to inform future expansion of the current network of protected areas or selection of new sites could be an effective conservation strategy under climate change. More importantly, due to the expansion of human land uses, efficiency of the network of protected areas could be largely improved by maintaining key corridors that facilitate ungulate movements in response to climate change.
In this study, employing species distribution and connectivity modelling approaches, we predicted the impact of future climate change on geographic distribution and corridors for the three vulnerable ungulates of goitered gazelle, wild sheep and wild goat in dry and semi-dry regions of central Iran. Our aims were to the following: (a) predict how climate change would impact on future distribution of the target ungulates in central Iran; (b) identify highly suitable corridors for the ungulates likely to persist under climate change that could assist their future movements in response to climate change; and (c) evaluate future changes in the extent of suitable habitats within the protected areas and identify those sites with remained efficiency under climate change. Additionally, for the three ungulates, we hypothesized that as a result of climate change, their geographic distributions would display significant shifts towards higher elevations and latitudes in central Iran.

| Study area
The study area encompasses the central part of Iran within an extent of about 547,000 km 2 and elevational gradient of 518-4300 m ( Figure 1).
The region exhibits high diversity of topographic features which include a mixture of flat plains, rolling hills, mountain ranges and isolated mountains. Low relief areas are the dominant topographic features in F I G U R E 1 Geographic location of the study area in Iran and the network of protected areas in the region (map source: www. doe.ir). (See Table S1 for protection category and protected area code.) central Iran occupying mainly the eastern and southern parts. In contrast, western and the south-western parts of the landscape are topo-

| Collection and preparation of occurrence data
Efforts to collect occurrence data for the three ungulate species were mainly based on extensive fieldwork conducted during 2016-2018. As parts of the occurrence data had been previously collected by environmental guards (field observations recorded during annual total counts of ungulates by Departments of Environment (DoE) across Iran), we first evaluated the existing datasets to identify data gaps and conducted further fieldwork to complement the datasets (Table S1). Due to the very low density of ungulates outside the protected areas (and general absence of goitered gazelle in eastern parts), we mainly focused the field observations within the protected areas for which we obtained most of the observations. However, fieldwork was also extended beyond the sites across the unprotected landscape where there have been reports of observing target ungulates (i.e. wild sheep and wild goat). During fieldwork, occurrence data were obtained from direct observations of individuals and indirect sources including scat identification. In total, we complied 167, 593 and 697 occurrence points for goitered gazelle, wild goat and wild sheep, respectively. We then checked the datasets for spatial autocorrelation and applied the global Moran's test to evaluate the structural pattern of the occurrence data. We then extracted three spatially independent subsets of the occurrence data by setting minimum interval distance of 5 km between all the occurrence points. This resulted in a total of 87, 235 and 224 spatially independent occurrence points for goitered gazelle, wild sheep and wild goat, respectively, which were used in the following analyses ( Figure S1).

| Modelling species distributions
To predict distribution of the ungulates, we selected 26 environmental variables belonging to four categories of topography, bioclimate, vegetation and human impact (Table S2). For climate, we used 19 bioclimatic variables at 30-s spatial resolution from CHELSA database (Karger et al., 2017;CHELSA-climate.org). To describe the topography of the study landscape, we selected four variables of elevation, slope, aspect and ruggedness. The digital elevation model (DEM) of the study area was downloaded from the USGS database (USGS.org) at 30-m resolution which was resampled to 1-km resolution and used to produce rasters of slope and aspect using surface analysis tool in ArcGIS 10.3 (ESRI, 2014). For ruggedness, we calculated vector of ruggedness index (VRM) (Sappington, Longshore, & Thomson, 2007) using the DEM of the study area and 5 × 5 km moving window. We considered two vegetation-related variables of soil-adjusted vegetation index (SAVI) and vegetation cover in developing the distribution models. Important classes of vegetation cover for each ungulate were extracted from rangeland and forest map of the study area (www.frw.ir), and the associated map was converted to a raster of continuous values by calculating the proportion of each class within grids of 5 × 5 km using neighbourhood analysis in ArcMap. We also took the same approach to prepare maps of human impact using two land use categories of human settlements and cultivated lands extracted from rangeland and forest map of the study area. After preparing layers of the environmental variables, they were evaluated for collinearity using Spearman correlation coefficient and based on the threshold value of 0.8 (Elith et al., 2006). To include the most important variables into distribution modelling, we first inputted all the variables into the MaxEnt model (Phillips & Dudík, 2008) and ranked them based on the result of jacknife analysis. Then, among the most important linearly correlated variables, final variables were selected with respect to their ranking in predicting distribution of the target species (Table S2).
Potential current and future distributions of the target ungulates were predicted employing an ensemble modelling approach. For this reason, we selected five SDMs with high predictive performance including generalized linear model (GLM), generalized boosting model (GBM), random forest (RF), maximum entropy (MaxEnt) and multivariate adaptive regression splines (MARS) to develop the ensemble models (hereafter EMs) using BIOMOD2 package (Thuiller, Lafourcade, Engler, & Araújo, 2009) implemented in R 4.3.0 software (R Development Core Team, 2014). As a requirement for BIOMOD modelling framework, three sets of 5,000 background points were randomly selected for each species from the entire landscape of central Iran as pseudo-absence points. We used 75% of the occurrence data to fit the individual SDMs and used the remaining 25% to evaluate their predictive performance. The predictive performance of SDMs was evaluated using threshold-independent indices of area under the curve (AUC) of a receiver operating characteristic (ROC; DeLeo, 1993) and the true skill statistic (TSS; Allouche, Tsoar, & Kadmon, 2006). To obtain the consensus predictions, we used a weighted-averaging approach (Thuiller et al., 2009) & Rowell, 2015). In total, we developed five future EMs (1 future scenario × 5 GSMs) for each ungulate which were averaged, and a final EM representing the distribution of each species in 2070 was obtained. In the next stage, current and future EMs were converted into binary maps using the minimum suitability score at occurrence points and then overlaid. Accordingly, we quantified three indices of climate change impact on distribution of the ungulates including lost, gained and stable habitats. Overlaying these maps with the map of protected areas, we also evaluated future changes in the extent of suitable habitats within the individual sites and the whole protected network for each species. We also evaluated whether climate change would result in distribution of the ungulates to shift along elevational and latitudinal gradients. For elevational shift, we obtained mean elevation of current and future distributions for each species and calculated difference between them (Luo, Jiang, & Tang, 2015).
For latitudinal shift, centroids of the current and future distributions were identified, and the geographic distance between pairs of centroids was taken as a measure of latitudinal shift (Luo et al., 2015).

| Modelling Landscape connectivity and identifying movement corridors
To predict the most permeable corridors for the three ungulates, we   Predicted changes between current and future are depicted in yellow (suitable habitats predicted to become unsuitable) and pink (unsuitable habitats predicted to become suitable). Areas depicted in green represent suitable habitats predicted to remain suitable between current and future conditions. (See Table S1     were predicted to lose larger proportions of their distribution compared with goitered gazelle. This is because mountainous regions have been predicted to become warmer at a higher rate than low-lying areas (Rangwala, Sinsky, & Miller, 2013), a trend that could result in more widespread negative changes in the distribution of mountain species in relation to their lower elevational counterparts.

| Climate change and future distribution of the vulnerable ungulates in central Iran
Comparing current and future predictions showed that as a result of climate change, current distribution of the three ungulates is likely to change along both elevational and latitudinal gradients. Along

F I G U R E 3
High suitability corridors predicted based on the density of cumulative current for movement of goitered gazelle (a), wild sheep (b) and wild goat (c) between the protected areas in central Iran. Areas depicted in yellow and red have the highest density of flowing cumulative current, and represent areas with the highest probability of movement. (See Table S1 for protection category and protected area code.) However, a reduced area of potential habitats with increasing elevation will hamper future upward shifts as these habitats may not be large enough to sustain viable populations. Our models also did not predict any climate-related northern range expansion for the three ungulates. However, because of the stability of northern habitats and the large habitat loss in southern parts, the estimated change could be regarded as northward shifts which are expected to occur to compensate limited opportunities for elevational movements.

F I G U R E 4
Future changes in suitability of modelled corridors for goitered gazelle (a), wild sheep (b) and wild goat (c) predicted based on the RCP 8.5 emission in 2070. Areas depicted in yellow and pink represent connectivity areas predicted to become unsuitable and remain suitable under future conditions, respectively. (See Table S1 for protection category and protected area code.) In this study, our inference on the impact of climate change on the distribution of the target ungulates was obtained based on the consensus predictions of several SDMs. As the results of model evaluation demonstrated, all employed SDMs represented high performance in predicting the distribution of the three ungulates in our study. However, owning to different assumptions of these models, we could not say for sure which model provides the most accurate prediction of species distributions in the future (Thuiller, 2004). Thus, we integrated the projections produced by the individual SDMs through the ensemble modelling approach as it considers a range of possible predictions, which tends to reduce the associated uncertainties. In another conservative approach, we used minimum suitability value at occurrence points to produce binary distribution maps associated with current and future conditions. Using this threshold allows us to reduce the risk of misidentifying additional suitable habitats when planning for the future conservation of these ungulates (Sony, Sen, Kumar, Sen, & Jayahari, 2018). However, it may come with the cost of including habitats of lower suitability which entails explicitly assessing the reliability of these predictions through field investigations before being used in future conservation plans.
Capturing the full realized niche of species is a key assumption when projecting SDMs in space and time (Guisan, Thuiller, & Zimmermann, 2017). In this study, the SDMs developed for the tar-

| Combined role of protected areas and corridors in increasing adaptation of the ungulates to climate change
Using SDMs as the basis for connectivity modelling enabled us to evaluate the potential of the study landscape in facilitating future For wild sheep and wild goat, we identified impacted populations from all three vulnerability groups. Across the study region, there is more concern about the future impact of climate change on populations inhabiting low relief areas of the east and the south-east for two reasons: (a) existence of mountains with small elevational gradient within the protected areas that could limit future upward movements of the ungulates, and (b) low topographic heterogeneity and therefore low permeability of the surrounding landscape that may largely hamper species movements towards key protected areas in the future (Malakoutikhah, Fakheran, Hemami, Tarkesh, & Senn, 2018). As a result of these landscape characteristics, populations inhabiting these parts are expected to be particularly vulnerable to climate change. Across eastern deserts, Abbasabad WR (34) was predicted to retain the largest stable habitats (nearly 600 km 2 ) in the future, as it includes the highest mountainous areas among the protected areas in this part.
Accordingly, this WR could be considered as an important refugia buffering populations of wild sheep and wild goat under climate change. In this part of the landscape, the only corridor predicted to remain suitable is the one connecting Kharoo (31) and Yakhab (29) NHAs. This corridor is highly crucial for the populations in these NHAs, as it provides the only pathway towards key protected areas in the north (both species) and the west (only for wild sheep). In this region, Yakhab NHA (29) could largely contribute in facilitating and maintaining this connection, and thus function as a stepping stone for both species.
Across the west and the north-west, we also identified vulnerable populations of wild goat and wild sheep to climate change; however, high level of topographic heterogeneity across these parts could largely help buffering the negative effects of climate change on these populations (Gavin et al., 2014;Loarie et al., 2009).
Protected populations of wild sheep and wild goat in the west and the north-west were identified as belonging to the second group of vulnerability; that is, populations predicted to lose habitats, but remain connected to other suitable habitats. We predicted, on average, stability of 28.5% of currently suitable habitats for both species in the future. However, with respect to high topographic heterogeneity of the landscape in these parts, the area of available habitats in the future may be larger than what was estimated by our predictions. This is because high topographic heterogeneity of these areas could contribute to climatic heterogeneity and creation of microclimatic refugia (Ford, Ettinger, Lundquist, Raleigh, & Lambers, 2013

| Conservation implications
Iran has a relatively well-established network of protected areas areas. This is particularly the case for roads, a major issue for conservation of large mammals in Iran (Moqanaki & Cushman, 2017). As our results showed, in some areas the predicted stable corridors are already impacted by roads. Many populations of these ungulates occur inside protected areas away from roads and their impacts. However, these structures could be a serious threat in the future by interrupting ungulate movements when tracking suitable habitats or increase their risk of mortality through collision. Therefore, road networks and their current status/ future development with respect to the future need of these mammals should be seriously considered in future conservation programs.
Our study was a case example demonstrating that climate change could have significant impact on our target ungulate species across a proportion of their global range. Accordingly, further studies could provide better insights into the impact of climate change on these mammals, which could be used to inform conservation programs and assist in evaluating their conservation status. Based on our results, we suggest three main conservation strategies for the conservation of the three ungulates, which may also be applicable for their populations in other regions. These include the following: (a) enhancing protection status of the NHAs, PAs and WRs predicted to maintain efficient in future and those likely to function as stepping stones; (b) expanding the network of protected areas through enlargement/ establishment of new sites in areas predicted to remain potentially suitable (protecting more heterogeneous habitats for wild sheep and wild goat, and larger suitable areas for goitered gazelle); and (c) preventing future expansion of anthropogenic land uses across areas with high potential to facilitate future movements of the ungulates, particularly for goitered gazelle as its movement corridors have a higher chance of being disrupted by human land uses.

ACK N OWLED G EM ENTS
We would like to greatly thank the authorities in Departments of Environment within provinces for coordination and cooperation in conducting fieldwork. We are also grateful to the Swiss Federal Research Institute (WSL) for their support and the Iran National Science Foundation, Presidency of Islamic Republic of Iran, for providing financial support to this research (project number: 95849735).

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.