Using species distribution modelling to guide survey efforts of the snow leopard (Panthera uncia) in the Central Kyrgyz Ala-Too region

Using species distribution modelling to guide survey efforts of the Snow Leopard (Panthera uncia) in the Central Kyrgyz Ala-Too region. — V. Tytar, T. Asykulov, M. Hammer. — Listed as Vulnerable (IUCN 2017), the snow leopard is declining across much of its present range. One of the major reasons for the snow leopard population decline in the last two decades is a reduction in large prey species that are the cornerstone of the conservation of the snow leopard; in the Central Kyrgyz Ala-Too region such species is primarily the Siberian ibex (Capra sibirica). Understanding factors affecting basic requirement of ibex and shaping its distribution is essential for protecting the prey species snow leopards rely on the most. Using a niche modelling approach we explored which environmental features are best associated with ibex occurrence, how well do models predict ibex occurrence, and does the potential distribution of highly suitable ibex habitat correlate with records of snow leopard. A PC analysis was used to capture aspects of ibex ecology and niche. Results of such analysis agree with the herbivore character of the species and bioclimatic habitat requirements of the vegetation it feeds upon, richer in flatter areas, and where plants may benefit from more sunlight. The niche model based on maximum entropy (Maxent) had “useful” discrimination abilities (AUC = 0.746), enabling to produce a map, where a contour line is drawn around areas of highly predicted probability (> 0.5) of ibex occurrence. In terms of nature conservation planning and setting snow leopard research priorities these areas represent the most interest. With one outlier, most of snow leopard records made in the study area (n = 15) fell within the 10 percentile presence threshold (0.368). Predicted probability of ibex occurrence in places where records were made of snow leopard presence (pugmarks, scrapes etc.) was 0.559 expectedly suggesting areas of high ibex habitat suitability attract the predator.


Introduction
The snow leopard (Panthera uncia) is an icon for conservation in the mountain regions of Asia. As a top-order predator, its presence and survival is also an indicator of intact, "healthy" eco-region. Snow leopards are listed as Vulnerable in the IUCN Red List (IUCN 2017) and its abundance is declining across much of its present range.
One of the major reasons for the snow leopard population decline in the last two decades is a reduction in prey resource base. A recent investigation into prey preferences of the snow leopard revealed three key large prey species that are cornerstone of the conservation of the snow leopard globally; one of these species is the Siberian ibex (Capra sibirica) (Lyngdoh et al., 2014). According to these authors, Siberian ibex were killed by snow leopards wherever they occurred, meaning the species is a vital resource for the predator.
Siberian ibex range is spread across the mountains of Pakistan, China, India, Afghanistan, Kyrgyzstan, Kazakhstan, Uzbekistan, Mongolia, Russia and Tajikistan. Ibex mainly occupies rocky mountainous regions, both open meadows and cliffs, coming down to low elevations during winter (Fedosenko, Blank, 2001). The species avoids densely forested areas and prefers to remain near steep and escape terrain (Fedosenko, Blank, 2001). The ibex is crepuscular in feeding, foraging in evenings and mostly in early morning hours (Fedosenko, Blank, 2001). They come down from their steep habitats during late afternoon and evenings to the alpine meadows below to feed.
Ibex lives in small groups (6-30 animals) varying considerably in size, rarely in herds of >100 animals (Fedosenko, Blank, 2001). In the study area of the Central Kyrgyz Ala-Too groups hardly exceed 20 individuals. Any snow leopards in this area would depend on this species as a primary food source, as far as wild sheep (argali) seem to occur here very occasionally.
The presence of species depends upon the specific environmental conditions that enable it to survive and reproduce (Marzluff, Ewing, 2001). Understanding the factors influencing its existence is a basic requirement for the assessment of the species distribution and devising efficient species conservation strategies (Wein, 2002). This knowledge could help us to focus our efforts on protecting the prey species snow leopards rely on the most. Additionally, Siberian ibex suffer from poaching and trophy hunting, decreasing population sizes even further. In order to protect these animals, there is a need to identify and preserve important areas for wildlife that could address the issues of over-grazing and poaching.
Species conservation usually focused on the most suitable habitat for a species of concern, but the challenge is to identify high-quality habitat across large areas (Bellis et al., 2008). Among the various tools used in conservation planning to protect biodiversity, species distribution models (SDMs), also known as climate envelope models, habitat suitability models, and ecological niche models provide a way to identify the potential habitat of a species in an ecoregion and their applications have greatly increased. SDMs are based on the concept of the "ecological niche" (Hutchinson, 1957), which can be defined as the sum of the environmental factors that a species needs for its survival and reproduction. Many niche models are based on climate variables because these data are readily available, covering large spatial scales. SDMs predict the potential distribution of a species by interpolating identified relationships between presence/absence or presence-only data of a species on one hand and environmental predictors on the other hand across an area of interest. From the array of various applications, Maxent  stands out because it has been found to perform best among many different modeling methods (Elith et al., 2006). Maxent is a maximum entropy based machine learning program that estimates the probability distribution for a species' occurrence based on environmental constraints . It requires only species presence data and environmental variable (continuous or categorical) layers for the study area.
Here, we addressed the following questions: (1) which environmental features are best associated with ibex occurrence? (2) How well do models predict the occurrence of this species? (3) What is the potential distribution of highly suitable habitat for ibex in the study area and does it correlate with records of snow leopard?

Study area
The chosen study area is located in the southern Kyryz Ala-Too, away from the main cities in the north. Surveys were centered around the Karakol Mountain Pass (3,452 m) and encompassed areas in the upper reaches of the West and East Karakol rivers. The expedition's main access route into the area was the Suusamyr plateau, a high steppe plateau (2,200 m) that although only some 160 km from the capital city of Bishkek, is also one of the more remote and rarely visited regions of Kyrgyzstan. The base camp was located close to the Suusamyr-Kochkor road approximately in the middle of the planned study area (42.359535 o N, 74.737829 o E, 3002 m a.s.l.). From the base camp mostly one-day surveys, but also some two-day/one-night surveys were conducted to various portions of the Kyrgyz Ala-Too Range

Species records
Between 2014 and 2018, teams of citizen scientists recruited by Biosphere Expeditions and NABU (Kyrgyzstan) gathered 103 records of Siberian ibex in the study area according to direct sightings and camera trap results, and 15 records of snow leopard presence.
The extracted points were georeferenced using a Garmin eTrex 10 GPS receiver in combination with GoogleEarth. All coordinates were expressed in decimal degree and converted to a point vector file for modeling the distribution of the species. Only spatially unique ones, corresponding to a single environmental grid cell (resolution of 30 arc seconds, ~ 1 km) were used.

Environmental and bioclimatic data
To relate the occurrence records of ibex with abiotic conditions, we downloaded 19 bioclimatic variables for the current climate at a 30 arc second resolution and WGS84 coordinate reference system (Hijmans et al., 2005). These variables represent annual trends of temperature and precipitation, seasonality, and extreme or limiting environmental factors (e.g., temperature of the coldest and warmest month).
Temperature has long been recognized as an important environmental factor in ecosystems in regard to its pivotal role over biological (development, growth and reproduction), chemical, and physical properties. Precipitation regimes and variation of precipitation events have broad effects on ecosystem productivity, habitat structure, and ultimately on species' distribution.
For the study region of the Kyrgyz Ala-Too scientific data on range of environmental resources (other than bioclimatic) are limited which hinders sustainable management and nature conservation.
The need for update information has long been recognized and stimulated the use of earth data using remote sensing techniques, which has become a universal and familiar instrument for assessing natural resources (Philipson, Lindell, 2003). Information from low-altitude satellite sensors and remote sensing offer an optimal path for understanding pattern and process related to rangeland condition in the area. The multi-temporal and multi-spectral data acquired by various satellite sensors are used to identify, map and monitor rangelands, derive specific environmental variables.
We used a Landsat 8 satellite image (path 151⁄row 31) taken on 7th of August 2014, freely acquired from the U.S. Geological Survey georeferenced GeoTIFF files at a 30 m resolution via the EarthExplorer website (https://earthexplorer.usgs.gov). This image encompassing the study area was selected because of the minimum cloud coverage (0.16 %).
Candidate predictor variables most suitable to predict the ecological niche dominance within the landscape include tasselled cap transformation, enhanced vegetation index (EVI), and Land surface temperature (LST). EVI is a standardized vegetation index which allows us to generate an image showing the relative biomass. Landsat-8 thermal bands i.e., Band 10 and Band 11, were used to calculate the brightness temperature over the study area. This gives an assessment of the ground temperature that may be hotter than the ambient air temperature.
Tasselled cap transformations, originally developed to understand changes in agricultural lands, generate three orthogonal bands from the six-band Landsat composite (Huang et al., 2002). The three generated bands represent measurements of brightness (band 1, dominated by surface soils), greenness (band 2, dominated by vegetation and correlates with EVI), and wetness (band 3, includes interactions of soil, vegetation and moisture patterns) (Kauth, Thomas, 1976). A digital elevation model (DEM) was used as input for capturing topographic variables. The DEM was aggregated from the 30 seconds (~30 m) NASA Shuttle Radar Topographic Mission (SRTM) DEM (http://srtm.csi.org). The following terrain features were extracted: slope, aspect and terrain ruggedness. Slope is the steepness or the degree of incline of a surface. The direction a slope faces with respect to the sun (aspect) has a profound influence on vegetation and snowpack. We split aspect into two components: eastness = sin(aspect) and northness = cos(aspect) These indices of northness and eastness provide continuous measures (−1 to +1) describing orientation. The topographic ruggedness index (TRI) was developed to express the amount of elevation difference between adjacent cells of a DEM. These were selected because generally, terrain roughness and slope create a template of risk, in which herbivores have to trade off between resource acquisition (e.g. foraging in high quality habitats, finding mates) and predator avoidance (Schweiger et al., 2015). Ibex are very good climbers that find protection from predators and the possibility to overview large areas in predominantly rocky terrain with steep slopes.
The resolution or grain of Landsat images (30 m) and the DEM is finer than the accuracy by which we can record ibex presence in the field; for this reason all the considered environmental layers have been rescaled to a 30 arc second resolution (~ 1 km).
SAGA GIS software (v. 2.2.7) has been used for the preliminary data processing, extracting (clipping) images for the study area (Conrad et al., 2006).

Statistical modelling
Factor analysis in Statistica 10 was used to examine the contributions and the main patterns of inter-correlation among the potential environmental controls. Principal component (PC) was used as the extraction method. By rotating the factors a factor solution was found that is equal to that obtained in the initial extraction but which has the simplest interpretation, and for this purpose the Varimax normalized type of rotation was applied. Usually a solution that explains 75-80% of the variance is considered sufficient.
Species distribution models obtained from a large data set of associated environmental covariates often inherently result in multi-collinearity, a statistical problem defined as a high degree of correlation among covariates. Principal component analysis is among the statistical procedures proposed to solve or to reduce multi-collinearity because the obtained PCs are independent of each other, that is, they are orthogonal. Accordingly, this permits them to be used as independent, noncorrelated variables in analyses of modeling potential species distribution (Cruz-Cárdenas et al., 2014) instead of the raw data of environmental variables.

Maxent distribution model
We used the freely available Maxent software, version 3.3.3k, which generates an estimate of probability of presence of the species. The default settings of Maxent were used in this study. We ran models with 25 bootstrap replicates.
Model performance was assessed using the average AUC (area under the receiver operating curve) score to compare model performance. AUC values > 0.9 are considered to have 'very good', > 0.8 'good' and > 0.7 "useful" discrimination abilities (Swets, 1988). The logistic output format was used, because it is easily interpretable with logistic suitability values ranging from 0 (lowest suitability) to 1 (highest suitability). Applying a threshold is the last step of many species modelling approaches. We used 10 percentile training presence because this threshold value is considered to provide a better ecologically significant result when compared with more restricted thresholds values (Phillips, Dudík, 2008)

Results
The PC analysis provided a comprehensive way to analyze the niche of Siberian ibex in the study area and captures various aspects of ibex ecology and environmental requirements of the species. PC1-5 extracted from all the variables explained ~88 % of the variance.
The first two principal components heavily correlated with temperature/altitude and greenness (50.2 % contribution of PC1), and precipitation (21.5 % -of PC2); the following three components collectively explained 16.3 % of the variation and correlated with one or two parameters: PC3 was negatively associated with the inter-correlated slope and topographic ruggedness index (TRI), PC4with the mean temperature of the wettest quarter, and PC5 -negatively with eastness.
These findings are in good agreement with the herbivore character of the species and bioclimatic habitat requirements of the vegetation it feeds upon, richer in flatter areas, and where plants may benefit from more sunlight.
Gridded data layers produced in SAGA GIS representing the first five principal components were then used for modelling species distribution. From the 25 model runs, the average AUC was 0.746 (meaning "useful" discrimination abilities of the final model), with little variation in AUC between runs (SD = 0.016). Рис. 1. Топографічна карта дослідженої ділянки Киргизького Ала-Тоо; контурна лінія обведена навколо териорій, де прогнозована ймовірність перебування сибірського гірського козла складає понад 0,5; трикутники вказують місця реєстрацій перебування снігового барса.     The averaged output from these model runs is shown on the map (Fig. 1), where a contour line is drawn around areas of predicted probability of Siberian ibex occurrence exceeding 0.50. In terms of nature conservation planning and setting snow leopard research priorities these areas of high predicted probability of Siberian ibex occurrence represent the most interest.
With one exception, most of the snow leopard records made in the study area (n = 15) fell within the 10 percentile training presence logistic threshold of 0.368, and there were two outliers when the 0.50 threshold was applied. One of them was a record made at the Karakol Pass, which the animal had used for crossing from one mountain range to the other.
On average the predicted probability of Siberian ibex occurrence in places where records were made of snow leopard presence was 0.559, with a minimum of 0.526 for the bulk of the records, expectedly suggesting areas of high habitat suitability for Siberian ibex attract the predator (figs 2-6 are from such areas).

Conclusion
Further research is needed to confirm wider snow leopard presence and monitor snow leopard and prey population trends in the survey area. Presence/absence surveys will need to be repeated in the coming years, using camera traps from the very beginning of the survey. Finding a trail and/or relic scrape(s) is a high priority. If either more of these can be found, remote camera-trapping would be enhanced as a survey tool. These efforts can be guided by modelling exercises as above, showing places where basic requirements for Siberian ibex, upon which snow leopards rely the most, are met to a significant degree.