Study area
Mazandaran province (Fig. 1), located in the north of Iran, is surrounded by the Caspian Sea from the north and the Alborz Mountain range from the south. With an area of about \(\text{23,000}\hspace{0.17em}k{m}^{2}\) and elevations ranging from 20 meters below sea level to over 5000 meters. This region is home to numerous swamps and sea-level plains, vast rice fields, and mountains covered with dense forests and rivers. According to the national census conducted in 2016, three million people live in Mazandaran, of which three hundred thousand are engaged in agriculture, animal husbandry, and fishery (Iran 2019).
Studied species and species data
We utilized field observations of species between 2015 and 2018, conducted by the Department of Environment of Mazandaran, Iran. According to the field report, surveys were based on direct observation and sound recognition. We used grid sampling in order to resample data in a one-kilometer radius. Finally, 631 occurrence records were prepared and used in distribution modeling.
Current and future land cover and climate data
This research used nine variables to model the common pheasant geographic distribution, including four bioclimatic and land-cover parameters and a topographic parameter. We obtained land-cover data for the year 2020 from the ESA-CCI (European Space Agency - Climate Change Initiative) project and generalized that into five classes of forest, cropland, pasture, built-up areas, and others based on IPCC classes for change detection (ESA 2017). We used all but the last class (“others”) in modeling, as they are more relevant to this species distribution and its nesting and breeding (Coates et al. 2017; Perkins et al. 1997; Riley et al. 1998; Smith et al. 1999). Correlation analysis of parameters is given in Supplementary Table S1.
Three scenarios of RCP2.6, RCP4.5, and RCP8.5 were incorporated to project the possible climate variations in 2040–2060. The aforementioned “Representative Concentration Pathways” (RCPs) are global radiative forcing levels of 2.6, 4.5, and 8.5 \(W/{m}^{2}\) by the end of the century, and describe climatic scenarios leading to a very low forcing level (RCP2.6), a medium stabilization scenario (RCP4.5), and a very high baseline emission scenario (RCP8.5) (Van Vuuren et al. 2011). We used averages of multiple GCPs (Global Climate Models), namely CCSM4, CNRM-CM5, HadGEM2-ES, MIROC-ESM-CHEM, and MPI-ESM-LR as bioclimatic variables for each scenario.
We used the SSP scenarios to estimate how the land cover might change in the study area. These scenarios are based on five distinct narratives. Each scenario describes the earth's socio-economic trends and environmental changes in the twenty-first century. These scenarios aim to identify the socio-economic challenges for mitigation and adaptation purposes (Riahi et al. 2017). SSP1 and SSP3 describe scenarios, one with low and another with high challenges for both adaptation and mitigation. In addition, SSP4 is characterized to represent low adaptation and high mitigation narrative, and SSP5 vice versa. Ultimately, SSP2 is designed to portray a world with intermediary challenges for both adaptation and mitigation (Popp et al. 2017).
Six IAMs (Integrated Assessment Models) have been developed to quantify these scenarios based on different assumptions and uncertainties. Each SSP has a baseline scenario containing only current climate change policies to estimate future conditions. In addition to these baseline scenarios, other scenarios have been developed that use different climate change policies. We used baseline scenarios of the IMAGE model to project future land cover change for the Middle East and Africa (MAF) (Riahi et al. 2017; Van Vuuren et al. 2017). SSP scenario data is available at https://tntcat.iiasa.ac.at/SspDb/ (accessed on 2022/04/05). Alongside SSP scenarios, we calculated the land cover change between 2015 and 2020 and used linear regression to estimate the total change in 2050. We named this the Trend scenario.
We calculated the change ratio of the cropland, forest, pasture, and built-up land cover classes between 2020 and 2050. We then multiplied this by the area of these classes in our study area to project land cover changes for 2050. The summary of individual land cover changes is presented in Table 1.
Table 1
Land-cover changes between 2020 and 2050. Values are occupied areas for each land type (km2) and change concerning the Baseline
| Baseline | Trend | SSP1 | SSP2 | SSP3 | SSP4 | SSP5 |
Cropland | 8230 | 9019 (10%) | 9872 (20%) | 10575 (28%) | 11003 (34%) | 10997 (34%) | 11130 (35%) |
Forest | 7604 | 6520 (-14%) | 6594 (-13%) | 5310 (-30%) | 4657 (-39%) | 4813 (-37%) | 5267 (-31%) |
Pasture | 4475 | 4351 (-3%) | 4364 (-2%) | 4884 (9%) | 5067 (13%) | 4936 (10%) | 4994 (12%) |
Built-up | 392 | 1390 (255%) | 766 (95%) | 739 (88%) | 714 (82%) | 857 (119%) | 744 (90%) |
Other | 2856 | 2278 (-20%) | 1962 (-31%) | 2050 (-28%) | 2117 (-26%) | 1955 (-32%) | 1424 (-50%) |
Following Préau et al. (2019), we used Dyna-CLUE to spatially allocate land changes in the area. This program comprises two parts: the spatial allocation section and non-spatial aggregated demand. The spatial allocation includes land conversion matrix, conversion elasticity, and location characteristics. The conversion matrix represents which lands can be converted to each other. The conversion elasticities indicate the expense and difficulty of converting one land to another. We used logistic regression to determine the potential characteristics of each land type regarding environmental variables (Verburg and Overmars 2009). Logistic regression was deployed using the Sickit-Learn library in Python programming language. Eventually, the land-cover map for 2050 was simulated by entering the spatial parameters and land-cover change scenarios into the Dyna-CLUE (For detailed information, check Supplementary Table S3-S9).
Species Distribution Modeling
We used Maxent 3.4.1 (Phillips et al. 2017), a presence-only algorithm that estimates a species' geographic distribution by determining the distribution with the highest entropy given a set of features derived from environmental conditions. Maxent implements linear (L), quadratic (Q), product (P), threshold (T), and hinge (H) of provided environmental variables as features. Linear features are equal to continuous environment variables, quadratics are their squares, and products are products of pairs of continuous variables. Threshold makes a binary feature based on an arbitrary value. The hinge feature is similar to the linear feature, but it is constant below a certain threshold (Phillips et al. 2006).
Additionally, Maxent uses the LASSO penalty to reduce model overfitting, which can be adjusted by a multiplier constant (Merow et al. 2013). Parameter optimization in Maxent involves feature selection and regularization tuning (Li et al. 2020). Maxent was executed based on a set of predetermined parameters, and the parameters that reached the best AUC were selected. This procedure was accomplished using the ENMeval library in R programming language (Kass et al. 2021).
To reduce bias in modeling, we masked one county in which there was no conducted survey and implemented a kernel-density map of all occurrences as a bias file (Fourcade et al. 2014). We ran Maxent with the maximum iteration of 500, and data partitioning was based on five-fold cross-validation. Finally, we converted continuous habitat suitability to binary outputs to compare scenarios and the amount of lost and gained areas. We used the maximum test sensitivity plus specificity threshold for this task (Liu et al. 2016).