**1. Introduction**

The global ecosystem is being affected by human-induced climate change in unprecedented ways, e.g., phenological change, range shifts, community shifts, and so on [1–3]. For wildlife, climate change will produce a series of terrible outcomes that could further contribute to the sixth mass extinction, such as wildlife population declines and extinction, range distribution shifts, habitat fragmentation, and increased evolutionary pressure for all species [4–6]. According to Thomas, et al. [7], a quarter of terrestrial plants and animals may become extinct by 2050, under the 'middle emission scenario' of climate change. Spooner,

**Citation:** Zhuo, Y.; Wang, M.; Zhang, B.; Ruckstuhl, K.E.; Alves da Silva, A.; Yang, W.; Alves, J. Siberian Ibex *Capra sibirica* Respond to Climate Change by Shifting to Higher Latitudes in Eastern Pamir. *Diversity* **2022**, *14*, 750. https:// doi.org/10.3390/d14090750

Academic Editor: Attila D. Sándor

Received: 18 July 2022 Accepted: 8 September 2022 Published: 11 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Pearson and Freeman [4] predicted that mammal populations would decline by 1.46% to 1.76% annually under the scenario of 'representative concentration pathways 8.5'. In Asia, about 30% of terrestrial species will be at very high risk of extinction under climate change [8].

Habitat shifts and eventual loss caused by climate change are the most pervasive threats to wildlife causing species to migrate towards higher latitudes and/or higher elevations [9–11] and fragmenting suitable habitats [12,13]. For example, Hickling, et al. [14] predicted that 275 of 329 terrestrial species in Europe would move north in response to climate change. Alpine and plateau species are particularly sensitive to climate change as they are typically cold-adapted or cold-tolerant species [15]. In fact, studies have demonstrated that climate change has caused ibex (*Capra ibex*), chamois (*Rupicapra rupicapra*), and red deer (*Cervus elaphus*) in the Alps to move to significantly higher latitudes (i.e., to the north) [16]. Similarly, due to climate change, four antelope species on the Tibetan Plateau (*Procapra przewalskii*, *Gazella subgutturosa*, *P.picticaudata* and *Pantholops hodgsonii*) are predicted to lose up to 53.2% of their habitat and will be forced to move to higher latitudes [17]. In addition, in 90 years' time, the suitable summer range for Alpine ibex will diminish to 26.4% compared with 2011 under the RCP8.5 scenarios, and they will respond to high temperatures by moving to higher altitudes [18], as in the case of chamois who move upslope when it is hotter [19].

The alpine zone of the Pamirs is more sensitive to climate change because the temperature rise is much greater than that in the low-altitude areas [20,21]. Previous studies on ungulates living in the Pamir plateau have shown that climate change will reduce the range of these animals [22–24]. For example, the suitable habitat of Marco Polo sheep (*Ovis ammon polii*) in Eastern Tajikistan may be reduced by more than 60.6% and 63.2% by 2050 and 2070, respectively [22]. Because of climate change, the Chinese population of the same species is losing close to 40% of its suitable habitat, and this is mostly happening at low elevations [24].

The Siberian ibex, *Capra sibirica*, is a typical mountain ungulate, which inhabits the mountains of Southern Siberia, Mongolia, Tianshan, Himalaya, Karakorum, Altai, and Pamir [25,26]. Despite its wide distribution, the species remains poorly studied [27]. Siberian ibex are listed as "near threatened species" by the International Union for Conservation of Nature (IUCN) [25], and as "state second class protection animals" by China [28]. They usually occur at higher elevations from around 1500 m up to 6000 m above sea level. [29].These ibex mostly inhabit mountains, alpine meadows, exposed rock, precipitous terrain interspersed with cliffs, deep valleys, and steep areas with escape terrain [30]. The Siberian ibex mainly feeds on herbs and shrubs in the middle part of the Tianshan Mountains, which is one of seven of the largest mountain systems in the world, located in remote areas of Eurasia [31]. Throughout their distributional range, Siberian ibex show seasonal movement to higher ridgelines in the summer, descending to lower elevations in the winter [32,33]. Such movement can be entirely vertical or reach travel distances of 40 to 100 km [22]. However, Siberian ibex face disturbance from human activities and domestic competition (i.e., hunting and grazing), and climate change, with changes in temperature, wind, and precipitation listed as the main drivers affecting the survival of this species [25,29]. Therefore, it is urgent to study the potential responses of the Siberian ibex to climate change across its distribution range.

In this study, we analyzed the current distribution of Siberian ibex in the Taxkorgan Nature Reserve (TNR), which is located in the southwest of Xinjiang, China, and the factors affecting their distribution by using ensemble species distribution models (eSDMs), which are considered to improve predictive performance compared to individual species models [34,35]. We focused our efforts on determining whether or how the distribution of Siberian ibex will change in response to future global warming. We predict that the suitable habitat for Siberian ibex would decrease in the TNR in the future, and they would be forced to migrate to higher latitudes and/or elevations as a response to climate change. Our study will contribute to understanding the future of ungulates on the Chinese Pamir

plateau under climate change. Furthermore, the study provides scientific guidance for the TNR's conservation planning, which needs to take into account potential future changes in the distribution of this species when designing conservation policies and planning the functional zonation of this nature reserve.

### **2. Methods and Materials**

### *2.1. Study Area*

The study area is situated in the Taxkorgan Nature Reserve (TNR) (35.5◦ to 37.5◦ N, 74.4 to 77.1◦ E, Figure 1, an area of ~16,253 km2) in the eastern part of the Pamir, which is mainly distributed in Taxkorgan County, Xinjiang, China. As part of the Qinghai-Tibetan Plateau, the elevation of the TNR ranges from 2098 m to 8572 m above sea level, with an average elevation above 4000 m. The area has limited precipitation and a dry climate [36]. The average annual precipitation is less than 100 mm, the mean annual temperature is 3 ◦C, the highest temperature in July reaches 38 ◦C, and the mean temperature of the coldest quarter is −17 ◦C. The region is rich in wildlife. Besides Siberian ibex, sympatric ungulates include Marco Polo sheep and blue sheep (*Pseudois nayaur*, AKA bharal). All three species are preyed upon by the snow leopard (*Panthera uncia*). The vegetation is dominated by dwarf shrubs (mainly *Artemisia L.* and *Ceratoides*) and grasses (mainly *Stipa*) [37].

**Figure 1.** The study area is located in the west of China, which borders Tajikistan, Afghanistan, and Pakistan. (**a**) A map of China. It reflects the location of Xinjiang, a province in China. (**b**) Taxkorgan Nature Reserve (TNR) is located in the southwest of Xinjiang, China. (**c**) The map shows the borders of the Taxkorgan Nature Reserve, with elevation indicated on a green to orange scale. The Karakoram highway is indicated as a red line transect.

Most human settlements are located below 3500 m above sea level, i.e., along both sides of the Karakoram highway, which runs south from the city of Kashi up to the Taxkorgan Valley, through the Marco Polo sheep distribution area, and over the Khunjerab Pass into Pakistan. The southern part of this route lies in the TNR. The reserve is sparsely populated by nomadic Tajik and Kirgiz herding people and their livestock (about 23,000 in total), who live on both sides of the Karakoram highway and in the main valleys inside of the reserve [38].They move yearly between winter (lower elevation areas) and summer pastures (higher elevation areas with productive vegetation) [39,40]. In some areas of the TNR, Siberian ibex share their habitat with domestic animals, which compete with them for food and space. Since 2018, all of the coal and jade mines have been moved out of the reserve, so the herding people and their livestock are the main human disturbance to wild animals in the TNR.

### *2.2. Occurrence Data*

We conducted a survey of Siberian ibex in the TNR from June to November of 2018 to 2021 (for more details on the line transect survey time please see Supplementary Materials Table S1). After interviewing the reserve staff and herding people, we set up 23 line transects with a total length of 1332.74 km, which covers all the habitats occupied by Siberian ibex in the TNR. Each transect was visited once a year. We followed the transects on foot, by car, and horseback according to the terrain, and scanned the surroundings with binoculars (ZEISS 10 × 42, Oberkochen, Germany) at intervals of 2–3 km to spot Siberian ibex. Upon locating Siberian ibex, we used a telescope (ZEISS 20–60×, Oberkochen, Germany) and laser range finder (ZEISS T\*RF 10 × 54, Oberkochen, Germany) to observe and record their location (longitude and latitude), their elevation, and their angle and distance from our observation point. A total of 495 locations were set up along the line transects to search for Siberian ibex, and 123 occurrences of ibex were recorded, i.e., 26, 29, 36 and 32 distribution locations were obtained between 2018 and 2021. To reduce any potential spatial sampling bias (only one occurrence point is retained in each 1 km × 1 km grid), we performed spatial thinning of Siberian ibex occurrence data. The spatial thinning was performed using a randomization approach that maximizes the amount of useful information retained [41–43], while simultaneously reducing the sampling bias by removing as few records as possible using the 'thin' function of the "spThin" package (version 0.2.0) [44] in R (version 4.1.0) [45] with a thinning distance of 1 km. The spatial thinning process kept a total of 109 occurrence data points that were used to build the species distribution models.

### *2.3. Environmental Variables*

We considered different environmental factors that can influence the presence of Siberian ibex, including climate, topography, anthropogenic disturbance, and availability of resources (i.e., Normalized Difference Vegetation Index (NVDI)) (Table 1). Climate data was obtained from the WorldClim version2.1 dataset for 1970–2000 (www.worldclim. org (accessed on 20 September 2020)) [46], with a resolution of around 30 arcseconds (roughly 850 m at the study area latitude). To predict the influence of climate change on Siberian ibex, we selected two projection periods for our model: (1) 2050 (average for 2041–2060) and (2) 2070 (average for 2061–2080) and three greenhouse gas emission scenarios: (1) RCP2.6 (low greenhouse gas emissions [47]), (2) RCP4.5 (medium greenhouse gas emissions [48]), and (3) RCP8.5 (high greenhouse gas emissions [49]) based on the Representative Concentration Pathways (RCPs) of Coupled Model Intercomparison Project Phase 5 (CMIP5). We chose the BCC-CSM1-1 as the General Circulation Model (GCM), which more closely conforms to the climatic conditions of our study area [50]. The 30 m resolution elevation data was derived from the Geospatial Data Cloud site of the Chinese Academy of Sciences (www.gscloud.cn (accessed on 2 September 2021)), from which we also calculated the aspect (the range of 0–360), slope (the range of 0–85), and ruggedness of the area, using ArcGIS 10.6 (ESRI). The latest version (1995–2004) of Human Influence Index, covering population density, human infrastructure, and human access was obtained from the Socioeconomic Data and Applications Center (sedac.ciesin.columbia.edu (accessed on 25 August 2021)) [51]. We derived the global land cover data from GlobeLand30 (www.globeland30.org (accessed on 3 July 2022)) and the NDVI from 2018 to 2021 from the DATABANK Remote Sensing Data Engine (databank.casearth.cn (accessed on 3 July 2022)) [52]. NDVI and HII data used in the future scenarios' models remain unchanged from the current dataset. All factor layers were projected to WGS 1984 UTM Zone 43N and resampled to a resolution of 850 m. To avoid collinearity between the environmental variables, we used the "usdm" package (version 1.1-18) [53] in R (version 4.1.0) [45]. We calculated the variance inflation factor (VIF) between predictor variables for each scenario, and then used a stepwise procedure to eliminate the factors with VIF > 10 [54]. All the

variables that entered the ensemble species distribution models of the current and future scenarios were selected by the above processes when modeling the ibex habitat were labeled with an asterisk (Table 1).

**Table 1.** The environmental variables used in the ensemble species distribution models (eSDMs) of Siberian ibex in Taxkorgan Nature Reserve.


Note: The asterisk (\*) means this variable entered the ensemble species distribution models. The future environment variables are consistent with the current environment variables, both of which have been tested for collinearity, and they are all labeled variables in this table \*. \*\* For more detailed information about climate variables see www.worldclim.org (accessed on 20 September 2020).

### *2.4. Species Distribution Model*

The ensemble species distribution models (eSDMs) combine predictions across different modeling methods [55]. There have been studies showing that the performance of eSDMs is better than single models [34,35]. We used eight species distribution algorithms (generalized linear models, generalized boosting models, generalized additive models, artificial neural networks, flexible discriminant analysis, multivariate adaptive regression splines, random forest, and MaxEnt) [55] to build the eSDMs for Siberian ibex in the TNR. Firstly, we used a random strategy [56] to create 5000 pseudo-absence points [57]. We then split our occurrence data into two parts: one part with 80% of the data was used to calibrate the models, and the remaining part (20%) was used for model testing. All algorithms were run ten times and to obtain eighty single models. Then, we took the average TSS value of these eighty models as the threshold of whether to enter the eSDMs (above average TSS were entered in the final model) and built the ensemble model with a weighted average to project the habitat suitability map for Siberian ibex. The continuous probability of the presence of the eSDMs was transformed into binary values using a cut-off

value that maximized TSS to estimate species vulnerability [58,59]. The AUC value (the area under the receiver operating characteristic curve; a value greater than 0.9 is considered excellent) and TSS value (a value greater than 0.85 is considered excellent; a value between 0.7 and 0.85 is good) were used to evaluate the performance of the eSDMs [60,61]. All eSDM analyses were performed using the "biomod2" package (version 3.5.1) [62] in R, version 4.1.0 [45].

### **3. Results**

### *3.1. Current Habitat Distribution*

The eSDMs for Siberian ibex in the TNR showed excellent performance with a high AUC (0.94) and TSS (0.77) value. Among the potential factors affecting the current distribution of Siberian ibex, precipitation seasonality, precipitation of the wettest month, Human Influence Index (HII), NDVI, elevation, and temperature seasonality were the main determining factors (Supplementary Materials Table S2). The presence probability of Siberian ibex was positively correlated with NDVI, and elevations of 3000 to 4800 m (Supplementary Materials Figure S1).

The current suitable habitat for Siberian ibex in the TNR is mainly distributed in the northwest and a small part of the northeast (Figure 2). The suitable habitat area is around 2702.15 km2, accounting for 16.6% of the total area of the natural reserve.

**Figure 2.** The current habitat suitability of Siberian ibex in Taxkorgan Nature Reserve.
