Habitat assessment of Marco Polo sheep (Ovis ammon polii) in Eastern Tajikistan: Modeling the effects of climate change

Abstract Identifying the factors predicting the high‐elevation suitable habitats of Central Asian argali wild sheep and how these suitable habitats are affected by the changing climate regimes could help address conservation and management efforts and identify future critical habitat for the species in eastern Tajikistan. This study used environmental niche models (ENMs) to map and compare potential present and future distributions of suitable environmental conditions for Marco Polo argali. Argali occurrence points were collected during field surveys conducted from 2009 to 2016. Our models showed that terrain ruggedness and annual mean temperature had strong correlations on argali distribution. We then used two greenhouse gas concentration trajectories (RCP 4.5 and RCP 8.5) for two future time periods (2050 and 2070) to model the impacts of climate change on Marco Polo argali habitat. Results indicated a decline of suitable habitat with majority of losses observed at lower elevations (3,300–4,300 m). Models that considered all variables (climatic and nonclimatic) predicted losses of present suitable areas of 60.6% (6,928 km2) and 63.2% (7,219 km2) by 2050 and 2070, respectively. Results also showed averaged habitat gains of 46.2% (6,106 km2) at much higher elevations (4,500–6,900 m) and that elevational shifts of habitat use could occur in the future. Our results could provide information for conservation planning for this near threatened species in the region.

II and as Near Threatened in the IUCN Red List. Argali have declined in numbers and distribution during the last century (Harris & Reading, 2008;Valdez & Weinberg, 2011). Eastern Tajikistan has greater numbers of argali than any other country with a minimum of 24,000 in the Pamir region (Michel & Muratov, 2010;Valdez, Michel, Subbotin, & Klich, 2016). Due to the high elevation of the Pamirs, the region is extremely susceptible to climate-driven threats, such as glacier recession (Khromova, Osipova, Tsvetkov, Dyurgerov, & Barry, 2006), and decrease in water storage and supply (Finaev, Liu, Bao, & Li, 2016) that could harmfully impact wildlife habitats (BIOFOR 2001). Glaciers play a crucial role in the hydrological cycle of highaltitude regions (Nogués-Bravo, Araujo, Errea, & Martinez-Rica, 2007) as they, together with snow packs, provide freshwater and soil moisture necessary for the survival of vegetation communities.
Because argali sheep sightings occur nearest to water sources where available forage is abundant (Salas, Boykin, & Valdez, 2016), and riparian habitats have been shown to be the strongest predictors for argali occurrence in Tibet (Singh, Yoccoz, Bhatnagar, & Fox, 2009) and in the Pamirs of Tajikistan (Salas, Valdez, & Michel, 2017), increased temperatures in mountainous regions could decrease forage availability in proximity to riparian areas and hence the amount of suitable habitat for argali. IPCC (Intergovernmental Panel on Climate Change) (2007) reported that retreat of various large glaciers could accelerate through the 21st century, which would initially increase summer flow but, in the long term, could reduce water availability in regions supplied by meltwater from major mountain ranges in Central Asia. This increased water scarcity could cause extensive habitat loss and degradation of mountain ecosystems (Breu & Hurni, 2003). Hole et al. (2009) highlighted the need to project future habitat for species affected by climate change in order to aid in conservation and management efforts and identify critical habitats. Environmental niche models (ENMs) are valuable spatial ecological tools to better assess the relationship between species distributions and environmental factors, and understand future steps for species management and policy (Elith & Leathwick, 2009;Salas, Seamster, Boykin, Harings, & Dixon, 2017;Wiens, Seavy, & Jongsomjit, 2011). In this study, we used niche models (Elith & Leathwick, 2009;Peterson et al., 2011) to project availability of suitable environmental conditions for argali species using various climate projections derived from general circulation models (GCMs), and postprocessed via application of a simple statistical downscaling method. We contrasted future projected climate envelope suitability results produced from combinations of four GCMs and two greenhouse gas concentration trajectories for two future time periods. Our methodology used five niche models to predict environmental argali habitat suitability in contrast to all other wild sheep habitat climate change studies which used a single modeling technique. To account for the uncertainties of each statistical model and improve predictions of the current distribution of a species (Araújo, Whittaker, Ladle, & Erhard, 2005;Marmion, Parviainen, Luoto, Heikkinen, & Thuiller, 2009), we did an ensemble of five ENMs to assess the agreement or disagreement of their predictions. Combining ENMs within an ensemble could add information not shown in one single algorithm (Meller et al., 2014). This is the first modeling study in Tajikistan to assess the impacts of climate change on argali habitat and the first quantitative study of factors driving climatic distribution patterns of argali. This study is directed toward filling this knowledge gap by combining statistical models with field sightings of argali collected in eastern Tajikistan.
Our objective is threefold: to develop models of present-day and potential future distributions of suitable environmental conditions for argali in mountainous regions in eastern Tajikistan, to compare how climate change impacts the future availability of habitat for argali, and to provide information of the diverse climatic and nonclimatic variables that could potentially affect argali habitat. Lastly, we hypothesized that there would be elevational habitat shifting of argali under future climate conditions when projected to years 2050 and 2070.

| MATERIAL S AND ME THODS
The methodology was composed of three major steps: (i) processing of argali occurrence dataset and environmental variables, (ii) modeling of current environmental conditions that includes selection of ENMs and evaluation of current conditions, and (iii) modeling of future environmental conditions that includes selection of GCMs, representative concentration pathways (RCPs), and projection of current conditions to future conditions. We modeled argali habitat with two approaches, first using a model consisting of climatic and nonclimatic variables (all environmental variables), and a second model with only climatic variables.

| Study area
The study area ( Figure 1) is located in the eastern region of Tajikistan between latitudes 37°N to 40°N and longitudes 72°E to 76°E in the province of Gorno-Badakhshan, and comprises approximately 41,000 km 2 . About 93% of Tajikistan is mountainous, and over half of the country is situated above 3,000 m. The terrains in the east of the country form the highest mountain systems in Central Asia, the Pamirs. The Pamir region is known as the "roof of the world" with wide and grassy valley floors, and meandering rivers and streams (Breu & Hurni, 2003). In the study area, elevations in the northern region range from 4,500 to 7,300 m, while those in the southern region range from 3,300 to 5,200 m above mean sea level (Salas, Valdez, & Boykin, 2015). The January mean temperature in the eastern Pamirs varies from −15 to −20°C with extreme seasonal temperature variations (WHC 2013). At lower elevations, the average January temperature is −1 to 3°C (NCCA 2003). Because of the specific climate conditions and varied landscape in Tajikistan, the country is deemed the main glacial center of Central Asia. The Fedchenko Glacier (907 km 2 ), the largest glacier in Central Asia, is situated in the northern Pamir region of Tajikistan (Merzlyakova, 2002). In fact, the Tajik Pamirs alone provide approximately 60% of the freshwater reserves of Central Asian lowlands (Badenkov, 1992).
The highest elevations of the northern region of the study area are in the Tajik National Park. The park is one of the largest high mountain protected areas in the Palearctic realm (Haslinger, Breu, Hurni, & Maselli, 2007). The state agency of Natural Protected Areas supervises all management activities in the park (WHC 2013).
However, even in the Core Zone activities such as fuel wood cutting, livestock grazing, and hunting tourism take place (Weaver 2013), including a sport hunting concession primarily for Marco Polo sheep .
The only other wild ungulate species in the study area is the Siberian ibex (Capra sibirica). Wild predators include wolf (Canis lupus), red fox (Vulpes vulpes), brown bear (Ursus arctos isabellinus), and snow leopard (Panthera uncia). There is limited mining activity at the southern edges of the park, and only one paved road with little traffic, and some unpaved roads in the study area. A barbed wire fence forms a barrier to argali movements with few gaps along the Tajik-Chinese border. Most of the area becomes inaccessible in the winter because of heavy snow accumulation. The majority of the study area outside of the park is assigned to private businesses as hunting concessions. Some of the hunting concession areas are patrolled to minimize illegal hunting, while other areas are unprotected . Domestic ungulates include sheep (Ovis aries), goat (Capra hircus), yak (Bos grunniens), and cattle (Bos Taurus); domestic sheep are the most numerous followed by yak and few cattle and goats. Domestic animals, except for yak and few herds of sheep and goats, are moved to lower pastures during the fall, winter, and early spring (October-May) because of the harsh winter conditions at higher elevations.
F I G U R E 1 Location of study area in the eastern region of Tajikistan (inset) between latitudes 37°N to 40°N and longitudes 72°E to 76°E, and covers an area of approximately 41,000 km 2 . The Central Asian country of Tajikistan is bordered by Afghanistan, China, Kyrgyzstan, and Uzbekistan

| Argali data
Argali occurrences were based on observations from field surveys over multiple years. Using the GPS location of the observers, the distance (measured by range finder or roughly guessed) and azimuth (measured with electromagnetic compass) between observers and observed animals, sightings of argali herds were georeferenced (Michel & Muratov, 2010;Valdez et al., 2016). Field surveys were conducted in March 2015, July 2015, August 2010, 2016, September 2013, and December 2009. Because clustered point occurrences could introduce potential spatial bias, our algorithm removed multiple presence localities in a grid and analyzed a single occurrence per pixel (Salas, Seamster, et al., 2017). A total 917 of 976 locations were used in the model. Further, we divided all occurrences into 70%-30% training-testing subsets when running our models.

| Environmental variables
To build the models, we used 26 environmental variables (Table 1), of which 19 were bioclimatic variables from WorldClim datasets (Hijmans, Cameron, Parra, Jones, & Jarvis, 2005), and the rest were habitat and topographic variables that were previously identified as preferable for argali habitat (Salas, Valdez, et al., 2017).
WorldClim provides climate projections statistically downscaled using a "delta method" approach to a spatial resolution of 30 arcsec, roughly 900 m at the equator. All raster images were resampled to the spatial resolution of the WorldClim data. We derived the green vegetation cover and sparse vegetation cover from the Global Land Cover 2000 project (Tateishi, Zhu, & Sato, 2000) and used them to represent forage abundance. Previous studies considered vegetation distribution as an important predictor in modeling the argali habitat (Salas et al., 2015;Singh et al., 2009). For escape terrain, we created a continuous distance around polygon patches with slopes ≥ 30° (Smith, Flinders, & Winn, 1991;Turner, Douglas, Hallum, Krausman, & Ramey, 2004). The processed DEM with a 1 arc-second, or about 30 m (98 feet) resolution was sourced from NASA's Shuttle Radar Topography Mission (SRTM) digital elevation dataset that is available for download online (USGS 2016).
We computed the slope and the aspect from the DEM. To further capture the relief characteristics of the landscape terrain, the terrain ruggedness index (TRI) was calculated from the DEM based on Sappington, Longshore, and Thomson (2007). Turner et al. (2004) showed that TRI could be a better predictor than proximity to escape terrain when both are used in the same modeling set. The TRI shows the average change in elevation between a center pixel and its eight neighboring pixels in a 3-by-3 window. All candidate variables were clipped to the extent of the study area.

| Species distribution modeling
To model the potential distribution of argali, we used the following statistical methods: generalized linear model (GLM), Random Forest (RF) (Breiman, 2001;Liaw & Wiener, 2002), boosted regression tree (BRT) (Elith, Leathwick, & Hastie, 2008) (Leathwick, Elith, & Hastie, 2006). We selected ENMs based on their performance with presence-only data . GLM, MARS, and BRT could be used for count data under the assumption that a count response could be modeled as Poisson (Talbert, 2012). The GLM is a linear regression adapted to binary count data. The method uses stepwise procedure to select covariates in the model. The MARS nonparametric algorithm builds flexible models by fitting piecewise logistic regressions. Although it has similarities with GLM, MARS is better in accommodating nonlinear responses to predictors and at the same time lessens the effects of outlying observations. The model RF uses decision trees through random grouping of the covariates. Random Forest models both interactions of the variables and their nonlinear relationships, and does not split the data into training and testing as RF utilizes bootstrapping to fit individual trees (Breiman, 2001). Like the Random Forest, BRT also uses decision trees, but the method is robust to missing observations. BRT uses cross-validation by choosing models based on model comparisons of evaluation metrics (Elith et al., 2008). Maxent is best for presence-only modeling. While observed absence is valuable in modeling, data are often not available and hence using only presence data are unavoidable (Talbert, 2012). Models were implemented in the modeling tool Software for Assisted Habitat Modeling (SAHM) run within VisTrails (Morisette et al., 2013;Talbert, 2012).
The tool creates a workflow of the selected ENMs and develops models for present-day conditions. As species lacked absence data, the tool randomly generated background points (i.e., pseudoabsences [Phillips & Dudík, 2008]) within a 95% minimum convex polygon defined by the presence data.
We made the choices between variables based on the results of a species-specific literature search. In particular, we selected variables that were identified in one or more studies regarding the argali as having an effect on the argali's range or population dynamics. In cases where the results of the literature search could not differentiate between two highly correlated climatic variables, we used a qualitative assessment of the distribution of values of the variable at all presence points and of the relationship between the variable and species presence or pseudo-absence (Talbert, 2012). We provided an asterisk in Table 1 to denote the final variables that were selected by our algorithms when modeling the argali habitat.
We produced ensemble maps for the current distributions. The ensemble map is a summation of binary maps generated from probability surfaces from each statistical modeling algorithm (Liu, Berry, Dawson, & Pearson, 2005;Lobo, Jiménez-Valverde, & Real, 2008;Stohlgren et al., 2010). We optimized the threshold using specificity = sensitivity in discretizing the probability maps (Manel, Williams, & Ormerod, 2001). The final maps consisted of pixel values that represented the number of models in agreement to indicate that a particular pixel is suitable for argali. A pixel with a value of zero meant that none of the models identified bioclimatic suitability for the species at that location, while a value of 5 meant there was agreement across all five models. We used the current distributions estimated by the ensemble ENMs and projected each to the future.

| Model evaluation
We evaluated confidence in individual model results in terms of concordance among the different distribution models. We had higher confidence that environmental conditions were suitable for a species when three or more (at least 60% of) algorithms were in agreement (e.g., Rehfeldt, Crookston, Sáenz-Romero, & Campbell, 2012

| Model performance and variable importance
The AUC values for all models were above 0.80 (Table 2). Statistics for models with climate variables only are enclosed in parenthesis.

| Current distribution models
Using all variables, the five ENMs showed similar patterns of potential suitable conditions for argali in the northern and southeastern regions of the study area ( Figure 2). Visual assessment of the maps revealed some differences, however. The five statistical models projected differing size of suitable areas even though they were using the same dataset of occurrence points. Among the five models, only RF and BRT (Figure 2a (Sohl, 2014).

| D ISCUSS I ON
The effects of climate change are already being detected in the shifting and contraction of species' ranges (Raxworthy, 2008;Sexton et al., 2009;Su, Aryal, Nan, & Ji, 2015), variations in species' ecological interactions (IPCC 2007), and suboptimal movement behavior patterns (van Beest & Milner, 2013). This change could be deleterious (Lannoo, 2005), leading to reduction in species abundance as habitats deteriorate (Halpin, 1997). When global temperatures exceed 2.5°C, there could be negative con-

| Present and future projections
This is the first study to model climate change and map the current and future habitat suitability of argali in Tajikistan (and in Asia). This is also the first study to use an ensemble of modeling algorithms and climate models for a large mammal species endemic to the TA B L E 3 Areas (km 2 ) and percentages of potential habitat stability, gain, and loss of Marco Polo argali for projections to 2050 and 2070 under climate change. Suitable habitat is considered stable when present and future ensemble maps agree that the area is suitable for the species. Gain in suitable habitat is recorded when projections find suitable habitat beyond the extent of the present suitable areas. There is loss of suitable habitat when present suitable areas become unsuitable in the future Stable (km 2 ) (%) Gain (km 2 ) (%) Lost (km 2 ) (%) Stable (km 2 ) (%) Gain (km 2 ) (%) Lost (km 2 ) (%) mountainous regions of Central Asia. We modeled the argali suitable habitats under current and future climates using two datasets, one with only climatic variables and another set that included other environmental variables. Models of present conditions were controlled by annual mean temperature, suggesting that this important climatic variable was a major limiting factor for argali. The association of argali to low temperatures was also reported in Khan et al. (2015), where it was labeled a major factor influencing species distribution and habitat suitability for ungulates in the Karakoram-Pamir Landscape between China and Pakistan. Other studies also acknowledged the importance of temperature in modeling ungulates habitat suitability (Aryal, Raubenheimer, & Brunton, 2013;Forrest et al., 2012), although some used elevation as a substitute for temperature (Aryal, Brunton, Ji, et al., 2014). Apart from the annual mean temperature, we expected the precipitation data in the warmest quarter to rank high in the models, which our results confirmed.
Other modeling studies for argali such Khan et al. (2015) did not include precipitation in the warmest quarter as a predictor. Abundant precipitation in dry months is key to moist meadows and forage development. St-Louis and Côté (2014) found that forage quality can be a key factor determining habitat selection patterns for large herbivores. Likewise, the vegetation cover was found by our models to be an important contributor to the predicted habitat of Marco Polo argali. Salas et al. (2015) showed that the availability of forage and the distribution of argali are strongly correlated. Previous studies also reported cases where the normalized difference vegetation index (NDVI) (Guyot & Gu, 1994), used as an index of forage abundance (NDVI > 0.4) (Singh et al., 2009), was linked to the distribution of large herbivores (Pettorelli, 2014). The slope and elevation did not rank high on importance unlike previous studies showing them as major limiting factors for ungulates (Aryal, Brunton, Ji, et al., 2014;Forrest et al., 2012). However, terrain roughness was a significant predictor of habitat probably because the variable is, in part, defined by slope. Suitable areas for argali were found primarily on gentler slopes (0° to 15°), as also observed by Singh et al. (2009), Chetri andPokharel (2005), and Namgali, Fox, and Bhatnagar (2004).   Radi et al., 2014). Also, ungulates located in different regions of the mountain may have different habitat requirements and not equally sensitive to climate change (Chen, Hill, Ohlemuller, Roy, & Thomas, 2011).
Although our future models detected retention of suitable habitat in parts of the northern region, we still observed a shift of the species distribution northward (latitudinal) and possibly to higher elevations propelled by climate change, as the southern region experienced a major decline of suitable habitat. Previous research showed the spatial response to climate change for a variety of species-shifting distributions northward and upward (Alvarez, Salas, Harings, & Boykin, 2017;Ogawa-Onishi, Berry, & Tanaka, 2010;Popy, Bordignon, & Prodon, 2010;Root et al., 2003;Wilson et al., 2005). In our results, about 20% (2,500 km 2 ) of the average habitat gains in 2070 showed a shift northwards. In contrast to Luo et al.

| Uncertainty and robustness of future projections
Three have been objections by a handful of studies about the use of climatic variables alone to model the climate envelope of species (Geyer, 2011;Sax, 2007) because there are a wide range of climate change-related stresses that are at play that could affect population ecology and physiology. Climatic variables alone may not always correctly ascertain current or future species distribution over space (Salas, Seamster, et al., 2017). Despite the limitations of models using climatic variables solely, several proponents have praised their crucial role to provide broad insights on the future effects of climate change on species distribution (Guisan & Thuiller, 2005;Soberón, 2007 (Soberón, 2007), were not integrated into the modeling process due to data unavailability.
Furthermore, when we projected present conditions to future years, we assumed that the land cover and other topographic factors in the region would remain stable in the future. These shortcomings of our models could result to spatial mismatch between the projected future and real distributions of argali habitats.
Results of our study may not deliver the most accurate predicted distributions; however, they could still facilitate in Marco Polo argali conservation planning in the region and provide approximation of changes of species distribution in the context of climate change.
In terms of our modeling techniques, multiple limitations associated with the projection of species distributions into the future under different climate scenarios have been documented. Three broad categories of uncertainties affecting the climate variables used to drive the ENMs include (a) uncertainties in future greenhouse gas concentrations (Meinshausen, 2011), (b) limitations in the accuracy of GCM-simulated, large-scale physical climate responses to changing greenhouse gas levels (Knutti & Sedláček, 2013), and (c) shortcomings and assumptions inherent to statistical downscaling methods used to refine GCM results to a finer level of spatial detail (Barsugli, 2013;Dixon, 2016). By utilizing data products derived from four GCMs and two RCPs, this study partially explores two of these three sources of climate variable uncertainty. Stoklosa, Daly, Foster, Ashcroft, and Warton (2015) specifically discuss approaches to account for some uncertainties in the climate variables used to drive ENMs. Furthermore, several authors have shown variability among future projections of suitable climatic conditions when using different climate models applied to the same species occurrence dataset (Bakkenes, Alkemade, & Ihle, 2002).
Variability of results makes assessment of projections of future conditions a complex effort. First, it is not possible to determine which single ENM could provide the most accurate information for a species, although one could argue that the climate model with the highest accuracy in capturing the present-day suitable habitat conditions may produce more accurate future projections. However, Thuiller (2004) reasoned that even when a model has the highest AUC and K statistics, the model may not provide the best estimate of the future distribution of suitable conditions, especially as every model is based on different assumptions. It is most fitting to use an ensemble model of future projections, as this ensemble represents the areas of agreement among individual model projections. The reliability of future conditions produced by ensembles may be questioned, but ensemble results do incorporate the positive aspects of multiple models and provide a more conservative assessment of these conditions.

| CON CLUS ION
The main contribution of this study is the quantification of the possible impact of climate change regimes on the availability of future suitable habitats for Marco Polo argali. Results in all our models showed major losses of suitable habitat for the populations of argali in eastern Tajikistan. While our results are potential projections of future Marco Polo argali suitable habitat, this is the first quantitative assessment of habitat shifts of argali under future climate change in eastern Tajikistan, which could provide information for conservation planning for this near threatened species in the region. There is a need to conduct a similar study to include habitats in adjacent populations in Afghanistan, China, Pakistan, and Kyrgyzstan and to determine the extent of changes in critical habitats in climate scenarios. Ltd., who assisted in the surveys, and their hospitality are greatly appreciated.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N
Eric Ariel Salas conceived and designed the research; performed the modeling; analyzed and interpreted the data; contributed reagents, materials, analysis tools or data; and wrote the manuscript.
Raul Valdez interpreted the data and contributed reagents, materials, analysis tools, or data. Stefan Michel contributed reagents, materials, analysis tools, or data. Kenneth Boykin contributed reagents, materials, analysis tools, or data.

O R C I D
Eric Ariel L. Salas http://orcid.org/0000-0001-9019-730X Reducing uncertainty in projections of extinction risk from climate