Leopard Panthera pardus fusca Density in the Seasonally Dry , Subtropical Forest in the Bhabhar of Terai Arc , Nepal

1 Department of Fish and Wildlife Conservation, Virginia Tech, Blacksburg, Virginia, VA 24061, USA 2 Eastern Himalayas Program (WWF-US), C/O WWF-Canada, 245 Eglinton Avenue E, Toronto, ON, Canada M4P 3J1 3 Nepal Engineering College-Center for Post Graduate Studies, Kathmandu 44600, Nepal 4WWF Nepal, Kathmandu 44600, Nepal 5 National Trust for Nature Conservation, Lalitpur 44700, Nepal 6Department of National Park and Wildlife Conservation, Kathmandu 44600, Nepal


Introduction
The leopard (Panthera pardus fusca Meyer, 1794) is one of the most widely distributed felids across the forested landscapes of the Indian subcontinent [1,2].Being a habitat generalist [3], the leopard has a wider fundamental niche than its larger congener, the tiger (Panthera tigris tigris Linnaeus), in terms of the habitat and area it occupies [4], extending from alluvial floodplains, subtropical deciduous moist and dry habitat in lowlands and Siwaliks, temperate deciduous forest habitat in mid hills and high mountains, to dry alpine forest in the Himalayas [5].Leopards occur sympatrically with tigers in Nepal, India, and Bhutan [6][7][8].
The tiger, being an apex predator, appeals to the public and serves as a flagship species [9].Perceived as more tolerant of anthropogenic influences, the leopard on the other hand has received less attention from conservationists despite its important functional role within ecosystems, including its potential to cause trophic cascades [10], its impact on mesopredators [11,12], and its competitive role within its guild [13,14].In the human dominated landscape of today's Indian subcontinent, habitat destruction and fragmentation remain major threats to leopards [4] and leopard numbers are declining due to both direct mortality and decreases in prey populations [15].In contrast to tigers, whose population sizes and trends have been intensively studied in multiple sites across the Indian subcontinent [8,[16][17][18][19][20][21][22][23], fewer studies have estimated leopard population sizes across both protected [24][25][26] and nonprotected areas [27], representing an array of habitat types in India.In Nepal, even fewer studies are available on leopard demography only representing alluvial floodplains, grasslands, and deciduous forest from Chitwan [28,29] and Bardia National Parks [30].However, there are no estimates of leopard density from seasonally dry subtropical forest in the Bhabhar region.Bhabhar is the alluvial apron of sediments washed down from the Siwaliks [31] and represents a key physiographic region extending to the "Terai zone" across the Terai Arc Landscape (Terai Arc).Hence, lack of information from Bhabhar habitat has hindered an overall assessment of leopard conservation status.
We estimate the abundance and density of leopards from the protected area within Bhabhar using a photographic capture-recapture sampling framework [32][33][34].This approach is widely acknowledged as a robust tool for estimating population abundance of elusive wild cats with individually distinct pelage patterns [15,20].We use traditional techniques of estimating the density using the ad hoc method of adding a buffer area around the polygon formed by connecting outer camera trap locations to account for edge effects (animals that do not occur entirely within trapping grid but have home range overlapping the edge of the grid).These buffering methods include using an estimate of home range radius derived from GPS or radio telemetry [35,36] or the common technique of using half of the mean maximum distance moved (1/2 MMDM) by animals within the trap array as a surrogate for home range radius [18,37,38].However, these ad hoc approaches have come under scrutiny because they are heavily influenced by small sample size, camera spacing, and extent of sampling grid relative to the animal's home range [35,39,40].To address the shortcomings of the traditional approach in estimating leopard density, we also use recently developed spatially explicit capture recapture (SECR) techniques that use spatial information more directly in the density estimation process [40][41][42][43].We compare the maximum likelihood [44] and Bayesian [45] SECR approaches in estimating density without the need to estimate an effective trapping area using the traditional ad hoc buffering approach.We present the result from both traditional and SECR approaches allowing us to compare our leopard density estimates with other studies in South Asia (Nepal, India, and Bhutan).

Materials and Methods
of Nepal.Encompassing over 499 km 2 , PWR is the largest wildlife reserve in the country and is contiguous to Chitwan National Park to the west.PWR is made up mostly of Churia hills (the outermost foothills of Himalayas) and Bhabhar regions, a rugged and highly porous landscape largely comprised of coarse alluvial deposits where streams disappear into permeable sediments [31].PWR has a monsoonal humid climate with more than 85% of the annual precipitation (2180 mm) occurring between July and October.The dry season occurs for 8 months between November and June [46].
The vegetation can best be described as subtropical, dry, deciduous forest with colonizing Saccharum spontaneum and Imperata cylindrica on the dry riverbeds and the floodplains, to a climax Sal (Shorea robusta) forest on Bhabhar and hillsides [47].The reserve supports a diverse mammalian fauna in addition to leopards, including carnivores such as the tiger (Panthera tigris tigris), dhole (Cuon alpinus), striped hyena (Hyaena hyaena), golden jackal (Canis aureus), Indian fox (Vulpes bengalensis), and ratel (Mellivora capensis).The principle wild prey species of the leopard include large size animals (>50 kg): gaur (Bos gaurus), sambar (Rusa unicolor), and nilgai (Boselaphus tragocamelus); medium size animals (20-50 kg): chital (Axis axis), muntjac (Muntiacus muntjak), and wild pig (Sus scrofa); and small size animals (<20 kg): common langur (Semnopithecus entellus) and rhesus monkey (Macaca mulatta).The combined ungulate prey density is estimated to be 6.6 individuals (±SE 1.1) per km 2 [48].Two of the settlements comprising approximately 100 households in Rambhori and Bhata in the core area of the reserve have been recently relocated [49].This is an event that is expected to trigger the recovery of carnivores within the reserve, making our density estimation an important baseline study.Illegal livestock grazing along the buffer zone is believed to reduce forage for wild ungulates, and livestock has been found grazing inside the reserve as far as 5 km from the reserve boundary.Photographic evidence from camera trap pictures [50] suggests that illegal poaching of wild prey is a direct threat to the carnivore populations in the reserve.
2.1.Field Methods.We conducted a camera trap survey for 64 days across a 289 km 2 core area of the reserve between December 2008 and March 2009.We followed the standard study design approaches prescribed for large felids at sites that had intensive signs of their usage [18,19,51].We first carried out extensive sign surveys for leopards [2,52] across 42 transect routes spread across the core area of the reserve amounting to 702 km searched on foot.These transect surveys enabled us to choose key locations for installing camera traps and identifying survey blocks that covered a large area without leaving potential gaps in our survey grid.
We selected 88 camera-trap locations spread throughout the study area based on the presence of leopard tracks, scats, scrapes, and other signs of use.To maximize capture probability, we positioned our camera traps along forest roads, trails, and dry stream beds, the habitat features known as leopard travel routes [6,53].The spacing between cameratrap locations (Figure 1) was maintained at approximately 1.9 km (±SE 0.06) [2].At each location, we used 2 passive digital camera traps (Moultrie D50, Moultrie Feeders, USA) activated by animal movement and placed on either side of the trail to photograph both right and left flanks of leopards [2,32].Each camera trap was active for 24 h and was checked on alternate days for proper recording of capture events, date, time, and any possible malfunctions.
We divided the study area into four "trapping blocks" each measuring an average area of 70.5 km 2 (±SE 7.54).We placed 22 camera stations per block resulting in an average of 32.2 camera stations per 100 km 2 .Each location was sampled for 16 consecutive days resulting in 16 encounter occasions, each consisting of capture data drawn from one day's trapping in each of the blocks [2].After 16 days, cameras were moved to the next block until all 4 survey blocks were completed.

Individual Identification.
Two investigators independently identified the photos for individual identification to build consensus on individual leopard identity.Individual leopards were identified based on their unique rosette patterns on the flanks, limbs, and forequarters [54] and given unique identification numbers (Figure 2) as done in other studies [24,26,27,55].

Population Estimation.
We followed the traditional capture-recapture analytic techniques used for estimating population sizes of large felids from remote camera data [2,18,19].We constructed capture histories for each individual leopard from photographic captures and assigned them to the appropriate encounter occasions (refer to Otis et al. [33] and Karanth and Nichols [2] for detail).We used two separate approaches to statistically test the assumption of population closure.
We first used the closure test implemented in program CAPTURE [56].We then used the Stanley and Burnham [57] closure test that assumes only time variation in recapture probability using the Pradel model [58] in Program MARK v. 5.1 [59].The Pradel model evaluates geographic closure by estimating the site fidelity (C), immigration (), and recapture probability () with regard to entry and exit into or out of the sampling area under assumption of the closure for the leopard population over our 64-day sampling period.We used the closed population models [33] implemented in program CAPTURE for estimation of overall capture probability ( p) and abundance ( N), using several different models that can incorporate effects of ecological and sampling-related factors (for details refer to [2,7]).We also used the Huggins closed capture with heterogeneity modeling platform in program MARK [33,59,60] to calculate abundance estimates.These models use a maximum likelihood framework and we fit 8 models of Otis et al. [33], which allow capture probabilities to vary over time (() = ()), by individual's heterogeneity ((ℎ) = (ℎ)), due to a behavioral response (initial capture being different from recapture probabilities, (⋅) ⋅ (⋅)) along with null model (with no variation in capture probabilities, (⋅) = (⋅)), and combinations of the above factors.Model input includes the one time winter season capture histories with 16 encounter occasions.We ranked all the models using sample size-adjusted Akaike's Information Criterion (AIC) [61] and considered all models with ΔAIC c < 2 as competing models [62].

Density Estimation.
To convert our abundance estimates from programs CAPTURE and MARK to densities, we used the traditional 1/2 MMDM [18,37] and full MMDM [63] approaches to calculate the buffer strip surrounding our camera traps to determine the effective trap area (ETA).The 1/2 MMDM and full MMDM were calculated from photographed individuals trapped in more than one location and we buffered each camera trap and dissolved the overlapping areas to calculate the ETA.We then divided our population estimates from CR (capture-recapture) analysis by total ETA to determine density.We used the delta method to calculate the variance in density estimates [64].
In program DENSITY v. 5.0 [41], we first modeled to select the appropriate detection (observational) process as either half-normal, hazard rate, or negative exponential.Using the selected detection function, we then allowed  0 (the capture probability at the hypothetical center of an individual's home range) and sigma (a function of the scale of animal movement) to vary using 2-class, finite mixture (ℎ2) to represent heterogeneity and/or a behavioral response ().Thus, a hazard detection function model with constant  0 and 2-class finite mixture of sigma would be represented as HZ  0 (⋅)sig(ℎ2).We used the estimated log likelihood and root-pooled spatial variance (RPSV) of varying integration buffers [41,69] for determining the appropriate buffer size.We ranked all the models using sample size-adjusted Akaike's Information Criterion (AIC c ) and considered all models with ΔAIC c < 2 as competing models.We used the model averaging techniques to determine final density estimates [62].We report the unconditional variance estimates for the model average estimates.
For the SECR-B approach, we used program SPACECAP [66] implemented in R package v. 3.0.1 [70] for estimating leopard density [65].We buffered 15 km around the sampling area to represent the probable extent of leopard home range centers and generated a grid of hypothetical home range centers with equally spaced points ( = 8150), each 0.336 km apart.This resulted in an area of 1,389 km 2 of leopard habitat over which these activity centers were uniformly distributed, after removing the 674 km 2 area of settlements (villages and agriculture areas: Rambhori, Bhata, Nirmal basti and amlekhganj).We used three standard input data files (animal capture locations and dates, trap deployment dates and locations, and hypothetical activity centers) and we assumed the half normal detection function.We performed 52,000 iterations, of which the initial 2,000 were discarded as the burn-in period, a thinning rate was set at 20, and we used an augmentation value of 180 individuals (more than five times the expected number of animals).We evaluated results using the Geweke diagnostic [71] and -score statistics of |-score| more than 1.6 implying lack of convergence [66].We produced the pixelated density map showing the estimated leopard densities per pixel of size 0.336 km 2 using ArcGIS 10.1.

Comparison of Leopard Estimates.
We compiled information on leopard densities from protected areas across their range in Nepal, India, and Bhutan and present cross-site comparisons of density and habitat type and describe density estimates based on type of modeling approach used and their corresponding standard errors.

Sampling Effort and Number of Individual Leopards
Captured.After discarding 66 trap nights of camera malfunctions, we amassed 1,342 trap nights and obtained 120 identifiable photographs comprised of 64 right flanks and 56 left flanks of leopard photographs (Table 1).Two investigators independently examined the photos for individual identification, both agreed on the 92% of all the capture events ( = 42), and we built the capture histories based on consensus ( = 39).We excluded the events for which consensus on the individual identity could not be built due to low quality pictures.We identified 19 individuals visually assessed to be 1 year or older (5 males, 12 females, 2 of unknown sex).

Closure Assumption.
Program CAPTURE closure test results were consistent with the assumption that the leopard population was closed during the 64-day survey period (Table 2).Additionally, the Stanley and Burnham [57] test was also consistent with the assumption of geographic closure during the sampling period as the model constraining site fidelity (C) to 1.0 and immigration () to 0, and underlying the  value to be 1, showed substantially more support than other models (Table 3).

Abundance Estimates.
The discriminant function analysis in program CAPTURE indicated that the  ℎ model, incorporating individual heterogeneity, was the best fit to the data (Table 2).The  ℎ (jackknife estimator) is widely acknowledged as a robust model in estimating abundance of large felids from the photographic capture-recapture analyses [2,18].The average capture probability () was 0.08 with an abundance estimate N(SE( N)) of 28 (±SE 6.071) leopards (Table 2).In program MARK, three models (  Spatially explicit models produced much lower estimates than traditional 1/2 MMDM techniques, but they were comparable to the full MMDM methods.The estimated log likelihood and RPSV of 245.04 units and 3,441 m suggested the appropriate buffer size of 15,000 m for SECR-ML in program DENSITY.Model selection in SECR-ML supported the hazard rate detection function with detectability (at home range center,  0 ) and spatial scale (function of movement, ) both constant (null model) (Table 5).There was some support for heterogeneity in the spatial scale as this model was within 2ΔAIC c of the top model.The model averaged density estimate for PWR based on SECR-ML was 3.78 (±SE 0.85) leopards per 100 km 2 .
Using the SECR-B analysis (SPACECAP program), all model parameters converged based on Geweke diagnostic statistics with  scores not more than 1.6 (Table 6).We obtained the posterior density estimates of the 3.48 (±SE 0.83) leopards per 100 km 2 and pixelated density map showing the relative animal densities over the animal home range centers (Figure 3).

Comparison with Leopard Densities in Other
Areas.The density point estimates in PWR via traditional methods (5.61 and 6.08) were nearly double the SECR estimates (3.48 and 3.78).However, even the lowest SECR estimate at 3.48 leopards per km 2 is highest reported so far in Nepal and is comparable to 3.45 leopards per 100 km 2 in the adjacent Chitwan National Park in Nepal.At the regional scale however, our estimates were lower in PWR than the estimates from India (Table 7).Using the traditional 1/2 MMDM approach, which we suspect overestimates density, leopard density in Bhutan was found to be low (1.01 leopards per 100 km 2 ) in comparison to 5.41 leopards per 100 km 2 .

Discussion
Our approach of conducting a sign survey [2] prior to the camera trap survey aided in fine tuning the survey protocol to maximize detection probability in PWR.In addition to relatively higher probability of capturing leopards present in the sampled area (67%:  +1 / N) and capturing 95% of the leopards in the first 12 days, our approach also demonstrated that operating cameras for short periods (16 days) per block and moving them sequentially to new blocks can provide reliable estimates of leopard abundance and density without violating the closure assumption during the total sampling period of 64 days [5,18,72].
PWR is the only protected area within the forest landscape of the Terai Arc that offers relatively pristine habitats for predators in the Bhabhar region.This study constitutes the first attempt to quantify leopard abundance and density from protected areas in the seasonally dry, subtropical deciduous habitat type and Bhabhar physiographic region in the Terai Arc.
Because of the prevailing controversy surrounding estimating animal densities, we used four different approaches to estimate leopard density in the study area.Of the four estimators, leopard densities from SECR-ML and SECR-B models were similar [66] but much lower in comparison to buffer strip method using 1/2 MMDM (Figure 4).Previous studies also indicated that traditional 1/2 MMDM methods overestimate density compared to SECR methods [22,25,26,40,67].While we found estimates using the full MMDM to be similar to the SECR models, there is no theoretical basis to justify this inference [40,73] since the 1/2 MMDM is meant to represent a home range radius for buffering camera traps.However, it is possible that the extent of the camera trap grid is simply too small to encompass the true leopard home range radius for multiple leopards and hence the 1/2 MMDM underestimates distance moved.For example, a study that used GPS collars and camera traps simultaneously on jaguars showed that actual home range radius was much larger than the 1/2 MMDM determined from camera traps [35].Therefore, the use of the SECR models may be better justified than traditional buffer strip methods in the absence of telemetry/GPS data.Understanding the potential pitfalls of traditional techniques, which may overestimate density, is crucial for managers as they require reliable estimates to make objective assessments of the conservation status of a species at risk.
Program DENSITY allows for selection from multiple possible observation models using an AIC framework and is computationally efficient and relatively quick [41] unlike SPACECAP [74].We did not evaluate the effect of "sex" as a covariate in density estimates due to inability to determine sex for some of the photographed individual leopards and the inability of program DENSITY to accept missing values.In addition, no feature in the current version of SPACECAP (GUI) allows us to incorporate missing data even though  it is computationally possible with Bayesian analysis [75,76].We reported both SECR density estimates in order to compare density estimates across studies that used these techniques [73,74,77].Also, we wanted to evaluate the tradeoff of being able to select among spatial detection models using an information theoretic approach in program DENSITY compared to choosing a spatial detection function in SPACECAP [74].We expected the posterior estimates from program SPACECAP to be superior as it also reports interval estimates of density as direct probabilities without problematic asymptotic assumptions [27,45].The benefits of SPACECAP include that its nonasymptotic assumptions appear to fit the low sample size camera trap data better [45] and it can be extended to open models in future studies.Ultimately we found the estimates from the two SECR methods to be very similar with similar precision.
Our results show that PWR harbors high leopard density similar to other protected areas, namely, Bardia and Chitwan National Park, areas that are considered better habitat.High prey biomass may be the primary reason for unexpectedly high leopard densities in PWR [15], where ungulate prey occurred at 6.6 (±SE 1.1) individuals per km 2 [48].The current densities of medium-and large-sized prey in PWR are probably not sufficient to support a high density tiger population, assuming an annual kill rate of 3000 kg per year per tiger and annual biomass cropping rate of 10% [6].As a result, tiger density in the reserve is low, at 0.87 tigers per 100 km 2 [48], and leopards are, perhaps, able to fully exploit their preferred prey, medium-and small-sized prey, perhaps explaining the high density of leopards in PWR.At the regional scale, results of this study demonstrate that SECR leopard densities of 3.48 per 100 km 2 in seasonally dry deciduous forest were lower than SECR estimates in similar sites in India.Prey availability and the broad niche width of leopards in forested landscapes could be the reason for their variable densities in the subcontinent.
Baseline population estimates for leopards are essential for monitoring the effectiveness of conservation initiatives [78] and for guiding conservation decisions.In this regard, our findings provide useful insights for leopard conservation both at the local and regional levels.At the local level, our results serve as the benchmark data to assess conservation impacts of the relocation of ca 100 households from the core area of PWR.At the regional scale, leopard demographic assessment across this large conservation complex of the Chitwan-Parsa-Valmiki Tiger conservation unit [79,80], which was established for tigers, can simultaneously provide reliable information on leopard density for the transboundary conservation landscape, following the model provided by the transboundary Manas conservation landscape across India and Bhutan [81].Hence, long-term monitoring of leopards with camera trap studies, as employed in our study, would be useful to understand conservation status and variation in the population sizes over time and from local to regional scales, thereby informing decision makers in implementing sound conservation management recommendations.
The current Nepalese legislation on wildlife does not recognize the leopard on the protected animal list [82], yet it is one of the most heavily traded species for its skin [83].The antipoaching strategy designed to protect the tiger and rhinoceros as flagship species in the region [84,85] should also aid in leopard conservation within protected areas.However, the large swathe of unprotected areas in the human dominated landscape of the Terai Arc [86][87][88] represents a considerable challenge in terms of the impact of habitat fragmentation [4] and increasing human wildlife conflict in protecting this iconic mesopredator.

Figure 1 :
Figure 1: Study areas showing the spatial location of camera traps in the Bhabhar region and the effective sampling area formed by drawing a minimum convex polygon surrounding the outermost camera trap locations.The area of MCP (minimum convex polygon) is 289 km 2 .

Figure 2 :
Figure 2: Identical pelage patterns from the same leopard in PWR as displayed in photographic captures from different camera traps.

Figure 3 :
Figure 3: A pixelated density map showing relative leopard densities per pixel of size 0.336 km 2 .The area of MCP (minimum convex polygon) is 289 km 2 .

Figure 4 :
Figure 4: Comparative density estimates of leopards in PWR using traditional density estimators that buffer trapping grid by 1/2 or full mean maximum distance moved (MMDM) among camera traps and recently developed spatially explicit capture recapture (SECR) models in a maximum likelihood (SECR-ML) or Bayesian (SECR-B) framework.

Table 1 :
Summary statistics for photographic capture-recapture data on leopards in PWR.

Table 2 :
[33]ure test and model selection (after Otis et al.[33]) for leopard population size estimation using photographic capture-recapture data from Parsa Wildlife Reserve in Program CAPTURE.DFS is the discriminant function score.The   shows the Behavioral effect, while  ℎ shows the heterogeneity effect and the   shows the constant null model.

Table 3 :
Model selection summary for geographic closure for the leopard population in Parsa Wildlife Reserve in program MARK.Phi represents the site fidelity. is the recapture probability and  is immigration onto the study site.Parameter with "⋅" indicates a constant value and with "1" and/or "0" indicates parameter is fixed.AIC c is Akaike's information criterion corrected for small sample size and difference in AIC value between top model and th model is represented by ΔAIC c .Weight of support for each model is AIC c weights, .

Table 4 :
[59]ard abundance from Parsa Wildlife Reserve was estimated using Huggins closed captures with heterogeneity models in program MARK(Pledger [60], White and Burnham[59]).[(⋅) ⋅ (⋅)] is the behavioral effect.ℎ [(ℎ) = (ℎ)] is the heterogeneity effect, while   [(⋅) = (⋅)]is the constant null model.AIC c represents Akaike's information criterion corrected for the small sample size.Difference in AIC value between top model and th model (ΔAIC c ). Weight of support for each model (AIC c weights, ).N represents abundance estimate, while SE represents standard errors.Model averaged density estimates and standard errors are given in bold.* represents model selected for inferring parameter estimates.

Table 4
3.4.Density Estimates.Six individual leopards were captured more than once resulting in a MMDM of 5.2 km and an effective trapping area of 473 km 2 using the buffer strip

Table 5 :
Model selection results from leopard density estimates using photographic capture-recapture data from Parsa Wildlife Reserve in program DENSITY using hazard rate detection function. 0 is the capture probability at home range center. is the spatial scale parameter of capture function.ℎ2 is the 2-class finite mixture probability for heterogeneity.Akaike's information criterion adjusted for small sample size by AIC c .  represents Akaike weight, while   is estimated density (per 100 km 2 ) and SE represents its standard error.Model averaged density estimates and standard errors are given in bold.

Table 6 :
[66]posterior summaries from Bayesian spatially explicit capture-recapture (SECR-B) of the model parameters including the leopard density estimates from Parsa Wildlife Reserve implemented in SPACECAP[66]along with Geweke diagnostic parameters.Sigma is the range parameter of the species.lam0 is the intercept of expected encounter frequency.psi is the ratio of the number of animals present within the state space, , to the maximum allowable number.Nsuper is the number of activity centers in .Density (per 100 km 2 ) is Nsuper divided by  and | score| greater than 1.6 implies lack of convergence.

Table 7 :
Leopard (Panthera pardus fusca) density estimates (per 100 km 2 ) across the array of habitat types from study areas in South Asia based on traditional mean maximum distance moved (MMDM) approaches and spatially explicit capture-recapture (SECR) using maximum likelihood (SECR-ML) and Bayesian (SECR-B) analytical methods.1/2 MMDM and full MMDM approaches were used in population estimates from program CAPTURE.
# based on combination of telemetry studies, camera trapping, and leopard tracks.Estimate is representative of southern most areas of Bardia National Park.$ A nonprotected area.§ Population estimates based on program MARK for density estimates.