Design of Protected Area by Tracking and Excluding the Effects of Climate and Landscape Change: A Case Study Using Neurergus derjugini

This study aimed to use the applications of Ensemble Species Distribution Modelling (eSDM), Geographical Information Systems (GISs), and Multi-Criteria Decision Analysis (MCDA) for the design of a protected area (PA) for the critically endangered yellow-spotted mountain newt, Neurergus derjugini, by tracking and excluding the effects of climate and landscape changes in western Iran and northeastern Iraq. Potential recent and future distributions (2050 and 2070) were reconstructed by eSDM using eight algorithms with MRI-CGCM3 and CCSM4 models. The GIS-based MCDA siting procedure was followed inside habitats with high eSDM suitability by eliminating the main roads, cities, high village density, dams, poor vegetation, low stream density, agricultural lands and high ridge density. Then, within the remaining relevant areas, 10 polygons were created as “nominations” for PAs (NPAs). Finally, for 10 different NPAs, the suitability score was ranked based on ratings and weights (analytical hierarchy process) of the number of newt localities, NPA connectivity, NPA shape, NPA habitat suitability in 2070, NPA size, genetic diversity, village density and distance to nearest PAs, cities, and main roads. This research could serve as a modern realistic approach for environmental management to plan conservation areas using a cost-effective and affordable technique.


Introduction
Amphibians are valuable components of biodiversity that have declined significantly in many regions of the world [1][2][3]. Extrinsic reasons that have contributed to these declines include several biotic and abiotic factors that act synergistically, such as habitat loss and fragmentation [4], chemical pollutions [5], UV radiation [6], alien species [7], direct exploitation [8], disease [9] and climate change [10]. For all of these reasons, the in situ conservation of viable populations in natural ecosystems is an essential strategy for managing and monitoring amphibian communities [11]. However, amphibians are often neglected for conservation efforts, and only ≤5% of 42% of threatened species are found in the protected areas (PAs hereafter) [12,13]. Although, the establishment and management of global PAs to conserve biodiversity at the lowest possible cost are accepted according to conservation policies [14].
A review of various studies has shown that climate change [15][16][17][18][19][20] and landscape change [21][22][23][24][25][26][27][28] are two critical factors that affect PAs but, to date, few studies have specifically considered how to design a new PA by excluding these emerging threats. Changes in habitat extent, as well as landscape and habitat structure due to human activities such as urbanization [29,30], intensive agriculture [31,32] and roads [33,34], have negative consequences on the abundance and distribution of amphibians [35]. Landscape fragmentation reduces connectivity and genetic diversity among populations [36], which is followed by decreasing fitness and increased risk of local extinction [37].

Environmental Data
To reconstruct recent and future (2050 and 2070) climate projections under optimistic (RCP26) and pessimistic (RCP85) scenarios, 19 bioclimatic variables with 30 s spatial resolution were downloaded from the WorldClim-Global Climate data (https://www. worldclim.org; accessed on 5 March 2021). Climatic data for 2050 and 2070 were derived from two atmospheric circulation models (ACMs): Community Climate System Model Version 4 (CCSM4; see Gent et al. [69]) and Meteorological Research Institute CGCM Version 3 (MRI-CGCM3; see Yukimoto et al. [70]). The CCSM4 is widely used for predicting climate change effects on fauna distribution in Iran, for example, [71][72][73][74]. The MRI-CGCM3 has also been suggested as one of the best models from 37 Coupled Model Intercomparison Project phase 5 (CMIP5) General Circulation Models (GCMs) of regional climate change projections in Iran [75]. Due to the high correlation between bioclimatic variables and preventing any potential collinearity problems, the Pearson correlations for 19 climatic variables were evaluated using SPSS software (v 22.0), and only variables with a correlation less than (r < 0.75) were used for analysis. The models were eventually run using six bioclimatic variables: annual mean temperature (BIO1); mean diurnal range (mean of monthly (max temp-min temp)) (BIO2); temperature seasonality (standard deviation ×100) (BIO4); annual precipitation (BIO12); precipitation of wettest month (BIO13) and precipitation of driest month (BIO14).

Ensemble Species Distribution Modelling (eSDM)
BIOMOD2 package in R v 4.0.30 was used for Ensemble Species Distribution Modelling (eSDM) [76]. Eight algorithms were run: Generalized Linear Models (GLMs), Generalized Boosted Models (GBMs), Random Forest (RF), Classification Tree Analysis (CTA), Artificial Neural Networks (ANNs), Surface Range Envelopes (SREs), Flexible Discriminant Analysis (FDA) and Multivariate Adaptive Regression Splines (MARSs). These models are based on the presence-absence algorithms, and since the absence records were not available, pseudo-absence records, with a number equal to the records of presence, were randomly generated for each model [77,78]. For each model, 80% (training set) of the data were randomly assigned for model calibration and 20% (validation set) for the performance of the algorithms [79]. Every model algorithm was run ten times to eliminate bias caused by the splitting of the total records [80]. Each model was evaluated by the true skill statistic (TSS), receiver operating characteristic curve (ROC or = AUC), and Cohen's Kappa (KAPPA) metrics (see Supplementary Materials Scripts S1 for more details). However, TSS ≥ 0.8 was kept to build the final ensemble [79]. Species range change (SRC) was also used to quantify and represent geographical change over time (Scripts S1) [79].

Data Type and Sources
Iran and Iraq digital elevation model (DEM) with a 30 s spatial resolution (~1 km) was downloaded from WorldClim. An Iran land cover map, including data of urban, water, forest, farmland, rangeland, shrubland, uncovered plain (bareland) with a spatial resolution of 10 m, was downloaded from Google Earth Engine (GEE) (https://earthengine. google.com; accessed on 5 March 2021) [81]. Iran's main roads, villages and protected areas with a 30 m resolution were obtained from the Iranian Forest, Rangeland, and Watershed Management Organization [82]. Iran and Iraq stream layers were created by the DEM layer in GIS. Data of Iraq villages were created by Google Earth Pro. Iraq land cover and the main rod with a 30 s spatial resolution (~1 km) were downloaded from the United States Geological Survey (USGS: https://www.usgs.gov; accessed on 5 March 2021) and Diva-GIS (http://www.diva-gis.org; accessed on 5 March 2021).

Geographic Information System (GIS)
The siting procedure was followed in three stages. After obtaining the eSDM results, in the first stage, based on habitat suitability (HS hereafter), N. derjugini habitats are classified into four classes including habitats with extremely high suitability (75-100%), high suitability (50-75%), medium suitability (25-50%) and low suitability (0-25%) for the present and the years 2050 and 2070. For the current habitat suitability, a polygon has been created around habitats with high and very high (50 to 100%) suitability ( Figure 3). This polygon was used as a study area for the multi-criteria decision analysis (MCDA) based on the GIS siting procedure (Figure 3).
The exclusionary criteria were then applied. The exclusionary criteria were considered as significant factors and were assumed to damage the abundance of N. derjugini if they are within or close to a PA. Supplementary Materials Table S2 listed the exclusionary criteria as well as a brief explanation of why they are important to choose. Therefore, within the study area, all main roads, all cities, all dams, areas with poor vegetation, farmlands (concretely irrigated lands), and elevations above 2500 m were excluded. The streams were classified into five classes with very low, low, medium, high, and very high densities and the very low density areas of the streams (≤0.03 km 2 ) were excluded from the study area. The villages were also classified into two classes with low and high densities and the high density areas of the villages (≥0.72 km 2 ) were excluded from the study area. After deleting the exclusionary criteria from the study area, 10 polygons were created in the remaining areas, especially in the areas where N. derjugini were present. These polygons were considered as "nominations" for the protected areas (NPAs hereafter) in the MCDA analysis. It should be noted that dry farming and gardens are the key activities of the people in these regions, and they covered the majority of the study areas. As a result, although they are depicted in Figure 4, they were not considered when designing the PA and only irrigated agricultural lands were taken into account.   RCP26-2070, genetic diversity, distance to nearest PAs, villages density, distance to cities and distance to the main roads were used in the process of MCDA analysis.
NPA connectivity was measured using the nearest neighbor distance (km) [83]. The size of the NPAs was obtained by measuring the area (km 2 ) in the GIS environment. NPA shape was calculated from the ratio of the boundary length of the NPAs divided by the area [84]: Values near 1 have a circular shape and are relatively unfragmented and compact. The percent of habitat suitability in 2070 under optimistic (RCP26) and pessimistic (RCP85) scenarios from MRI-CGCM3 and CCSM4 models were used for MCDA analysis. To measure nucleotide diversity (Nd hereafter), sequences were obtained from GeneBank from three genes: NADH dehydrogenase subunit 4 (ND4), control region (D-loop), and NADH dehydrogenase subunit 2 (ND2), with accession numbers MN995079 to MN995069, MK098476 to MK098471, and MK035726 to MK035716, respectively. Distance to nearest PAs (km); villages density; distance to cities (km); distance to the main roads (km) were also measured by GIS. Supplementary Materials Table S3 lists the non-exclusionary criteria as well as a short explanation of why they are important to choose. ArcMap v 10.6.1 was used for all the spatial analyses.

Analytical Hierarchy Process (AHP)
The analytical hierarchy process (AHP), as one of the MCDA methods, was used to obtain weights for non-exclusionary criteria [85]. A matrix cell was used to assess each pair-wise comparison according to the following values: 1: relative to the column variable, the row variable is less important; 3: relative to the column variable, the row variable is moderately more important; 5: relative to the column variable, the row variable is equally more important; 7: relative to the column variable, the row variable is strongly more important; 9: relative to the column variable, the row variable is extremely more important [64].
Finally, all NPAs (polygons) were evaluated using the Suitability Index (SI): The suitability index determines the suitability score for NPAs attained by nth alternative. The Wi is the weight of the factor calculated using the pair-wise comparison between various criteria and Xi is a value obtained from rating curves [64].
In particular, the loss of habitats at low altitudes and the gaining of new habitats at high altitudes were observed. Habitats with 75 to 100% suitability are also more established at higher altitudes or latitude due to the future climate projections ( Figure 5), especially in the CCSM4 model.  Table S4. According to all scenarios and models, the area of habitats with 75-100% suitability will be decreased in 2050 and 2070 (Table S4). Habitats with 50-75% suitability, on the other hand, will expand in size (Table S4).

Siting Procedure
All landscaping features that could damage the N. derjugini were excluded from the study area (Figure 4). These exclusionary criteria including all main roads, all cities, high village density (≥0.72 km 2 ), all dams, the area with very low stream density (≤0.03 km 2 ), areas with poor vegetation, farmlands (concretely irrigated lands), and elevations above 2500 m ( Figure 2). As shown in Figure 4, the results of removing these harmful landscaping criteria from the study area are to create suitable locations where the N. derjugini may be safe. Then, within these areas, with the focus on areas contains N. derjugini habitats, 10 polygons were created as NPAs so that they could accommodate the maximum number of N. derjugini habitats (Figure 4). After that, 13 non-exclusionary criteria were examined, including the number of newt localities, NPA connectivity, NPA shape (edge effect) Figure 5), genetic diversity, distance to nearest PAs, villages density, distance to cities and distance to the main roads. Relationships between the quantities of the non-exclusionary criterion obtained from rating curves as values are shown in Supplementary Materials Table S5. The pair-wise comparison was used to the assigned weight of every non-exclusionary criterion (Supplementary Materials Table S6). The relative importance (weight) of 13 non-exclusionary criteria used to assess the final suitability score of each NPA is shown in Table S6. Finally, the suitability score was ranked for 10 different NPAs based on ratings and weights of the non-exclusionary criteria ( Table 3). Among these NPAs, polygon no. 6 (NPA6: Figure 4) with scores of 0.94 can be the most suitable habitat to introduce as a PA for N. derjugini (Table 3).

Discussion
Protected areas (PAs), as one of the leading conservation tools for biodiversity, are mainly associated with potential failure to protect species against climate and landscape change [19,[86][87][88][89][90][91][92][93][94][95][96]. While some studies have shown that PAs remained suitable climates, facilitated the range expansion and colonization of species, and decreased historical habitat loss [89,[97][98][99][100]. In this study, for the first time, a PA was designed for critically endangered yellow-spotted N. derjugini by integrating eSDM, GIS, and MCDA methods to map and exclude the impact of climate and landscape change. The SDM, by predicting the range of species distribution based on the present condition and the future projection, allows calculating the gain or loss of habitats [101]. The GIS-based MCDA approach allows land management by decision-makers the a landscape scale based on a set of criteria, indicators, and preferences for the area [102]. The methodology proposed in this study can easily be used as an extensible model to introduce a new PA for other animal and plant species and it has already proved to be a powerful tool in natural system development and conservation. This research, on the other hand, can be generalized to other conservation efforts, such as identifying a suitable reintroduction site.

Climate Change
Recent findings highlighted that many amphibians, including salamanders, may face threats as a result of future climatic and landscape change [103][104][105][106][107][108][109][110][111][112][113][114]. Based on modelling results, most of Iran will be exposed to drought within the next 30 years if global warming continues [115][116][117], and this climate change will have an effect on the region's small native amphibians [95]. According to projections by the Intergovernmental Panel on Climate Change (IPCC) for 2100, based on the assumption of continued increase in greenhouse gas emissions in the twenty-first century, Iran may be exposed to temperatures of 1.5 to 4.5 • C under the pessimistic scenarios [118]. In a study by Daneshvar et al. [119], the temperature increase is predicted to range from 1.12 to 7.87 • C by 2100. Studies on the future climatic conditions in Iran show an unpleasant picture, including repeated periods of extreme humidity and drought throughout the country, particularly along the Zagros mountains [115,118]. In the study by Vaghefi, Keykhai, Jahanbakhshi, Sheikholeslami, Ahmadi, Yang and Abbaspour [115], a temperature increase of 1.1 to 2.75 • C has been reported in Iran under four climatic scenarios (RCP2.6, RCP4.5, RCP6.0, RCP8.5), with a temperature increase of 2.25 to 2.5 • C for Zagros under the pessimistic scenario. The analysis of rainfall indicates that the humid regions of the country are becoming wetter, while the arid regions are becoming drier, so rainfall in some areas of the Zagros, especially in the southwestern parts, is expected to increase by 25 to 50 mm [115,120].
The eSDM projections of this study revealed that future climate change will have a greater impact on the edges of species distribution range's ( Figures S1 and S2). According to the MRI-CGCM3 model, in 2050, under both optimistic (RCP26) and pessimistic (RCP85) scenarios, habitats suitable for the N. derjugini will be more lost at low altitudes, but new habitats will be gained at higher altitudes, especially on the southern and southeastern edges of the distribution range ( Figure S1A,C), whereas in 2070, under both scenarios, especially the pessimistic one (RCP85), the edges of the species distribution range will be lost either at low or high altitudes, especially in the northeastern part ( Figure S1B,D). However, under an optimistic (RCP26) scenario, the new habitat would be gained in the southern part of the distribution range at higher altitudes ( Figure S1B). Based on the CCSM4 model, both in 2050 and in 2070, the edges of the distribution area (except for the northwest of the distribution area) also suffered a loss of habitats ( Figure S2). At higher altitudes or latitudes, new habitats were gained in the northwest and particularly in the northeast of the distribution range in this model ( Figure S2). Compared to closely related mountain newts, N. kaiseri, it seems that N. derjugini is more resistant to future climate change. Ashrafzadeh, Naghipour, Haidarian, Kusza and Pilliod [73] found that suitable habitats for N. kaiseri under the influence of future climate change will be reduced by 56% and 96% by 2050 and 2070, respectively. It was also predicted that the species would be shifted to higher altitudes [73].
Low dispersal ability in some species, especially N. derjugini, with an average migration distance of 49.19 ± 71.75 m, see [68], may increase the probability of local adaptation and decrease their distribution range to track favorable climate conditions [121][122][123][124]. The ability of these species to survive in highly variable climates usually depends on their ability to find climate refugia that are relatively buffered from contemporary climate change over time [125,126]. The results of this study are consistent with previous findings of N. derjugini, which suggest that glacial and interglacial cycles affect the species' distribution pattern [127,128]. Climate conditions seem to have been more favorable for the N. derjugini during the Last Glacial Maximum (LGM) in a broader area than today [127]. In contrast, the distribution range of N. derjugini has retreated during the Holocene to current climate conditions [127]. According to the present result, this decline will continue in the future (2050 and 2070). Cold-adapted species usually follow this pattern, expanding their distribution during the last glacial period and contracting it during post-glacial warming [127,[129][130][131].
In a study conducted by Malekoutian, Sharifi and Vaissi [128], three glacial refugia in the southern, central, and northern parts of the N. derjugini distribution range were shown using genetic data. Additionally, in the study conducted by Afroosheh, Rödder, Mikulicek, Akmali, Vaissi, Fleck, Schneider and Sharifi [127], the existence of glacial refugia in the central and southern parts of the distribution range was predicted by SDM. The present study also showed that due to potential future climate change, the range of species distribution would be maintained in the southern, central, and northern parts of the species distribution (except for the margins). Therefore, the results of this study and previous studies [127,128] confirm that the Zagros mountains may act as climatic refugia for N. derjugini. According to recent studies, identifying and subsequently protecting climate refugia, where climates are likely to remain suitable, may be a solution for species conservation [132][133][134].

Landscape Change
In addition to climate change, the negative impact of human activities in and outside PAs, such as agricultural activities, urbanization, roads construction, and industrial and mining activities, as well as negative effects of noise related to these activities on the performance and protection of these areas is well known and has been the subject of much research [135][136][137][138]. In Iran, insufficient management of natural resources, pollution, unchecked urbanization, dam construction, draining of wetlands, deforestation, excessive irrigation, poaching, and lack of scientific and financial support are a severe concern for the loss of biodiversity as a result of landscape changes [139][140][141][142]. These negative changes have been increasing rapidly in recent decades [143,144] and may have an effect on the natural habitat of fauna and flora, animal and plant species, and presentation of ecosystem services [145,146].
The oak forests of the Zagros in western Iran have been used for agricultural activities, grazing, purposes and livestock breeding since about 5 millennia years ago [147,148]. Integration of rapid development and urban expansion with disturbance and traditional livestock grazing as well as severe drought, especially in recent years, are the major issues that have led to changes in the vertical structure or deforestation, configuration, and composition of Zagros mountains forests [149]. In this study, it was not possible to remove some degree of human resources such as some villages, dry farmland, and orchards due to their recognition as enclaves for places with decades or centuries of human inhabitation. However, this study tried to minimize and remove all factors that may damage and decline amphibians inside the introduced PA (Tables S2 and S3). This policy of incorporating PAs in the broader cultural landscape has recently been seen as a vital issue in both conservation and development [150]. However, landscape change must be monitored from the viewpoint of the effectiveness of mixed-use PAs over time [25].

Conclusions
The current study developed a practical approach to integrating eSDM, GISs, and MCDA in order to provide a set of ranked areas by tracking and excluding the effect of climate and landscape changes, which helped managers and other stakeholders in creating PAs. However, more detailed studies on actually selected sites may be needed and can be carried out by field experts. This study can also act as a model for other species that require the establishment of PAs or reintroduction sites.
Supplementary Materials: The following are available at https://www.mdpi.com/article/10.3 390/su13105645/s1; Table S1. The data set used in this study for ensemble species distribution modelling. The information includes the name of localities, geographical coordinates and origin; Table S2. A brief description of the exclusionary criteria used in siting procedure; Table S3. A brief description of the non-exclusionary criteria used in siting procedure; Table S4. Habitat areas (Km2) of the Yellow-spotted mountain newt, Neurergus derjugini, according to different classifications and recent and future (2050 and 2070) periods under two optimistic (RCP 2.6) and pessimistic (RCP 8.5) scenarios within the MRI-CGCM3 and CCSM4 models in western Iran and northeastern Iraq; Table  S5. The no-exclusionary criteria and their relative suitability obtained from 13 rating curves; Table S6. Weight of 13 non-exclusionary criteria result of the pair-wise comparison; Figure S1. Species range change of Neurergus derjugini in currently suitable habitats (gain/loss) by 2050 and 2070 under two optimistic (RCP26) and pessimistic (RCP85) scenarios within the MRI-CGCM3 model in western Iran and northeastern Iraq; Figure S2. Species range change of Neurergus derjugini in currently suitable habitats (gain/loss) by 2050 and 2070 under two optimistic (RCP26) and pessimistic (RCP85) scenarios within the CCSM4 model in western Iran and northeastern Iraq; Scripts S1. The scripts used in the biomod2 settings.  [81]. The Iran main roads, villages and protected areas with a 30 m resolution are available in the Iranian Forest, Rangeland, and Watershed Management Organization [82]. Iraq land cover and the main rod with a 30 s spatial resolution (~1 km) are available in the United States Geological Survey (USGS) and Diva-GIS (http://www.diva-gis.org; accessed on 5 March 2021). The NADH dehydrogenase subunit