On the snow leopard Trails: Occupancy pattern and implications for management in the Pamir.

The snow leopard (Panthera uncia) inhabits one of the most challenging environments on Earth, referred to as the 'third pole'. Only a fraction of its vast range has been explored thus far, owing to myriad of barriers inflicted by the remote terrain and socio-ecological realities of the landscapes. Understanding distribution patterns of species is essential to devise practical management measures. This study aimed to understand the distribution pattern and factors influencing occupancy of snow leopard in the Pamir Mountain range through sign-based occupancy modelling. Our study confirmed that the Pamir range is a snow leopard stronghold, with occupancy estimated at 0.57 ± 0.02. The topographic features positively influenced the detection probability (p = 0.37 ± 0.005) of snow leopards. Occupancy was influenced by mean annual precipitation (β = -6.12 ± 1.8), density of roads (β = -1.61 ± 0.6) and water sources (β = 0.74 ± 0.4). Our findings underpin that sign-based distribution surveys provide vigorous scientific knowledge about elusive species and merit replication being used for other species. We propose to redefine the protected area boundaries based on ecological knowledge and encourage transboundary cooperation to safeguard snow leopards at a landscape scale.

Introduction

Conservation of elusive species and associated complex landscapes largely depends on knowing where they occur. Patterns of habitat use by species explain the interactions between individuals of the same or different species, the overall population stability, and the functioning of the ecosystem (Guo et al., 2017). However, wildlife habitats are deteriorating and rapidly changing due to the increasing conventional i.e., anthropogenic and emerging i.e., climate change, threats. Consequently, distributions and habitat use patterns of species are altering, sometimes in unpredictable ways (Sanderson et al., 2002). Given this scenario, knowledge about the species distribution and habitat use is crucial (MacKenzie et al., 2017) as such efforts provide the logical and scientific basis to formulate informed management strategies and action plans (Liu et al., 2013).

Among wildlife species, large carnivores are more vulnerable to the waves of change in ecosystems, owing to their low densities, and large home ranges (Ceballos et al., 2005). Climatic and human induced changes coupled with the thin population and loss of suitable habitat have affected global distribution of many carnivore species and many of them are at the risk of extinction (Din et al., 2019, Michalski and Peres, 2005, Simberloff, 1998).

The snow leopard (Panthera uncia) is a keystone species in high-altitude ecosystems of Central and South Asia (Snow Leopard Network, 2014a). The species occurs across 12 countries of South and Central Asia, over a predicted area of 1.8 to 3.26 million km2 (Nyhus et al., 2016, Snow Leopard Working Secretariat, 2013) in high mountain ranges of Sayan, Altai, Tien Shan, Kunlun, Pamir, Hindu Kush, Karakoram, and Himalayas (Fox and Chundawat, 2016, Jackson et al., 2010, Thompson, 2013). However, ecological and geopolitical barriers such as vast and remote range, elusive nature, logistic constraints, political instability, trans-border issues, and limited provisions for research have limited our understanding of the snow leopard ecology (Hunter et al., 2016). Consequently, as of today, we are certain about snow leopard’s occurrence in only 27% of its total range (McCarthy et al., 2016), and population estimates are available for 2% of the distribution range (Snow Leopard Working Secretariat, 2013). Site occupancy modelling using presence/absence data has emerged as a go-to method in monitoring of difficult-to-detect species such as snow leopards (Fleishman et al., 2001, MacKenzie et al., 2002). However, snow leopard specific studies are meagre and have been conducted rather sporadically within its vast range including China (Alexander et al., 2016, Li et al., 2013, Wolf and Ale, 2009), Nepal (Aryal et al., 2016), India (Bhatnagar et al., 2016, Watts et al., 2019), and Russia (Kalashnikova et al., 2019).

Understanding of the ecology of snow leopards at landscape scale is, fundamentally, still in its nascent stages and the available literature largely focusses on the human-snow leopard conflicts, diet and climate change (Robinson and Weckworth, 2016). The recent controversial (Ale and Mishra, 2018) down listing of snow leopards by the IUCN (McCarthy et al., 2017) is likely to hamper securing of resources required for further explorations. The snow leopard range in Pakistan spreads across the Hindu Kush, Karakoram, Himalaya and Pamir mountain ranges of northern Pakistan, over an area of 80,000 km2 (Ahmad et al., 2016, Din and Nawaz, 2011, Hussain, 2003). About 60% of the total range falls in Gilgit-Baltistan and potential habitat is concentrated in and around the Pamir and Karakoram mountain ranges(Snow Leopard Network, 2014a). The trans-border region of the Pamir mountain range is considered as a biodiversity hotspot and an important corridor for snow leopards (Din et al., 2017; Din et al., 2019, Hameed et al., 2020). This region is a candidate site for transboundary cooperation between Afghanistan, China, Pakistan, and Tajikistan to manage biodiversity at a regional scale (Moheb and Paley, 2016).

This study focused on the Pamir range, aiming at assessing occupancy of snow leopards at the large landscape scale (>10,000 km2). Since, snow leopards are solitary and territorial animals and mark their territory through scrapes and scent marks (Jackson and Hunter, 1996), we considered these traits to assess the site use and detection rate of snow leopards by conducting grid-based sign surveys, in this study. We also investigated environmental and anthropogenic factors that influence occupancy of snow leopards in northern Pakistan.

Materials & methods

Study area

This study area falls in the Pamir mountain range, one of the most magnificent mountain ranges in the world and verged by the mighty mountain ranges of Hindukush and Karakoram in Pakistan and is fondly known as the ‘roof of the world’ (Bliss, 2006; Din et al., 2017). Due to its exceptional geographic range, socio-cultural and ecological diversity, the Pamir is considered a biodiversity hotspot (Din et al., 2019) and listed among the world’s 200 eco-regions (Olson et al., 2001). Being situated in the alpine zone, the winter is harsh, autumn is mild and summers are pleasant with annual temperature and precipitation ranges from −18 °C to 23 °C and 50 mm to 300 mm, respectively (Khan et al., 2014). Our study sites encompassed three high altitude National Parks (NPs) viz; the Khunjerab National Park (KNP), Qurumbar National Park (QNP), Broghil National Park (BNP) and Community Managed Conservation Areas (CMCAs) of Misgar-Chipursan Valleys falling in the Karakorum-Pamir and Hindu Kush-Pamir mountain ranges, respectively (Fig. 1).

Fig. 1

Map of the study site showing overall study area, survey blocks and occupancy points surveyed.

The KNP (74° 52′- 76° 02́E, 36° 56′- 36°13́N) covers an area of 4,455 km2 and ecologically continues to Taxkorgan Nature Reserve, China. It was notified as a National Park in 1974 to protect the Marco Polo sheep (Ovis ammon polii) and snow leopard (Panthera uncia). The KNP and its buffer zone also supports a good population of ibex (Capra sibirica) and blue sheep (Pseudois nayaur), where a community-based trophy hunting program is being implemented (Khan et al., 2014).

The BNP spread across 1,348 km2 in Chitral District and lies in between 73° 13́- 73° 52́E and 36° 42́- 36°55́N and was notified as a National Park in 2010 to safeguard high altitude (3,300–4300 m) wetlands and associated wildlife. It borders the North and West with the historical Wakhan corridor of Afghanistan, in the East and South borders with Yasin valley and Qurumbar National Park of Gilgit-Baltistan (GB), respectively (Shah, 2012).

The QNP (73°40́-73°05́E, 36°40́-36°55́N) lies in the Ghizer District and borders with the Wakhan Corridor of Afghanistan in the North, the Broghil National Park in the West, the Hunza Valley in the East and Ghizar River in the South. The park extends over 740 km2 area and was notified in 2011 to conserve the natural flora and fauna, wetlands and improve livelihood of the local people (WWF-Pakistan, 2016). The QNP and its buffer zone, the community managed conservation area falls in the mountain desert ecosystem where the average rainfall rarely exceeds 150 mm.

The CMCAs of Misgar and Chipursan Valleys (74° 0́-74° 5́E, 37° 0́-37° 5́) spread across 2,248 km2 in the Pamir and are surrounded by KNP in the Northeast, Afghan and Tajik Pamir in the North, and QNP in the West. The Pamir and the Karakoram Mountain ranges meet at Kilik Pass of Misgar Valley on the borders of China, Afghanistan, and Pakistan. The Kilik and Mintika watersheds used to be potential habitat of Marco Polo sheep (Ali et al., 2019) which constitutes one of the major wild prey of snow leopard in the region.

Study design and data collection

Snow leopards are elusive, solitary and territorial animals and mark their territory through scrapes and scent marks. We considered these traits to assess the site use and detection rate of snow leopard in this study, undertaken between the years 2010–12. Due to the logistic challenges imposed by the remote terrain, we divided the landscape into four (KNP, Misgar-Chipursan, KNP and BNP) larger study units (survey blocks). Next, we further split our study area into smaller sampling units by overlaying 5 × 5 km grid cells (sites) using ArcGIS v.10.7. We surveyed each survey block separately. In each accessible site, a team of two observers surveyed the most likely micro-habitats (points) such as ridgelines, saddles, gorges, cliff-bases where snow leopards usually tend to leave signs. Since, the snow leopard landscape restricts free movement, we confined each point surveyed to a rectangular route of 5 × 100 m to ensure uniformity of the survey effort across the study sites. Only one fresh sign, scrape or scent spray site (Jackson, 1996) was recorded at each point and a maximum of 12 points were recorded in each site. Each point was marked with GPS and habitat attributes (survey covariates) were noted.

We surveyed 109 accessible grid cells (KNP = 21, BNP = 14, QNP = 24, CMCAs = 50) covering 2,725 km2 area, on altitude ranging from 2,257 m to 4,837 m (average = 3,737 m). About half of the grid cells (n = 533) in the study area fall in the high altitude (>5,000 m) cold desert glaciated terrain and hence, were therefore not accessibly by people.

Sampling of covariates

Based on the field observations and previous research (Ghoshal et al., 2019, Henschel et al., 2016), we hypothesized 11 site and 7 survey covariates as predictors of snow leopard occupancy in the study area (Table 3). Topographic covariates such as elevation was downloaded from the website (www.earthexplorer.com) while, slope was derived from the elevation. Density of rivers (water) in the study sites was calculated using Spatial Analyst Tool in ArcGIS v.10.7. (ESRI, Redlands. USA). Similarly, ruggedness was extracted using the geo-processing script (Sappington, Longshore & Thompson, 2007). The climatic covariates (annual mean temperature and annual mean precipitation) were downloaded from www.worldclim.org developed by Hijmans et al. (2005). We extracted the anthropogenic covariates (density of roads and settlements) following a previously reported study (Luo et al., 2019). We collected survey covariates based on the pre-categorization of micro habitats and observer’s experience. We categorized “topographic habitat variables” into valley, ridgeline, cliff base, plateau; “physiognomic vegetation variables” into pasture, scrub, barren; and “terrain brokenness” into flat, slightly, moderately, and very broken, respectively.

Table 3

Estimates of β coefficient values for covariates (derived from the top model) assumed to influence snow leopard site use.

Coefficients	β(s)	SE
ψ (Intercept)	−1.06	0.60
ψ (Mean annual precipitation)	−6.12	1.80
ψ (Mean density of roads)	−1.61	0.57
ψ (Mean density of rivers)	0.74	0.42
ψ (Mean ruggedness)	0.68	0.47
p (Intercept)	−1.25	0.24
p (Barren)	0.78	0.29
p (Scrub)	1.76	0.30

Data analysis

We standardized all the covariates to mean zero (0) and variance of one (1) using psycho package in R (R Core Team, 2019), before running the occupancy models. We also ran Pearson correlation test to exclude highly correlated (r > 0.70) variables (Schober and Schwarte, 2018) in the analysis (Alexander et al., 2016).

Having the data consolidated, we developed detection histories as “1s” and “0s” in a matrix of sites vs. points (MacKenzie et al., 2002) and analyzed using software PRESENCE v.12.37 (Hines, 2006). Since, we had surveyed each of the four study blocks once and interested to estimate probability of site use rather than true occupancy of snow leopard, we used the single species-single seasons occupancy model (MacKenzie et al., 2002). On the onset, we ran the null model by keeping both detection (p) and occupancy (ψ) constant as we were interested to model covariates assumed to influence detection probability first, while, keeping occupancy constant and to use best detection model based on the AIC weight in further modelling aimed to investigate significance of occupancy covariates (Henschel et al., 2016, Karanth et al., 2011). We continued with the univariate, bivariate and then multivariate modelling, unless finding the model that best explained the variation in detection and markability, and occupancy (site use) of snow leopard at the site level (Henschel et al., 2016). The best fitting model was determined using the Akaike Information Criteria (AIC). The model that has the best fit (likelihood) obtains the minimum value of AIC value (Burnham and Anderson, 2002). Finally, we used the output of the best model to predict the probability of unconditional occupancy (site use) at the sites (girds) not sampled (Kshettry et al., 2017) to have holistic overview of the study area. We pooled the occupancy values obtained for each site to calculate the average occupancy estimates with standard errors across the study blocks and study area.

Results

We surveyed 109 accessible grid cells (KNP = 21, BNP = 14, QNP = 24, CMCAs = 50) covering 2,725 km2 area, on altitude ranging from 2,257 m to 4,837 m (average = 3,737 m) About half of the grid cells (n = 533) in the study area fall in the high altitude (>5,000 m) cold desert glaciated terrain and hence, were humanly inaccessible.

We surveyed 822 sites, an average of 7 sites/grid (Fig. 1). Among them, snow leopards were detected in half sites (4 1 3), which corresponded to a naïve occupancy estimate (proportion of sampled units in which signs were detected) of 0.39. Accounting for imperfect detection resulted in a mean detection probability of 0.37 (SE = 0.005) and an average occupancy estimate of 0.57 (SE = 0.02), which is 18% higher than the naïve estimate. Snow leopard occupancy was higher in KNP, followed by QNP, Misgar-Chipursan (MC) Valleys and BNP, respectively (Fig. 2 and Fig. 3).

Fig. 2

Average occupancy estimates and detection probability in different blocks of the study site (BNP = Broghil National Park; KNP = Khunjerab National Park; MC = Misgar-Chipursan; QNP = Qurumbar National Park).

Fig. 3

Distribution of snow leopard occupancy probabilities in the study sites.

Our final model retained four survey (mean annual precipitation, road density, river density and mean ruggedness) and two site (barren and scrub) covariates (Table 1). Snow leopard occupancy at site level reduced with increase in mean annual precipitation and density of roads, whereas increased with a corresponding increase in the availability of water sources. Similarly, the detection probability was higher in barren and scrub areas as compared to pastures and plateaus (Table 2).

Table 1

Eleven site and seven survey covariates hypothesized to affect snow leopard occupancy and detection probability across the study sites.

Covariate	Data source	Likely effect on occupancy	Published evidence of the effect of the variable
Site Covariates
Precipitation	www.worldclim.org (Hijmans et al. 2005)	−	Bai et al., 2018
Road	Spatial Analyst Tool in ArcGIS v.10.7. (ESRI, Redlands. USA)	−	Luo et al., 2019
Ruggedness	Geoprocessing script (Sappington, Longshore & Thompson, 2007)	+	Taubmann et al., 2016
River	Spatial Analyst Tool in ArcGIS v.10.7. (ESRI, Redlands. USA)	+	Aryal et al., 2016
NDVI	Spatial Analyst Tool in ArcGIS v.10.7. (ESRI, Redlands. USA)	+	Taubmann et al., 2016
Elevation/Altitude	www.earthexplorer.com	−	Taubmann et al., 2016
Temperature	www.worldclim.org (Hijmans et al. 2005)	−	Farrington and Li, 2016
Slope	www.earthexplorer.com (Derived from elevation)	+	Alexander et al., 2016
Settlements	Spatial Analyst Tool in ArcGIS v.10.7. (ESRI, Redlands. USA)	−	Luo et al., 2019
Prey	Snow Leopard Foundation	+	Ghoshal et al., 2019
Staff	Collected in this study (Wildlife Department GB)	+	Henschel et al., 2016
Survey Covariates
Pasture	Assessment of microhabitats in this study	−	Fox and Chundawat, 2016
Scrub	Assessment of microhabitats in this study	+	Fox and Chundawat, 2016
Barren	Assessment of microhabitats in this study	+	Snow leopard Network, 2014
Terrain	Assessment of microhabitats in this study	+	Snow leopard Network, 2014
Ridge	Assessment of microhabitats in this study	+	Snow leopard Network, 2014
Cliff base	Assessment of microhabitats in this study	+	Snow leopard Network, 2014
Plateau	Assessment of microhabitats in this study	−	Snow leopard Network, 2014

Table 2

Results of multivariate model selection for estimating snow leopard occupancy (ψ) in the study sites.

Model	K	AICc	ΔAICc	AIC wt	Moddel Likelihood	logLik
Ψ (Precipitation + Road + River + Ruggedness), p(Scrub + Barren)	8	558.78	0	0.28	1.00	542.78
Ψ (Precipitation + Road), p(Scrub + Barren)	6	559.68	0.9	0.18	0.64	547.68
Ψ (Precipitation + Road + Ruggedness), p(Scrub + Barren)	7	560.24	1.46	0.14	0.48	546.24
Ψ (Precipitation + Road + NDVI), p(Scrub + Barren)	7	561.34	2.56	0.08	0.28	547.34
Ψ (Precipitation + Road + Aspect), p(Scrub + Barren)	7	561.45	2.67	0.07	0.26	547.45
Ψ (Precipitation + Road + Altitude), p(Scrub + Barren)	7	561.6	2.82	0.07	0.24	547.6
Ψ (Precipitation + Road + Temperature), p(Scrub + Barren)	7	561.65	2.87	0.07	0.24	547.65
Ψ (Precipitation + Road + Settlement), p(Scrub + Barren)	7	561.67	2.89	0.07	0.24	547.67

Discussion

This study was the first ever attempt to ascertain the distribution and factors influencing the occupancy and detection of snow leopard in the vast landscape of Pamir using robust scientific approaches. Our understanding of the snow leopard’s ecology is still in its nascency (Robinson and Weckworth, 2016) and hardly 30% (3.26 million square kilometers) of potential range has been assessed so far (McCarthy et al., 2016). Large and remote habitat, the elusive nature of the species coupled with the geopolitical and bureaucratic barriers (Hunter et al., 2016) are some factors limiting ecological research across the snow leopard range. Our study confirms the occurrence of snow leopard in myriad of protected area network in the Pamir bordering with China and Afghanistan, and thus adds to the limited body of knowledge on snow leopard’s ecology.

Among the study sites, the occupancy of snow leopard was higher in KNP, followed in MC Valleys, QNP and BNP. KNP is the oldest National Park in the snow leopard range in Pakistan, established in 1975 to protect Marco Polo sheep and snow leopard as an outcome of the exploratory surveys by Schaller in 1970 s (Schaller, 1980).The habitat features, diversity of wild prey (Khan et al., 2016) and protection measures in place may have contributed to higher occupancy estimates of snow leopard in the national park. The MC valleys are contiguous to the KNP and share similar habitat features such as rugged and broken rocky outcrops and scrub patches, preferred by snow leopards (Fox and Chundawat, 2016), while QNP and BNP represent high-altitude desert and wetland ecosystems.

The MC Valleys were recently notified as Community Managed Conservation Areas (CMCAs), a lower IUCN PA category (VI) than the national park (Dudley et al., 2013); primarily to promote trophy hunting of the ibex. Recent studies (Din et al., 2017; Din et al., 2019) revealed high intensity of human-snow leopard conflict in the MC Valleys and the current PA status cannot safeguard the survival of snow leopard in these valleys. We recommend the denotification of the CMCAs and their annexure with KNP by expanding the boundary of the park to ensure survival of snow leopards in this region. The spillover effect resulting from the proposed management measures will likely enhance the distribution of snow leopards in the adjacent national parks of QNP and BNP as well.

Our results reinstate that the habitat features such as ruggedness and availability of water resources best explain the distribution of snow leopards in the study area. These findings are consistent with the studies conducted in other parts of the snow leopard range (Bai et al., 2018, Liu et al., 2013, McCarthy et al., 2005, Wolf and Ale, 2009) affirming that habitat features of snow leopard remain constant across its range. Moreover, precipitation and anthropogenic activities were found to have negative impact on the habitat use by snow leopards. Recent studies (Bai et al., 2018, Li et al., 2016) on snow leopards have also reported similar findings. The remote mountain communities in the study sites practice transhumance and with the passage of time the developmental needs are increasing posing threat to snow leopards, its wild prey and habitat (Khan et al., 2014; Din et al., 2017). We suggest that the shared ecosystem services of the snow leopard landscape (Din et al., 2020) should be realized in formulating developmental agenda for this region. Significant portion of Pamir Mountain Range falls in the 20 model landscapes identified under the Global Snow leopard and Ecosystem Protection Program (GSLEP) to achieve the global goal of “securing 20 by 2020” (Snow Leopard Working Secretariat, 2013). The findings of this study may be useful in the zonation of the landscapes into the priority conservaiton areas and the multiple use areas while, formualting the landscape management plans.

Our study site falls in the Pamir trans-border region and reflects promising habitat with considerable presence of snow leopards. Since, snow leopards do not recognize international boundaries, we recommend undertaking similar surveys in the adjacent landscapes falling in China, Afghanistan, and Tajikistan to complete occupancy map across the whole range of Pamirs. This will help foster regional cooperation for conservation and rationalize establishment of transboundary protected area as envisaged in the international workshop held in Urumqi, China in 2006 (Moheb and Paley, 2016).

Conclusions

We conclude that the sign-based distribution surveys provide vigorous scientific knowledge with comparatively minimal resources involved and merit replication and further corroboration. Establishment of new and expanding on the existing PAs is highly suggestive together with transboundary cooperation in the Pamir to secure both national and range wide agenda of snow leopard conservation(Snow Leopard Working Secretariat, 2013). We also suggest to conduct similar surveys for wild and domestic prey to devise a holistic landscape conservation strategy that can be implemented far and wide.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.