Precise estimates of abundance and density are crucial for species conservation. For secretive felids, such as leopard cat Prionailurus bengalensis, acquiring such estimates based on conventional methods is difficult. We demonstrated the possibility of individual identification of leopard cats using coat patterns, and estimated their density using photographic capture-recaptures in a small watershed (182 km2) of Khangchendzonga Biosphere Reserve, Sikkim, India. On comparing the different body parts, we found that hind-quarter had the maximum usability (83.9%) for individual identification. The overall photo-capture rate representing an index of leopard cat relative abundance was calculated as 3.7 ± 1.27 captures/100 trap days. We used both non-spatial and spatially explicit capture- recapture (SECR) approaches to estimate leopard cat abundance and density. Our spatially explicit models estimated leopard cat density as 17 ± 5.33 (maximum-likelihood based approach) and 17.52 ± 5.52 (Bayesian approach with Markov chain Monte Carlo simulations) individuals/100 km2, while in our non-spatial model, density estimates varied from 18.01 to 22.25 individuals/100 km2. Camera trap results also indicated that the leopard cat used temperate and subtropical habitats to a large extent. Our study validated the applicability of camera trap based capture-recapture techniques to estimate the density of leopard cat. Therefore, we recommend the use of this technique with appropriate site-specific modifications for population estimation and monitoring of this species throughout its distribution range.
The leopard cat Prionailurus bengalensis is a small wild cat species distributed throughout Asia and adapted to a variety of habitats (Green 1991, Nowell & Jackson 1996, Sunquist & Sunquist 2002). Owing to its extensive distribution and variation in colouration and size, it was formerly considered to be several different species (Guggisberg 1975, Sunquist & Sunquist 2002). Later, based on rigorous morphological analysis, this concept got reconstructed and modified to 12 subspecies, widely differing in appearance (Groves 1997, Sanderson et al. 2008, Wilson & Mittermeier 2009). Contrary to other Asian small cats, the ecology of leopard cat has been extensively studied in terms of its behaviour, movement and spatial organisation, diet and prey selection and activity patterns (Inoue 1972, Rabinowitz 1990, Izawa et al. 1991, Grassman 2000, Rajaratnam 2000, Austin 2002, Khan 2004, Grassman et al. 2005, Rajaratnam et al. 2007, Watanabe 2009, Bashir et al. 2013). But, there is no standardised technique available for estimating abundance and density of leopard cat, and hence no reliable estimates from any part of its range are available.
Accurate and unbiased estimates of population size are desirable for species conservation and management (Silveira et al. 2003). Census tools must therefore be accurate, reliable, cost-effective and reasonably easy to apply (Jackson et al. 2006). However, for small elusive felids, such as leopard cat, such exercises are challenging. Estimates based on track observations are failure prone and unreliable, while radio-telemetry is costly and constrained to a small number of individuals (Karanth 1995, 1999). Karanth (1995) suggested an application of natural variation in fur markings for identifying secretive mammals. Camera trapping in combination with capture-recapture (CR) statistical modelling (non-spatial (Otis et al. 1978, Karanth & Nichols 1998) and spatially explicit capture-recapture (SECR; Efford 2004, Royle & Young 2008)) has been successfully used to reliably estimate densities for nocturnal, elusive felids with distinct coat patterns, such as tigers Panthera tigris (Karanth 1995, Karanth & Nichols 1998), leopards Panthera pardus (Harihar et al. 2009), jaguars P. onca (Kelly 2003, Silver et al. 2004, Sollmann et al. 2011), ocelots Leopardus pardalis (Trolle & Kery 2003), Geoffrey's cats L. geoffroyi (Cuéller et al. 2006), snow leopards P. uncia (Jackson et al. 2006) and Eurasian lynx Lynx lynx (Blanc et al. 2012, Weingarth et al. 2012). Cheyne & Macdonald (2011) used natural variation in fur markings for identifying leopard cat individuals, but could not estimate their density due to limited sample size. Now that SECR density estimation techniques have immerged as a robust tool to deal with low sample size, the leopard cat population density can be inferred. These models first determine an individual's activity centre by using the spatial location of captures and then estimate the density of these activity centres across a precisely defined polygon containing the trap array (Gardner et al. 2009, Royle et al. 2009b), thereby avoiding the issue of estimating the effective area sampled (Sollmann et al. 2011). In addition, SECR methods also deal with uncertain edge effects and spatially heterogeneous detection probability caused by movement in conventional animal trapping (Efford 2004, Borchers & Efford 2008).
In our study, we aimed at developing a protocol for identifying leopard cats from their coat patterns and estimating their abundance and density. This ecological information on leopard cat from high- altitude eastern Himalayan landscape is rare and hence generates a basis for its future conservation and management.
Material and methods
Study area
Our study was carried out in the Prek ‘chu’ catchment of the 2,620 km2 Khangchendzonga Biosphere Reserve (National Park (core zone) = 1,784 km2, buffer zone = 836 km2), located in a small state of Sikkim in India (27°30‘-27°55‘N, 88°02‘-88°37‘E; Fig. 1) which is positioned at the convergence of Palaearctic, Africo-tropical and Indo-Malayan biogeographic realms (Mani 1974) and also included in the eastern Himalayan global biodiversity hotspot (Myers et al. 2000). The Khangchendzonga Biosphere Reserve (BR) encompasses a sharp elevation gradient of 1,220 to 8,586 m a.s.l. varying within an aerial distance of just 42 km with about 90% area above 3,000 m a.s.l. and 70% area above 4,000 m a.s.l. (Tambe 2007). The area of Khangchendzonga BR is divided into seven catchments or river subsystems: Lhonak, Jemu, Lachen, Rangyong, Rangit, Prek and Churong. Among these, Prek catchment (27°21‘-27°37‘N, 88°12‘-88°17‘E; 182 km2) was selected for intensive studies because it represents all the elevation, habitat, slope and aspect classes of the BR in similar proportions (Sathyakumar et al. 2011). Our study was conducted for a period of about two years from April 2008 to December 2009. Due to the topography and remoteness of the area, all field activities were carried out in the form of expeditions, i.e. trekking, camping and sampling different areas of catchments.
Reconnaissance surveys
Reconnaissance surveys were done from April 2008 to February 2009 in different catchments of Khangchendzonga BR (Sathyakumar et al. 2009) to generate baseline information on the presence and distribution of leopard cats based on trail sampling, sign surveys and local interviews. This information was also used to identify the most adequate site for intensive camera trapping.
Survey design
For estimating the population size of leopard cat in the Prek catchment, grid-based sampling was done for a short period of two months (October-November 2009) to assure no significant changes in the population size to be sampled and hence maintain geographical and demographic closure (Karanth 1995). The entire subtropical, temperate and lower part of subalpine zone of the Prek catchment (< 3,500 m a.s.l. elevation; considering the species' highest recorded occurrence in literature and non-detection above 2,700 m a.s.l. during reconnaissance surveys) was divided into 2 × 2 km grids. Within all accessible grids, 23 infra-red triggered camera traps (18 Stealthcam, LLC, Grand Prairie, Texas, USA and five Moultrie, Moultrie Feeders, Alabaster, Alabama, USA) were deployed along natural trails and junctions (Fig. 2). Accessibility of the grids being a limiting factor, camera traps were deployed to warrant minimum possible inter-trap distance based on the average home-range size of the species (3-14 km2; Grassman et al. 2005, Rajaratnam et al. 2007). Such coherence between the home-range size and the camera trap deployment design is a prerequisite for non-spatial CR models (Karanth & Nichols 1998), but not for SECR models (Sollmann et al. 2012). Camera traps were placed at 15-25 cm above ground, attached to a rock or tree trunk at 3-5 m from a trail or point where animal movement might be expected. Cameras were set with one minute delay between successive activations, were operational for 24 hour-monitoring during the entire sampling session of 60 days and were checked at 10-15 day intervals.
Identification of leopard cats from spot patterns
We evaluated the reliability of identifying leopard cats from their spot patterns following a double-blind observer identification protocol designed for pumas Puma concolor by Kelly et al (2008) and later applied to striped hyenas Hyaena hyaena by Harihar et al. 2010. To implement this protocol, the photographs were numbered and three sets were distributed among three investigators. Each investigator independently analysed each photograph. Based on the investigator's confidence in identifying each body part (head, neck, fore-, mid-, hind-quarters and tail) in a photograph, and by comparing its characteristics (dots, stripes, tail rings and other marks) with that of other photographs (Fig. 3), the three investigators independently assigned identities to each photograph and categorised them into different individuals. All photographs hence assigned to an individual were grouped and represented its capture history. Based on this procedure, each investigator ended up with a data set (capture histories of individuals identified by each investigator). In addition, each investigator had kept a record (for each photograph) about the body parts used to identify a photograph before assigning identity. These records were compared and pooled across investigators for all photographs. A Kruskal-Wallis test was performed to test the differences in the use of various body parts by the investigators to assign individual identities to the photographs.
Statistical analysis
Photo-capture rate
Photo-capture rate as an index of relative abundance was calculated for the entire trapping session and for each habitat type. The photo-capture rate was calculated as the number of photographs divided by the number of effective trap days per site (Carbone et al. 2001) and represented as relative abundance index (RAI) per 100 trap days (O'Brien et al. 2003).
Non-spatial density estimation
The entire sampling session was divided into 20 occasions of three days each and tested for closure assumption using Close-Test (Stanley & Burnham 1999). Non-spatial estimation of abundance (N) and capture probability (p) was done using full closed CR model in program MARK version 6.1 (White & Burnham 2000). MARK offers four main population size estimators that differ in their assumptions about capture probability as constant (Mo) or varying as a function of time (Mt), behaviour (Mb) and individual heterogeneity (Mh). Density estimation was done by dividing the abundance estimate (N) by the effective sampled area. For calculating the effective sampled area, a minimum convex polygon (MCP = 62.28 km2) was created over the trapping locations and a buffer width equal to half mean maximum distance moved (½ MMDM) by a recaptured individual was added over the trap polygon (Karanth & Nichols 1998). We then subtracted the unsuitable habitat (beyond subalpine, i.e. > 3,500 m a.s.l.) covered under the buffer and MCP using the spatial analyst tool of Arc GIS 9.3, and used only area covered by suitable habitat to obtain the effective area sampled and subsequently density.
Spatially explicit density estimation
Considering the fact that space and movement have no explicit manifestation in classical MMDM based CR models (Royle et al. 2009a), we tried to address this issue through SECR models. A usual framework for developing spatial models is based on point process models (Efford 2004, Borchers & Efford 2008, Royle & Young 2008), which assume that each individual i in the population has a fixed point associated with it considered as its centre of activity si (two-dimensional coordinate representing a point in space about which its movements are concentrated). SECR models hence operate as generalised linear models with random effects (i.e. GLMMs) with many unknown variables (random effects) in the model (Royle et al. 2009a, Sollmann et al. 2011) particularly the activity centres of each individual (si) and the number of such activity centres (i.e. the population size N). We used two different calculation techniques to estimate SECR density of leopard cat. The maximum likelihood-based method using program DENSITY 4.4.1.2 (Efford 2007) which directly estimates density by fitting spatial detection functions to CR data from arrays of passive detectors such as camera traps (Efford 2004). Here, the probability density functions for detections of animals based on distance from activity centres are modelled using hazard rate, half-normal or exponential detection functions (Efford et al. 2008). These models can be considered as mixture models where the mixture is over the distribution of animal locations, and the estimators are based on the marginal distribution acquired by integrating the joint likelihood over the distribution of unobserved locations (Borchers & Efford 2008). Maximum likelihood-density (MLDens) was hence estimated using the estimator that best explained the individual capture probability in MARK and based on minimum model AIC value. In order to estimate the unknown parameters (including the number of individuals (N) and their activity centres (si) and hence the density), we also used a Bayesian framework based on data augmentation (Royle et al. 2007) which can provide valid inferences even with small sample sizes (Sollmann et al. 2011). We created a 5-km large buffer around the MCP of the trap array to ensure inclusion of all individual home ranges within a reach of cameras. These were described by 1,492 equally spaced pixels, each representing an area of 0.2025 km2 within a total area of 302.94 km2. The area judged as non-suitable habitat (beyond subalpine, i.e. > 3,500 m a.s.l.) was excluded from MCP and 5-km buffer using the spatial analyst tool of Arc GIS 9.3, and the remaining suitable habitat was hence represented by only 902 pixels. The activity centres were therefore assumed to be uniformly distributed over this discrete space of 902 pixels, an area S of approximately 182.665 km2. To implement data augmentation, we supplemented the n observed encounter histories with some large number of ‘all-zero’ histories say M-n (Royle et al. 2007, Sollmann et al. 2011). We assumed that this M includes the actual N (population size) as a subset and hence choose M sufficiently large (larger than the largest possible population size in the area) so that the posterior of N is not truncated (Royle et al. 2009a). “This can be achieved by trial and error with no philosophical or practical consequence” (Royle & Young 2008, Royle et al. 2009a). This reformulation of the model, based on data augmentation, is a zero-inflated binomial mixture where the number of activity centres N in area S are estimated as a fraction of M (Sollmann et al. 2011). We implemented a Markov chain Monte Carlo (MCMC) simulation method which estimates the joint posterior distribution of the unknown quantities in a statistical model in R version 2.14.1 (R Development Core Team 2009) with R package SPACECAP 1.0.3 beta (Gopalaswamy et al. 2011). The software uses an uninformative prior to estimate the posterior distribution which is recommended, and considered to be appropriate, while dealing with small data sets and in cases with no reliable information on the priors (Royle et al. 2009a). MCMC chains are started at arbitrary parameter values and are needed to converge to the stationary distribution (single peak) within an acceptable error. Since the successive iterations depend on the outcome of the preceding iteration, the effect of the beginning value may be reflected in a number of initial iterations which hence need to be discarded (the burn-in). This characteristic can also result in autocorrelation of successive iterations which can be reduced by specifying a thinning rate as every ith iteration used in the characterisation of the posterior distribution of the parameters (Sollmann et al. 2011). Hence, for model analysis using an uninformative prior, we ran MCMC chains (trial and error) with half-normal and negative exponential detection functions by setting different iterations, burn-in and data augmentation values. The chain finally got stabilised at 50,000 iterations, a burn-in of 5,000 and a thinning rate of 10 and data augmentation value 10 times the individuals with known capture histories using a half-normal detection function in a Bernoulli encounter model to characterise the posterior distributions.
Results
Among the total camera-trap samplings of 1,380 days, we used 1,109 operational trap days for analysis due to temporary malfunctioning of cameras. A total of 45 photo-captures (31 right and 14 left flanks) were recorded, out of which 27 and 10, respectively, were categorised as usable by the investigators. The numbers of leopard cat individuals identified by investigators ranged from 13-14 using right flank and five using left flank photo-captures. Therefore, only right flank photo-captures were used for further CR analysis. This differential identification of individuals by investigators resulted in three different data sets, as mentioned above, which were analysed separately. On comparing the different body parts, we found that hind-quarter was utilised more successfully in classifying 83.9% of photographs, followed by mid-quarter (66.7%) and tail (64.2%; Table 1).
The overall RAI/100 trap days for leopard cat was 3.72 ± 1.27 for the entire trapping session, while habitat-wise RAI was 3.55 ± 1.91 for subtropical (elevation < 2,000 m a.s.l.) and 5.81 ± 2.30 for temperate habitat (2,000-3,000 m a.s.l. elevation). No photo-captures were recorded beyond 2,700 m a.s.l. elevation. The results of the closure test did not provide enough support to indicate violation in the closure assumption for all three data sets. The non-spatial ½ MMDM based model estimated leopard cat densities varying from 18.01 to 22.25 individuals/100 km2 under the best selected Null (Mo) model estimator with no significant difference across estimates evident from overlapping standard errors between investigators (Table 2). For further SECR analysis, the data set with 13 individuals using 27 captures was used based on its higher capture probability (0.094). The maximum likelihood SECR model, selected based on minimum AICc value, described the detection function with a hazard rate function. The density estimate was 17 ± 5.33 individuals/100 km2. As regards the SECR model in a Bayesian framework, the density estimate reached 17.52 ± 5.52 individuals/100 km2 while using a half-normal function for the detection function. Additionally, MCMC chains using negative exponential detection function poorly converged/stabilised resulting in less reliable estimates (Table 3). Posterior pixel densities for each activity centre of 0.2025 km2, estimated through point process of the Bayesian model, ranged from 0.002 to 0.117 individuals (Fig. 4). Scaling these figures to 100 km2 yielded a density range between 1.097 and 58.163 leopard cats/100 km2.
Discussion
Our study demonstrates the possibility of identifying individual leopard cats from their spot patterns, and established that hind-quarters contain maximum information and, thus, should be given priority for identification in future camera-trap studies. Robust estimates of abundance and density of leopard cats were not available from any part of the species' range prior to our study. However, during the same study period, another study carried out in Sabah, Malaysia, reported density of leopard cats ranging from 9.6 to 16.5 individuals/100 km2 based on camera trapping (Mohamed et al 2013). Our results based on photo-capture index reestablishes leopard cat as the most abundant felid in the temperate habitats of the Khangchendzonga BR (Sathyakumar et al. 2011). Comparisons of our estimated RAI of leopard cat (3.72 ± 1.27) with different studies across its range suggested that leopard cats were more abundant in Khangchendzonga BR compared to secondary forests of Peninsular Malaysia (Azlan & Sharma 2006; RAI = 1.44) and Sabangau peat-swamp forest of Indonesian Borneo (Cheyne & Macdonald 2011; RAI = 2.45). Photo-capture rate being also considered as a relative index of animal's spatial use (Carbone et al. 2001), our results indicate that leopard cats used more temperate and subtropical habitat than subalpine and alpine habitats.
Our results showed that the densities estimated under different models (non-spatial and spatial) were almost similar. Generally, non-spatial ad hoc models are known to overestimate densities as these underestimate individual movements (Royle & Young 2008). The possible explanation for such resemblance in density estimates under different models may be due to differences in movement patterns of individual leopard cats, i.e. restricted for some individuals, whereas very wide-ranging for others probably engaged in establishing territories. Estimates of the MMDM are constrained by the size of the sampling grid, as camera traps do not capture any movements beyond it (Sollmann et al. 2011), and hence the estimates are less reliable and cannot be extrapolated for a larger area. On the contrary, SECR models are more comprehensive as these address the movement patterns more explicitly which was even evident from comparatively high standard errors in their density estimates. While comparing among SECR models, although classical maximum-likelihood based inference procedures are asymptotic and establish unbiasedness for large sample sizes (Efford 2011), their relevance to small sample situations has been inferred questionable (Royle & Young 2008). Conversely, Bayesian inferences do not rely on asymptotic arguments and are valid, regardless of the sample size (Royle et al. 2009a). But, in reference to identification of some bias and poor coverage of credible intervals for Bayesian estimates of density attributed in part to small samples by Marques et al (2011) and Efford (2011), the problem continues to be a topic of debate. However, in our study, the Bayesian model validated the appropriateness of using a half-normal detection function in explaining the posterior distribution in a real ecological sense. The predicted posterior pixel densities by the Bayesian method indicated that high‐density areas were situated in dense temperate and subtropical forests of low elevation. It corroborates the results of relative index of animal's spatial use deduced from photo-capture rates mentioned above.
Our study also demonstrates the protocol for identifying leopard cat individuals from their coat patterns, presents photographic CR sampling as a promising approach in this context, and finally elucidates the applicability of SECR models to reliably estimate its abundance and density, even with small sample sizes. In order to use this technique for leopard cat, sign surveys are a prerequisite for identification of sites for camera trap deployment. Future studies are required to optimise camera placement for leopard cats in terms of distance from the path, camera orientation with respect to path in case of narrow trails (in mountainous terrain) and height of deployment in order to overcome the sample size constraints of limited captures and poor image quality. For instance, in a steep and narrow trail, it would be better to deploy the camera in such a manner that the field of view is maximised by orienting the camera unit to a suitable angle with reference to the trail. We recommend this technique for estimating leopard cat densities and for long-term monitoring. In order to improve and optimise this monitoring protocol, we suggest that intensive sampling using a smaller grid size (1 × 1 km) in a double flank mode and using lures would increase capture probability, and consequently reduce confidence limits in the estimates. Such a precise estimate would be extremely useful to deduct changes in status of leopard cat in an area. Further, this technique would also help in estimating other population/demographic parameters such as population trend, vital rates, recruitment and survivorship which are crucial for management planning. Based on the population trends observed through regular monitoring, reappraisal of the status of leopard cat could be done and appropriate conservation strategies could be develop throughout the species distribution range.
Acknowledgements
We are grateful to the Department of Forests, Environment and Wildlife Management, Government of Sikkim for granting us permission to work in Sikkim. We thank the Wildlife Institute of India, Dehradun for providing us the grants and support. We thank Mr. Sukhbahadur and Mr. Sukhraj for their help during field work. We thank the two reviewers for their valuable comments that helped in improving this manuscript.
References
Figure 1.
Location of Khangchendzonga Biosphere Reserve in Sikkim showing the different watersheds including Prek ‘chu’, the intensive study area.
Figure 3.
Example of identification of separate leopard cat individuals based on pelage and spot patterns. Each photograph was sectioned into head, neck, fore-, mid-, hind-quarters and tail to compare quality of body parts (separated by broken lines) for individual identification. Photos A), B) and C) represent different individuals, whereas C) and D) represent the same individual as indicated by similar patterns.
Figure 4.
Map of posterior density of leopard cat activity centres in 450 × 450 m pixel area. Gray scale indicates different density levels, dark indicates higher densities and light corresponds to lower densities. Plus signs denotes cameras where leopard cats were photo-captured.
Table 1.
Percentage of exact matches among the three investigators based on photo-captures with respective usability for each body part and significance of difference in usability categorisation of each part across investigators.
Table 2.
Abundance and density estimates for leopard cat using conventional (non-spatial) MMDM based capture-recapture method under the best selected Null (Mo) estimator, based on the number of individuals identified by different investigators in Prek ‘chu’ catchment of Khangchendzonga BR. M(t+1) = Number of individuals captured, C+R = Number of captures and recaptures, χ2 = chi square statistics, P = significance level, ½ MMDM = Half mean maximum distance moved (km), ETA = Effective trapping area (km2), P hat = Capture probability, N = Population size, SE = Standard error, D = Density of individuals/100 km2.
Table 3.
Abundance and density estimates for leopard cat using SECR method with maximum-likelihood (A) and Bayesian (B) approach in Prek ‘chu’ catchment of Khangchendzonga BR. AICc = Akaike Information Criteria for small sample sizes, SE = Standard error, Psi = Data augmentation parameter, σ = Spatial scale parameter, g0 = Detection probability (frequentist), λ0 = Expected encounter frequency at trap location considered as home-range centre, N super = Population size of individuals having their activity centres within the effective trapping area with 5-km buffer (i.e. 182.665 km2), D = Density of individuals/100 km2, CI = Confidence interval, LL = Lower limit, UL = Upper limit.
