Density estimation of tiger and leopard using spatially explicit capture–recapture framework

Study area

Bandhavgarh Tiger Reserve (BTR) is located between 23° 27′ 00″ to 23° 59′ 50″North latitude and 80° 47′ 75″ to 81° 15′ 45″ East longitude in Central India. BTR consists of two Protected areas: Bandhavgarh National Park (BNP) in the south and the Panpatha Wildlife Sanctuary (PWS) in the north (Fig. 1). BNP and PWS constitute the core zone of the reserve, having a combined area of 716 km². The surrounding buffer zone has an area of 820 km², adding the total area of the reserve to 1,536 km². A more detailed account of the study area is available in Rather, Kumar & Khan (2020). The vegetation comprises moist peninsular low-level sal forest, northern dry mixed deciduous forest, dry deciduous scrub, dry grassland, and west Gangetic moist mixed deciduous forest (Champion & Seth, 1968). Sal and sal mixed forests occur in the major portion of the reserve. Besides tiger and leopard, sloth bear (Melursus ursinus), Indian wolf (Canis lupus), Asiatic wild dogs (Cuon alpinus), and striped hyaena (Hyaena hyaena) are notable carnivore species occurring within the reserve. Major prey species include chital (Axis axis), sambar (Rusa unicolor), Indian guar (Bos gaurus), barking deer (Muntiacus munjtak), Indian gazelle (Gazella bennetti), four-horned antelope (Tetracerus quadricornis), and Indian blue bull (Boselaphus tragocamelus). Most of the reserve has a flat to the gentle slope with an average elevation of 570 m (asl). The annual average rainfall in reserve is reported to be 1,173 mm, most of which occur in the monsoon season. The reserve has a fair road network presence, allowing easy and adequate camera traps placement along roads and trails.

Camera trap survey

We conducted extensive reconnaissance within the study area to identify the most optimal locations for setting up camera traps to maximize the probability of detecting as many individuals and obtain as many photo-captures of the individual tigers and leopards as possible. We selected trap location for camera traps based on secondary evidence such as pugmarks, scats, scratches, marking on trees, etc. The average inter-trap distance was 856 meters to avoid leaving any large gaps between the successive camera traps (Supplemental Information S1). At any time, 20–25 pairs of camera traps Cuddeback^™ model C1 (total 35 pairs) remained active within the study area yielding 2,211 trap nights at 220 sites. Camera traps were placed along roads and trails at a distance of 5–10 m from the center of roads or trails. At each location, a pair of camera trap was tied to the trees on either side at the height of 30–40 cm from the ground. The delay between the successive photo-captures was set to 5 seconds to increase the chances of photo-capturing animals traveling either in groups or mothers accompanying cubs. Camera traps were kept active 24 h at the stations where the theft risk was minimal and checked daily. No baits were used to lure the animals towards camera traps. The study area was divided into three blocks corresponding to the three buffer ranges of the reserve. Camera traps were placed within each sampling block for eight consecutive days before being moved to new sites. Camera traps were set for two months in each sampling blocks between 2016 and 2017.

Habitat covariates

We derived a set of eight habitat covariates to model the density of the tiger and leopard (Supplemental Information S2). The habitat variables were based on a recent study predicting tiger and leopard occurrence in the study area (Rather, Kumar & Khan, 2020). We used land use land cover map of the study area prepared by the Indian Institute of Remote Sensing (IIRS; www.iirs.gov.in) and derived eight habitat variables at the spatial resolution of 1km. We used a moving window analysis in ArcGIS (version 10.3) and extracted the habitat variables for each camera trap location. We used a circular buffer of 1,000 m radius around each camera trap location and calculated the percentage of available habitat types. We used road and river density calculated within the radius of 1,000 m instead of distance to roads and rivers due to the high concentration of roads within the study area. The density of roads, percentage of human settlements, and percentage of degraded forests (mosaic of croplands and natural vegetation) within the spatial extent of 1,000 m were considered as a proxy to human disturbance.

Individual identification of tigers and leopards

Each photograph was carefully examined for the shape and pattern of natural markings; stripes in tigers and rosettes in leopards using both right and left flanks, limbs, tails, and forequarters (Schaller, 1967; McDougal, 1977; Karanth, 1995; Franklin et al., 1999). Tigers and leopards were individually identified and given a unique identity. Subsequent identification of new individuals was achieved by comparing and matching each photograph with the known and uniquely identified tigers and leopards. Respective capture histories of identified tigers and leopards and their trap locations were constructed in the form of two separate text delimited files (Supplemental Informations S3, S4 and S5). The analysis was done using the package ‘secr’ version 4.3.1 (Efford, 2020) implemented in R (R Core Team, 2018).

Buffer width and detection functions

Since there was a reasonably suitable habitat around camera traps, so it was reasonable to pay attention to the area immediately around the camera traps by specifying a habitat buffer. In SECR, the buffer width is not critical as long as it is wide enough so that the animals occupying edges have virtually zero chances of appearing in the sample. The buffer widths were based on the average inter-trap distances (initial sigma) moved by animals. In SECR, the buffer widths are multiplies of initial sigma, and the fit is assessed when the density estimates do not change with an increase in buffer widths. In this analysis, the average inter-trap distance (initial sigma) moved by tigers and leopards was 1,581 and 1,574 m. The buffer widths were calculated by multiplying the initial sigma values for each detection function until no change in density was observed.

Density models

We fitted both the null models (D ~ 1, lambda0 ~ 1, sigma ~ 1) and the models with habitat covariates (D ~ cov, lambda0 ~ 1, sigma ~ 1) to model the distribution of animal activity centres within the defined state-space S. Among null models we tested three models each for tiger and leopard as: (D ~ 1, lambda0 ~ 1, sigma ~ 1), (D ~ 1, lambda0 ~ b, sigma ~ 1), (D ~ 1, lambda0 ~ bk, sigma ~ 1). We tested different hypotheses by studying the effects of habitat covariates on the tiger and leopard density. First, we defined a null model (D ~ 1, lambda0 ~ 1, sigma ~ 1). Subsequently, we assessed the effect of the following covariates on the density of tiger and leopard: (i) sal dominated forests with thick vegetation cover, (ii) moist deciduous forests, (iii) dry deciduous forests, (iv) degraded forests, (v) road density, (vi) river density, (vii) human settlements and (viii) terrain ruggedness. We expected tiger and leopard density to increase with sal-dominated forests, moist deciduous forests, and dry deciduous forests as these forest types are associated with high prey presence. We also expected density to be positively correlated with rugged terrain due to the low anthropogenic influence in rugged topographical areas. We also expected tiger and leopard densities to be negatively correlated with disturbed habitats (degraded forests). We set a maximum buffer of 8 km based on initial sigma around each trap site based on rigged density estimates. We calculated the Akaike Information Criterion (AIC) for each competing model and selected the best model based on the difference between AIC scores (ΔAIC).

Population estimates

We obtained 95 detections of 33 individual tigers and 74 detections of 38 individual leopards over 2,211 trap nights in three sampling blocks. Tigers and leopards were detected in 87 and 70 camera trap stations, respectively. The overall capture probability (M (t+1)/n) of leopards (0.56) was slightly higher than tigers (0.45). The cumulative number of detected individuals (tigers and leopards) on each occasion (t) was pooled together for all sampling blocks to estimate sampling adequacy (Supplemental Information S6). The camera trap arrangement for eight consecutive days at each site before moving to new sites used in this study was adequate for sampling tigers. The unique individuals captured did not reach asymptote as promptly as tigers in the case of leopards (Supplemental Information S6). Sharma et al. (2010) found that camera traps failed to detect all known tigers in Kanha Tiger Reserve, even after prolonged sampling. Thus it is challenging to achieve total counts in natural populations using camera traps even after intensive sampling efforts (Sharma et al., 2010). We achieved a sampling coverage of over 500 km² in the reserve buffer zone (836 km²). Thus the trapping area met the size requirement as recommended by Silver et al. (2004), Maffei & Noss (2008) and Sharma et al. (2010) for estimating the densities of tigers and leopards.

About 40% of trap-revealed movements in tigers were within the range of 1,000–2,000 m, and 25% of trap-revealed movements in leopards were within the range of 3,000 m; thus, there is a high probability that peripheral individuals had a good chance of being photo-captured even if their home ranges were centered outside the plotted area.

Fitting SECR models and density estimations of tiger and leopard

The models were fitted using the buffer W of 7,908 and 7,870 meters for HN, EX and HR detection functions under null models to estimate the densities of tigers and leopards, respectively. We estimated densities using three model specifications: a null and permanent global learned response model (b) and a permanent detector-specific learned response (bk) model. The Null model assumes that detection of all individuals at all camera traps is governed by the same detection versus distance curve and is represented as:

(g0 ~1, sigma ~ 1)

The learned response model or behavioral response model (b) assumes that detectors (camera traps) incite a behavioral response in individuals after their first encounter with detectors. Thus the probability of capturing an individual on any later occasion is affected (Otis et al., 1978). Moreover, the learned response may be specific to the detector location (bk) rather than generally applying across all detectors. We considered the possibility of either sort of induced behavioral response in tigers and leopards.

The density of the tiger and leopard under the null model was estimated to be 3.28 ± 0.65 (Mean ± SE) per 100 km² and 4.04 ± 0.83 (Mean ± SE) (Table 1). Formal model comparisons based on AIC placed the detection functions with longer tails (EX) ahead of (HN) in tiger (Table 2). However, the density estimates from the HN detection function reached a plateau fairly promptly with increasing buffer W (Supplemental Information S7) compared to negative exponential (EX) and hazard rate (HR). The AIC differences between the models using different detection functions (null models) in the tiger were also very small (Table 2). The density estimates from the null model (D~1 g0~1 sigma~1) and behavioral response models (D~1 g0~b sigma~1; D~1 g0~bk sigma~1) were not noticeably larger than each other (Table 1). The null models (D~1 g0~b sigma~1) in the tiger and leopard had the lowest AIC value in comparison to the behavioral response models (b and bk) (Table 3). The density estimations are slightly higher in leopards from the EX detection function under the null model (Table 4). The formal AIC comparison indicated the HN detection curve model as the best model (Table 2). The density of the leopard was estimated slightly lower under the behavioral response model (b) compared to the null model and (bk) model (Table 1).

The density estimates based on the models with covariates were slightly higher 4.71 ± 1.20 (3.05–5.11) per 100 km² for tiger and lower for leopards 3.03 ± 0.78 (1.85–4.99) compared to null models (Tables 5 and 6). The models, including degraded forest type and moist deciduous forest type, were the most parsimonious models based on AIC weight for tiger and leopard (Tables 5 and 6).

Correlates of tiger and leopard density

The density of tigers and leopards were calculated using maximum likelihood-based SECR approaches under two broad model categories. The density estimates from the first category models (null and behavioral response models) in the tiger were not noticeably larger than each other. The AIC differences between the models using different detection functions (null model) in the tiger were also very small. No obvious inferences could be drawn based on null models. The leopard density was estimated slightly lower under the behavioral response model (b) than the null and (bk) model. The null model category performed all equal, and no meaningful conclusions could be achieved.

Thus, we created a second model category where the density was modeled with habitat covariates hypothesized to influence tiger and leopard densities in human-dominated habitats. We expected densities to be positively correlated with undisturbed habitats such as sal-dominated habitats and moist and dry deciduous forests. The undisturbed habitats are associated with high prey abundance (wild prey) and low anthropogenic influence. The densities of tigers and leopards have been reported to depend on prey biomass (Karanth et al., 2004b; Carbone & Gittleman, 2002; Khorozyan, Malkhasyan & Abramov, 2008; Carbone, Pettorelli & Stephens, 2011). Although tiger density was positively correlated with sal dominated habitats; however, contrary to our expectations, the tiger density was also positively associated with disturbed habitats. The positive correlation of tiger density in degraded habitats may likely be due to livestock availability in such habitats. In a recent study, the density of Asiatic lions (Panthera leo persica) was found to be negatively correlated with their natural prey in the tourism zone (disturbed) due to the assured food provisioning (Gogoi et al., 2020). A recent diet analysis indicates that more than 40% of the seasonal biomass consumed by tigers in the buffer zone of BTR was represented by livestock (Rather, 2020) compared to only 5% biomass consumption of livestock in the core zone of the reserve (Navaneethan et al., 2019). The density of leopards was positively influenced by road density, rugged terrain, and degraded habitats. The AIC comparison indicated the model with a moist deciduous forest to be the best density estimation model of leopard (Table 6). The high density of leopards relative to road density indicates the frequent use of forest roads as regular travel routes by leopards.

Contrary to the study of leopards in Thailand (Ngoprasert, Lynam & Gale, 2017) and Africa (Havmøller et al., 2019), our results indicate the frequent use of roads and disturbed habitats by leopards. Our results suggest that tigers and leopards in human-dominated areas frequently occur in degraded habitats relying partially on domestic prey species. Athreya et al. (2013) found a similar pattern of density estimates of leopard (4.8 ± 1.2) in human-dominated regions where leopards relied on domestic prey intake. The density estimates of leopard in this study are higher than the lowest estimates of leopard density (0.74 ± 0.18) recorded in India using the secr framework (Noor et al., 2020). Overall, the density estimates of tigers and leopards were lower compared to other studies in India. Karanth & Nichols (1998) estimated tiger densities to be (11.7 ± 1.93) in Kanha, (16.8 ± 2.96) in Kaziranga, (11.5 ± 1.70) in Nagarahole and (4.1 ± 1.31) in Pench Tiger Reserve. However, these estimates were calculated from the core zones that are relatively free from human interference.

Moreover, the estimates were calculated using traditional non-spatial CR methods that depend on extra buffer W around trapping grids. Recently, Sharma et al. (2010) compared the tiger densities using spatial likelihood methods in program density and using MMDM, 1/2MMDM as buffer strip, and based on home range radius calculated from telemetry, and found that densities were overestimated using 1/2MMDM. It is argued that in the absence of telemetry based home ranges, the spatial capture histories of camera traps should be used in likelihood-based density estimation methods (Borchers & Efford, 2008; Efford, Borchers & Byrom, 2008; Royle et al., 2008). The density estimation methods based on spatial likelihood approaches (SECR) do not depend on buffer widths around camera trap arrays. Thus the density estimates are least biased (Efford, 2004; Foster & Harmsen, 2012).

The observed density of tigers and leopards in the buffer zone may be regarded as moderate and not too low or high for a disturbed habitat. A literature survey indicated a huge disparity in the density estimates of leopards in India, varying from as low as 0.74 ± 0.18 (Noor et al., 2020) to the 28.9 ± 7.2 per 100 km² reported in Mudumalai (Kalle et al., 2011). The densities of tigers in India have been found to correlate with the prey densities (Karanth & Nichols, 1998; Karanth et al., 2004a). The observed density estimates in the buffer zones are likely. Leopards have been reported to occur relatively at a higher density in human-dominated landscapes (Athreya et al., 2007, 2013) in the absence of tigers. Athreya et al. (2013) reported that leopards occurred at a density of 4.8 ± 1.2 individuals per 100 km² in western India’s human-dominated landscapes. Kalle et al. (2011) reported a comparatively high density of leopards using SECR maximum likelihood (13.17 ± 3.15) and Bayesian models (13.01 ± 2.31) per 100 km² in Mudumalai Tiger Reserve. Leopards consume a wide variety of prey species, and even in prey deficient habitats are reported to occur at moderate densities by consuming small to medium-sized prey species (Hayward et al., 2006; Athreya et al., 2013).

Sampling adequacy

A maximum of 25 camera trap pairs were placed for eight consecutive days at any sampling occasion in 2 × 2 km grids before moving them to new locations. The presence of agricultural fields and human habitation within the buffer zone limited our inter-trap distance to 850 m. It is reported that in case of low tiger density (<2 per 100 km²), estimates of N can be achieved by increasing trap density. However, the estimates are likely to lack precision (Sharma et al., 2010). Silver et al. (2004) recommended the systematic approach of trap placement in medium to high tiger density areas with at least four trap nights km² and greater than 50 camera traps per 100 km². They also recommended a minimum sampling effort of over eight trap nights km². While studying ocelots, Maffei & Noss (2008) suggested that camera trap grids should cover at least a minimum of three to four average home ranges. The placement of camera traps in 4 km² grids spread over more than 500 km² of the trapping area and thus would have provided a fair chance for all tigers and leopards to be photo-captured. The 95 detections of 33 individual tigers and 74 detections of 38 individual leopards over 2,211 trap nights within the study area also indicate the area’s adequate sampling.