Introduction

The monitoring of large carnivores living in mountainous ecosystems is a formidable challenge for conservationists and managers. For example, carnivores inhabiting mountainous landscapes tend to persist at lower densities than carnivores in lowland habitats (Alexander et al. 2015; Ghoddousi et al. 2010; Kachel et al. 2017) and have large movements and home ranges (Cheraghi et al. 2019; Farhadinia et al. 2018b; Johansson et al. 2016). Spatially explicit capture recapture (SECR) models are widely used to monitor populations of large carnivores (Sharma et al. 2014; Sollmann et al. 2013). SECR models incorporate spatial locations of captures within a unified model to provide reliable estimates of density (Borchers and Efford 2008; Efford 2004). They estimate animal density from a set of individual animal detections made at capture locations, for example by means of motion-detector camera traps, nested within a broader network of potential home-range centers (Efford 2004).

However, when the population size is small, several challenges arise for SECR density estimates (Gerber et al. 2014). Most importantly, infrequent detections (Mohamed et al. 2019; Rostro-García et al. 2018), the low number of individuals detected (Hearn et al. 2017; Sharma et al. 2014) and variability in detection among individuals (Alexander et al. 2015; Gerber et al. 2014) can compromise SECR models’ precision. The SECR method also requires sampling to be conducted relatively densely to obtain sufficient recaptures of individuals on multiple cameras; this is required for adequate estimates of the spatial scale parameter (σ) (Sollmann et al. 2013; Wilton et al. 2014).

Nonetheless, SECR models perform well for species with large home ranges that are present at low densities if some realistic requirements are met (Wilton et al. 2014; Zimmermann and Foresti 2016). First, the extent of the detector array has to be similar or larger than the extent of individual movement (Efford 2011; Sollmann et al. 2012). Second, the optimal distance between traps is twice as large as the σ (Sollmann et al. 2012; Sun et al. 2014). Third, repeat detections should exceed 20 (Efford 2011), and finally, each individual should be detected on average at least 2.5 times (Gerber et al. 2014). Importantly, sharing detection information across time and space (Gerber et al. 2014; Morehouse and Boyce 2016) can increase the precision of density estimates. However, when landscape features restrict the spatial arrangement of detectors, notably in mountainous areas, biologists inevitably select small-sized study areas, usually in good quality habitats, which can result in violating some of the above requirements (Suryawanshi et al. 2019).

Asian mountains harbor large felids with large spatial requirements, such as common leopards (Panthera pardus) and snow leopard (P. uncia) (Farhadinia et al. 2018b; Johansson et al. 2016). Hence, most population estimates of large felids in mountainous areas, based on SECR photographic data, come from study areas that are small relative to the species home range (Alexander et al. 2015; Farhadinia et al. 2019; Kachel et al. 2017; Suryawanshi et al. 2019), usually with small populations of fewer than ten individuals (Ghoddousi et al. 2010; Kachel et al. 2017; McCarthy et al. 2008). Their rarity and potential vulnerability of these populations to demographic stochasticity underscores the critical and time-sensitive need to develop efficient analytical methods to track their populations robustly over time (Karanth et al. 2006).

In this study, we employed multi-session SECR modeling to understand the population density and composition of Persian leopards (P. p. saxicolor) in a mountainous desert environment in central Iran. We expected that leopards would show inter-seasonal variation in density parameters, with larger movements in winter than summer. Our study provides the first population density estimate of leopards from the driest area in which they have ever been studied globally.

Materials and methods

Study area

We implemented our study in the Bafq Protected Area (hereafter Bafq) of central Iran (31.5023° N, 55.7134° E). Bafq was designated a protected area in 1996, near the city of Bafq in central Iran. With an area of 885 km2, Bafq is an arid mountainous region typified by sparse plains and rolling hills (Fig. 1). The altitude range varies between 1060 and 2860 m above sea level and mean annual precipitation is 70 mm (Sohrabinia and Hosseini-Zaverei 2010). Springs, wells, air pumps and small dams constructed at high altitude are the main water supplies in the region. There are four main villages inside the area. Anti-poaching law enforcement is carried out by 13 rangers.

Fig. 1
figure 1

Spatial configuration of study areas and locations of camera trap stations across the two sessions in Bafq Protected Area, central Iran. The map inset shows locations of the study area in Iran

Full size image

Sampling design

We deployed camera trap stations (Panthera® IV and V (New York, NY 10,018, USA) and Cuddeback Capture Model 1125 (Non Typical, Inc., Park Falls, WI, USA), all working with white Xenon flashes) on park-wide 4 × 4 km grids to ensure even coverage of the whole study area (Fig. 1) for 88 and 91 days in winter 2011–2012 and summer–autumn 2016 within Bafq, respectively (Table 1). Short sampling periods (≤ 3 months) are necessary to avoid violating the assumption of demographic closure (Zimmermann and Foresti 2016). In winter 2011–2012, camera traps were placed either along trails or dirt roads as these are regularly traversed by leopards (Farhadinia et al. 2019). In 2016, when the survey was conducted during the driest months of the year (August–October), we equipped 38 water resources (springs or artificial waterholes, 58% of stations) with a camera trap, additional to the trails (n = 28, 42% of stations), predominantly along valley bottoms or trails. Camera stations were placed at a mean spacing of 1620 and 1409 m in the two consecutive sessions to simultaneously achieve the twin objectives of maximizing the number of individuals caught and adequately recapturing individuals at different camera traps, as required in SECR designs. A total of 22 and 26 grids were sampled during the 2 sessions, respectively (Fig. 1).

Table 1 Details of sampling design for the multi-session spatially explicit capture–recapture framework in Bafq Protected Area (2011–2012 and 2016)
Full size table

Data analysis

We considered each 24-h day as a sampling occasion. All consecutive photographs of the same individual taken no more than 0.5 h apart were defined as a single event (Braczkowski et al. 2016). The identity of leopards was determined by the unique rosette patterns on their pelage, independently by two researchers (PB and MSF), and individuals were only included in the analysis if consensus on individuals was achieved by both observers. The sex of individual leopards was determined where possible from sex-specific cues, such as visible genitalia, the presence of dependent individuals. We only included independent individuals, including adults and sub-adults, in the analysis. Dependent cubs were excluded from the SECR modeling.

We estimated the density of leopards using a maximum-likelihood SECR approach (Borchers and Efford 2008) using the package ‘secr’ version 3.2 (Efford 2019) in the R software version 3.3.3. (R Development Core Team 2013). The estimates of population density produced by SECR are less sensitive to the edge effects of traditional non-spatial estimators, incomplete detection and heterogeneous capture probabilities, and eliminate the need for an ad hoc estimation of the sampling area (Efford 2004). The secr package also allowed us to evaluate the effect of sex on density parameters.

We combined the 2 years into a single model framework and considered each year as an independent session. For each ‘session’, secr modeling assumes a new realization of the underlying population process (Efford 2019). The multi-session analysis enables the fitting of models with parameter values that apply across sessions, and the data are then effectively pooled with respect to those parameters. Our two sessions did not overlap temporally and only one individual (F2) was observed in both sessions, ensuring that the independence assumption between the sessions was not violated.

We used a hierarchical model composed of an explicit state-space process model and an observation model (Efford 2004). The animal population size and their respective central locations (“home-range centers”) constitute the state-space process, assuming a Poisson distribution (Borchers and Efford 2008; Efford 2004). The observation model describes the probability of encounter as a function of an individual’s location at the time of sample, and a probability of “count” detector parameter (Efford, 2017). The half-normal detection function contains two parameters: g0, which is the baseline detection rate when the distance between the animal’s activity centers and the camera traps is zero; and σ, which is the spatial scale parameter (with the unit of meters) of the encounter probability model (Borchers and Efford 2008; Efford 2004).

We defined the area of integration (i.e., state-space model) by equally spaced points in a regular grid, with a mesh spacing of 1 km2. A buffer was plotted around the detector array to incorporate individuals with activity centers outside of the trapping area, but whose movement range extends into the sampling area (Borchers and Efford 2008; Efford 2004). We used a mask with a buffer of 120 km to define the outer limit of the state-space area. This buffer size (3–4 times σ) corresponds to an area beyond which the animals would have a low probability of being detected, and thus was unlikely to influence density estimates (Efford 2004). As leopards are unlikely to have their home-range centers in these areas of low productivity or increased conflict (Braczkowski et al. 2016; Farhadinia et al. 2019), areas such as villages, open deserts and sand dunes were masked out from the effective area, based on an ensemble suitability map developed by Ahmadi et al. (2020).

For our models of density estimation, we fitted 14 a priori models with varying effects on g0 and σ in a multi-session framework (Table 2). We included a trap-specific behavioral response in the baseline detection model because we expected that the leopard behavior to change after being detected at a specific trap for the duration of the session (bk). Large terrestrial carnivores typically feature differences in their home-range sizes, movement patterns and capture probability (Farhadinia et al. 2019; Sharma et al. 2014). This can affect the observation process in capture–recapture exercises (Sollmann et al. 2013). SECR is considered a robust method for calculating animal densities because it accounts for individual heterogeneity due to space usage (whereas non-spatial capture–recapture does not). To account for heterogeneous capture probability in leopards, we also included a sex-specific covariate in the observation process, which would account for different capture probabilities amongst males and females (Farhadinia et al. 2019). We included a sex-effect for the encounter rate (g0) and spatial scale (σ), both of which are widely seen in leopard studies (Braczkowski et al. 2016; Goldberg et al. 2015; Rostro-García et al. 2018). We finally fitted session-stratified estimates, meaning that all parameters vary across sessions (i.e., years), by maximizing the likelihood. Akaike’s Information Criterion corrected for the small samples (AICc) was used to identify the most parsimonious model (i.e., the lowest AICc score; (Burnham and Anderson 2002). The highest ranked model was used to estimate leopard density (D), population size, detection rate at home-range center (g0) and σ. The population size was the number of individuals within the region for the current realization of the process (Efford 2017).

Table 2 Model selection results for 2 sets of 14 fitted multi-session models (totaling 28 models) ranked by Akaike’s Information Criterion corrected for the small samples (AICc) for 2012 and 2016 Persian leopard density Bafq Protected Area, central Iran
Full size table

Results

Photographic captures

During a sampling period of 2823 trap nights in winter 2011–2012, we obtained a total of 49 independent leopard detections (events) of 8 independent leopard individuals. Five leopards were photographed on both flanks, two male leopards were photographed only on their right flank and one female leopard accompanied by a single cub was photographed only on its left flank. For analyses, we used eight different leopards, equally from each sex (Fig. 2). In summer–autumn 2016, we obtained 141 independent detections of leopards during 3901 trap nights (Table 3). Since 58% of camera traps were deployed at water resources in the summer of 2016, leopards regularly turned their bodies which enabled us to photograph both flanks. From these detections, we were able to identify five independent individuals. They included only a single adult male and four adult females (Fig. 2). In each session, we also detected two leopard families with a total of three cubs (ranging 1–2 per family), resulting in a cumulative number of four families with six cubs for the study period.

Fig. 2
figure 2

Comparison of detection frequency for all demographic classes between winter 2011–2012 and summer–autumn 2016 sessions in Bafq Protected Area, central Iran. Each code on the x-axis refers to a single individual leopard. Only one individual (F2) was detected in both sessions

Full size image
Table 3 Details of baseline information on leopard based on systematic camera trapping across two sessions of winter 2011–2012 and summer–autumn 2016
Full size table

Leopard detections were approximately three times more frequent at water-based cameras compared with trail-based stations comparing to proportions of each camera placement (105 vs. 36 independent detections; X2 = 5.2, df = 1, P = 0.02). The four independent females observed during the second session were detected only 1–4 times at trail-based camera traps, while they were more frequently seen at water resources, varying between 5 and 24 detections. During the winter session, same-trap recapture was rare (10% of total detections) whereas it increased to 27% of detections during the summer session. Male leopards showed large-scale movements relative to the size of the detector array.

Density estimation

There was strong support for one model based on the AICc (Table 2), defined as g0 ~ bk + session + sex σ ~ session + sex. We found sex and session-specific variability in spatial scale (σ). Males had larger σ than females in each session, yielding an inter-sexual ratio of 1.6 (Table 4), indicating large spatial movements of male leopards throughout Bafq. Likewise, although σ estimates were overlapping in their 95% confidence intervals between the two seasons for each sex, the mean estimates of σ for the winter session were 1.5 larger than summer–autumn session, regardless of the sex (Table 2). This suggests less mobility of leopards during warm months.

Table 4 Density estimates of Persian leopards with standard error (SE) and 95% confidence interval (lower and upper) of parameters for spatially explicit capture–recapture models fit to camera trapping data from Bafq Protected Area, central Iran
Full size table

The top model also supported sex and session-specific baseline detection rates (g0). Thus, the g0 estimates showed substantial inter-sessional difference for each sex, increasing by a magnitude of 21.4–30.0 from the winter to summer–autumn sessions (Table 2). For all sex groups within a session, g0 bk=1 was larger than g0 bk=0 with a ratio of 3.2–4.3 between the means. The bk effect tested the hypothesis that leopard behavior changes after being detected at a specific site for the duration of the survey (trap response). Therefore, when a leopard was detected in a specific camera trap site, the probability of a subsequent encounter for the entire survey was increased, i.e., the individual became ‘trap happy’ at the population level (beta coefficient g0 bk = 1.4, 95% CI = 0.9–1.8).

The best performing model estimated a density of 1.6 (95% CI = 0.9–2.9) and 1.0 (95% CI = 0.6–1.6) independent leopards/100 km2 for the first and second sessions, respectively (Table 4). A total of 10.9 ± SE 3.4 independent leopards were estimated for the winter 2011–2012 while the population size was calculated as 6.9 ± SE 1.8 independent leopards for the 2016 session.

Discussion

We estimated the density of leopards from the driest known study site of the species global distribution, with an average annual precipitation of 70 mm (Sohrabinia and Hosseini-Zaverei 2010). Although data on leopard density estimates in arid regions are sparse, they, including the present data, indicate some of the lowest densities across the leopard global range, all below 2 adults per 100 km2 (Edwards et al. 2016; Ghoddousi et al. 2010; Mann et al. 2020; Stein et al. 2011).

Density patterns of leopards over time

Temporal changes in population density and composition in Bafq’s leopards support two types of variation. First, substantial overlap in confidence intervals between the two sessions suggested that the different density point estimates could be due to sampling variation, i.e., variation in a statistic from sample to sample, which is commonly seen in multi-session SECR studies of leopards (Rogan et al. 2019; Rosenblatt et al. 2016). Second, detection of only one individual in both sessions (4 years apart; Fig. 2) suggested that despite the observed overlap between the large confidence intervals, density variation, i.e., real population differences, cannot be totally rejected.

Accordingly, there are three possible reasons to explain observed variations between the two sessions. First, prey density as the main bottom-up driving process of carnivore populations can shape leopard demography (Ramesh et al. 2017; Strampelli et al. 2018). We obtained annual census data of ungulates in the study area for the period between 2011 and 2016, conducted by Yazd Provincial Office of the Department of the Environment every November using 10–15 groups of field staff (Table S1). These census data suggested that the population size for bezoar goat Capra aegagrus and urial Ovis vignei, two key prey for Persian leopards (Farhadinia et al. 2018a; Rezaei et al. 2016), was not decreasing.

Second, top–down regulation due to human persecution, particularly due to conflict with human communities can cause high anthropogenic mortality (Rosenblatt et al. 2016; Rostro-García et al. 2018; Sharma et al. 2014). A dietary investigation, based on fecal analysis, found only 7.54% of the total leopard diet (expressed as biomass consumed) to be composed of domestic animals in Bafq (Rezaei et al. 2016). Equally important, only one record of an individual poisoned between the two survey efforts was obtained (Seyed Jalal Mousavi, pers.obs). Unlike many parts of the leopard range in west Asia where conflict with human communities is a major source of concern (Babrgir et al. 2017; Khorozyan et al. 2020; Sharbafi et al. 2016), we found little evidence, based on poaching records and dietary analysis, that the Bafq leopard population suffers direct anthropogenic persecution. Nonetheless, almost no replication in leopard individuals between the two sessions suggests that the real anthropogenic mortality of leopards is likely to be higher than what was detected.

Finally, emigration can affect the demography of large felids (Fattebert et al. 2015; Sharma et al. 2014). Persian leopards can disperse up to 80 km (Farhadinia et al. 2018b), implying that many surrounding scattered mountains around Bafq can be potentially visited by leopards. We expect that emigrated individuals have much less chance of survival due to a lower level of law enforcement, if any, outside Bafq Protected Area.

In conclusion, we expect that additive effects of sampling variation (because of different seasons and sessions) and density variation (mostly due to human persecution and emigration) resulted in the observed density pattern over time. However, the small sample sizes and a lack of further replications requires us to view our findings as suggestive rather than conclusive, and to concede that further research is necessary to assess the demographic trends of leopards in central Iran.

Seasonality and conservation implications

Importantly, season had a strong effect on baseline detection rate (g0) and spatial scale (σ), both of which influence the density estimate. Nonetheless, season can have confounding effects with year. Larger σ during winter (2011–2012) suggested a higher mobility of leopards preceding and during the mating season, which peaks in mid-winter (Farhadinia et al. 2018b). In contrast, water resources are crucial for leopards in warm seasons, not only to meet their drinking requirements, but to optimize their encounter rate with prey (Farhadinia et al. 2018a). This results in smaller σ and larger g0 in the warmer months of the year. In addition, leopards exhibited an increased tendency to revisit camera trap stations, a result that is unsurprising when 58% of traps are at water resources and successive captures are not independent. In addition to season, leopard populations may exhibit marked inter-annual variation in σ (Rogan et al. 2019; Rosenblatt et al. 2016), possibly due to their dynamic spatial patterns and home-range sizes.

Population monitoring of large carnivores, particularly in remote montane landscapes, represents a considerable conservation challenge (Farhadinia et al. 2018c; Suryawanshi et al. 2019). The concentration of cameras around water resources in hot months can result in a higher recapture rate and consequently higher precision in density estimates. Future studies are encouraged to investigate the potential advantages and disadvantages of preferential water-based sampling on density estimation. Finally, season can have strong effects on parameters for SECR models in mountainous areas with remarkable seasonal pattern. A multi-session SECR framework can enhance accuracy of density parameters.

The persistence of small leopard populations can be adversely affected by demographic stochasticity through reducing survival and fecundity (Rostro-García et al. 2018; Williams et al. 2017). Importantly, controlling two threats, prey depletion and leopard persecution in nearby areas to support individuals dispersing from this breeding population are crucial for the future persistence of this population. In addition to single density estimates, temporal changes in population density and composition, similar to those we provide here are helpful to inform data‐driven decision‐making processes to define conservation action priorities. Most importantly, and to the best of our knowledge, no study prior to ours has assessed the demographic trends of leopards, and other large carnivores inhabiting the remote drylands of the Middle East. Our study has provided an important and first step in addressing this knowledge gap, and this is essential for the future management and conservation of leopards.