Estimating snow leopard population density using camera traps

We deployed Reconyx RM45 camera traps at 30 sites over an area of 953 km² (Minimum Convex Polygon joining the outermost trap locations) with an average inter-trap distance of 4035 m (SE = 374m) (Fig 1). The camera traps were deployed from October 2011 to January 2012 for 80 days with an overall trap density of 3 camera traps per 100 km² following recommendations of placing at least two traps per average home range [17] or at least two traps per average female home range [18]. The camera traps were deployed at sites where we encountered relatively high frequency of snow leopard signs such as scrapes, pugmarks, scats and scent marks, especially around terrain features that snow leopards are known to prefer for marking and movement such as ridgelines, cliffs and gully beds. We used a combination of single (n = 14) and double (n = 16) camera trap placement to optimise coverage and identification of individuals (Fig 1A). Double side camera traps were installed to enable capture of both side flanks of as many snow leopards as possible so they can be used to improve our ability to identify individuals with only single flank captures. Our cameras recorded snow leopards at 25 out of 30 sites without using any baits or scent lures.

Individual snow leopards captured in the images were identified based on their pelage patterns by two independent observers using at least three similarities or differences [5, 19]. There was no discrepancy in the identified individuals reported by the two observers. We only count each set of photographs as a new encounter when it was separated from another set of photographs from the same snow leopard by at least four hours. This was done to prevent the overdispersal of counts in using count detectors for our analysis. The mean time between consecutive encounters of the same animal on the same camera trap was 537 hours (95% CI: 409–665 hours) that ensures the validity of the count detector. We obtained a total of 2,830 snow leopard images from 124 encounters. A total of twelve encounters were discarded as the pictures were not good enough to identify the individuals. Using a mix of both side and right side only flanks, we obtained complete identification of 16 individual adult snow leopards. We assumed individuals moving about on their own (dispersed from their mother) to be mature individuals. There was no discrepancy between the two observers in individual identification of snow leopards. Following concerns raised by Johansson et al [20], we used the Snow Leopard Identification: Training and Evaluation Toolkit (https://camtraining.globalsnowleopard.org/leppe/login/) to test the skills of both of our observers in identifying snow leopards. Our observers scored 96.3% and 88.9% accuracy respectively in identifying snow leopards from 40 blind, independent trials, thus leaving us confident of identifying individuals with reasonable accuracy. Snow leopard capture histories were built using the standard count detector format of the `secr’ package in R [21] where each encounter of an identified cat was linked to a detector (camera trap site), whose location, period of operation and other relevant covariates were recorded in a separate table. We restricted the study period to 80 days and assumed that the population was closed and that there was no temporal effect on detection probability of snow leopards during the sampling period.

Typically, SCR models assume that the expected encounter rate depends on the Euclidean distance between detector and activity centre. This may not always be true in highly structured environments such as steep mountains. For example, we may record more encounters for a snow leopard in a distant trap than a closer trap if the habitat between the closer trap and activity centre has a higher resistance to movement (e.g. a deep gorge separating two detector locations). Royle et al. [22] and Sutherland et al. [23] proposed replacing Euclidian distance with a least-cost path distance (ecological distance) in which movement cost depends on the habitat. The method involves the estimation of movement cost parameter(s) simultaneously with other SCR parameters. Sutherland et al. [23] demonstrated that violations of the Euclidean distance assumption could bias estimates of density, and they suggest that least-cost distance be tested in highly structured landscapes.

We used the maximum likelihood-based SCR models [24] to estimate density while investigating the effect of least-cost path distance on movement, using package ‘secr’ [21] in R [25]. The method involves integration over a 2-dimensional region containing the possible (and unknown) locations of the activity centres of animals at risk of detection. The region of integration is based on a polygon extending a certain distance (the buffer width) beyond the outermost traps.

We used the inbuilt ‘suggest buffer’ function of ‘secr’ to arrive at a buffer of 24,000 meters assuming it to be wide enough to keep any bias in estimated densities as acceptably small (i.e. snow leopards with activity centres beyond 24 km from the outermost traps had a negligible probability of being captured in the detectors). Areas above 5,200 m from mean sea level were excluded from the set of possible activity center locations because they are mostly devoid of vegetation and prey species. We defined an integration area with a spacing of 500 m, resulting in 20,513 pixels for the entire integration area. We used the model with minimum AIC to estimate population size (N) and density (D) over the integration region [26], but use a model averaged density surface to present the distribution of the density of snow leopard activity centres.

Spatial capture-recapture models The spatial distribution model in SCR is a spatial Poisson process for animal activity centres whose intensity (expected number of animal activity centres per unit area) can be homogeneous (constant over space) or inhomogeneous (varying over space) [24]. We use the notation D(x;θ) for density, signifying that density is a function of activity centre location, x, which is a vector representing the x and y coordinates of an activity centre, and of parameters represented by the vector θ.

We fitted SCR models with various combinations of covariates defined a priori. A candidate model set was developed to investigate the effects of various covariates potentially influencing snow leopard behaviour, ecology and natural history. We investigated models with various combinations of covariates for the density model, the encounter function intercept model, and the encounter function range model. The general forms of the density model, and encounter function intercept and range models, respectively, are as follows:

where

x_d(s) is the dth spatially referenced covariate at location s that affects density (D), and θ₀ and θ_d are the density intercept parameter and dth regression parameter;

x_l is the lth covariate that affects expected encounter rate at distance zero (λ₀), and ϕ₀ and ϕ_l are the intercept parameter and lth regression parameter for expected encounter rate at distance zero;

x_i is the ith covariate that affects the encounter rate range parameter (σ), and β₀ and β_i are the range intercept parameter and ith regression parameter.

Half-normal encounter function forms were used, such that the expected number of encounters of an animal at a camera that is a distance d from its activity centre is E(n) = λ₀exp{−d²2σ²}.

For snow leopard density, we considered models in which the x_d(s)s were wild prey density, large and small livestock density, overall livestock biomass, least cost distance from settlement, terrain ruggedness and altitude at s. We investigated the effect of topographies (a factor with levels “ridgeline”, “cliff” or “gully bed”), and the effect of single and double camera traps on the encounter function intercept and range parameters. We also investigated models in which movement cost depends on altitude by replacing Euclidian distance with a least-cost path distance in which movement cost depended on altitude.

We modelled D(s) as a function of six spatial covariates (x_d(s)s) that could affect snow leopard density (Fig 1). These included terrain ruggedness (typical snow leopard habitats are steep and rugged) [27], altitude (snow leopard densities are known to be a function of altitude) [28], wild prey density (believed to be the main determinant of snow leopard population abundance) [29], stocking density of large-bodied livestock and small-bodied livestock (potential prey for snow leopards, source of disturbance, and competitors for wild prey) [27–30]. The terrain ruggedness across the integration region ranged from 12.63 to 74.55 (Mean = 40.97, SD = 9.45), the altitude ranged from 3298 meters to 5500 meters above mean sea level (Mean = 4775, SD = 509), the least cost distance (after estimating the associated parameter) from nearest village, using elevation as a cost function from 1.21 to 1580.82 (Mean = 599.81, SD = 415.62). The density of large livestock ranged from 0 to 14.11 per km² (Mean = 1.41, SD = 2.21) while that of the small bodied livestock ranged between 0 to 11.39 (Mean = 0.89, SD = 1.75) per km². The mean livestock biomass ranged between 0 and 5,539 kg (Mean = 230.71, SD = 471.06).

Terrain ruggedness was derived using the terrain ruggedness index [31] from a 30 × 30 meter Digital Elevation Model from Aster Global Digital Elevation Model data using the terrain analysis plugin in the Quantum GIS 3.16.3 software [32]. Livestock density was determined through a door to door censuses in 51 villages in the integration region. The pastures used by each village were mapped using Google Earth and livestock stocking densities for small- and large-bodied livestock were computed separately (as they are often herded separately [33] by dividing the total livestock heads using a pasture, by the area of the pasture in square kilometres. We used the average biomass of large-bodied and small-bodied livestock [34] to estimate the livestock biomass availability to snow leopards across the integration region. We smoothened the livestock density and ruggedness surfaces across the integration region by averaging over a moving window of 1.5km.

The abundance of wild prey, which primarily included blue sheep and ibex, was estimated using the double observer survey technique for the entire integration region [35] (Table 1) between April and June 2012. Four teams, comprising two observers each, carried out the surveys for eight days to cover the entire study area. Observers recorded the GPS coordinates of the sightings, the group size and age-sex classification of the groups encountered. The unique identity of each observed ungulate group was established through immediate post survey discussions between two observers using the age-sex classification, size and the location information of sightings [35]. The integration region was divided into seven blocks delineated based on natural topographic features in the landscape such as rivers and contours of prominent ridgelines. For each block, the cumulative number of wild ungulates encountered by the two observers were calculated. The relative density of wild ungulates for each block was estimated as total number of wild prey in a block divided by the size of the survey block (Table 2). The wild prey density surface was smoothened by averaging over a moving window of 5 km.

C is the number of groups seen in both surveys; S₁ is the number of groups seen in first survey only; S₂ is the number of groups seen in second survey only; G^ is the estimated number of groups; N is the naïve population estimate; N^ is the estimated population size; P₁ and P₂ are the means of the estimated detection probability for observers one and two, respectively.

LCL is the 95% interval lower limit, UCL is the 95% confidence interval upper limit.

We developed an a priori set of models that we anticipated would best explain the variation in the density of snow leopards. Our global (most complex) model included terrain ruggedness, linear and quadratic effect of altitude above mean sea level, density of wild prey, stocking density of small-bodied livestock, stocking density of large-bodied livestock, least-cost distance from settlements considering altitudinal gain as the added cost, and cumulative livestock biomass. We then fitted 20 candidate sub-models using subsets of the variables used in the global model. Each candidate sub-model represented a specific hypothesis about the relationship between snow leopard density and how snow leopards use space about their activity centres, and explanatory variables. We used Akaike’s Information Criterion (AIC) for model selection [36]. All data analysis was implemented using package secr [21] in program R [25].