Elusive cats in our backyards: persistence of the North Chinese leopard (Panthera pardus japonensis) in a human-dominated landscape in central China.

The North Chinese leopard (Panthera pardus japonensis), the least-known big cat, disappeared in most historical range for decades, following the development of modern civilization. Unfortunately, we have scarce knowledge about the status of this big cat so far, apart from anecdotal reports. In this study, we investigated density, distribution, and habitat use of the leopard, the apex predator, in a complex forest landscape in the Loess Plateau. We used a camera-trapping network to obtain population estimates for leopards over 2 years through spatially explicit capture-recapture models. Our results, based on maximum likelihood and Bayesian/MCMC methods, reveal that the largest wild population of the leopard was found widely distributed in remnant forests in central Loess plateau. The population is increasing in our study area, and the density of leopards (1.70 (SE = 0.48) - 2.40 (SE = 0.67)/100 km2 ) is higher than other areas of China. According to the analysis of 2 seasonal occupancy models, prey species drive partially the leopard habitat use, predicting that the big cat thrives from the recovery of prey community. However, human disturbances, especially oil wells, seem to have negative impacts on the habitat use of leopards. Specifically, it is necessary to have joint efforts by the government and researchers to improve human disturbances management and prey species population density, as well as strengthen the investment in research on the North Chinese leopard, which could all further strengthen protection ability and ensure the long-term survival of this species.

INTRODUCTION

Large terrestrial mammalian carnivores serve as keystone species and aid in the conservation of biodiversity (Harihar et al. 2011). Populations of such carnivores decline due to habitat loss, poaching, prey depletion resulting from human activities (Ripple et al. 2014), and high metabolic demands make them highly prone to conflict with humans (Treves & Karanth 2010). Population declines of such species below certain thresholds can lead to oversimplified ecosystems (Terborgh et al. 2001; Soulé et al. 2003). Strong inferences of the density and distribution of many terrestrial mammalian carnivores are crucial for their conservation (Karanth et al. 2004, 2006). However, large mammalian carnivores often occur at low densities and exhibit wide home ranges, which make their densities and distributions challenging to study (O'Connell et al. 2010).

The leopard (Panthera pardus) has the largest distribution of all wild cats across Asia and Africa (Nowell & Jackson 1996; Sunquist & Sunquist 2002; Stein et al. 2016). Leopards can tolerate a wide range of habitats, but suffer nevertheless from human activities and habitat loss (Athreya & Karanth 2011; Hebblewhite et al. 2011; Swanepoel et al. 2015; Stein et al. 2016), which has led them being listed as endangered, critically endangered or vulnerable on the IUCN Red List (Jacobson et al. 2016). According to phylogenetic analysis, the leopard is partitioned into nine subspecies: P. p. pardus, P. p. nimr (critically endangered), P. p. saxicolor (endangered), P. p. fusca, P. p. kotiya (endangered), P. p. delacouri, P. p. japonensis, P. p. orientalis (critically endangered), and P. p. melas (critically endangered) (Uphyrkina et al. 2001). Four subspecies are distributed in China: (i) P. p. orientalis, Northeast China; (ii) P. p. delacouri, southern China; (iii) P. p. fusca, Tibet, China; and (iv) P. p. japonensis, northern China (Miththapala et al. 2010; Song et al. 2014; Alice et al. 2015). Historically, leopards were distributed throughout China, but recent research has reported that leopards may only occur in about 19 (out of 34) provinces (Bao et al. 2010; Alice et al. 2015). However, we still lack information of the density and distribution of leopards in China, except for the Amur leopard (P. p. orientalis) (Wang et al. 2016, 2017, Vitkalova et al. 2018).

The North Chinese leopard (P. p. japonensis) is only distributed in North China (Miththapala et al. 2010; Song et al. 2014). Population study on the North Chinese leopard is very sparse, and almost all available data are derived from unproven information provided by local people (Gao et al. 2007; Liu et al. 2007; Alice et al. 2015). Specifically, there are 11 researches on wild population, among which 2 researches used camera trap technology and indicated how many individuals they identified, while others were based on line transect or access investigation (Song et al. 2014; Xie 2019). This subspecies has been ignored for a long time. Nevertheless, as the only apex predator in central Loess plateau, leopard is crucial for this weak ecosystem. Leopards are adapted to various environments (Nowell & Jackson 1996), due to wide range of food habits (Hayward et al. 2006). Precious studies suggested that leopard would be relatively more abundant because of their ability to survive on medium‐sized and smaller prey (Rabinowitz 1989; Seidensticker et al. 1990), and vice versa, if leopard become scarce in a given area, it proves that the environment is so fragile that even small animals cannot survive (Rabinowitz 1989; Seidensticker et al. 1990). Therefore, the lack of information on the density, distribution, and habitat use of North Chinese leopards means that protection status of ecosystem of central Loess Plateau cannot be effectively evaluated (Miththapala et al. 2010), which in turn precludes the protection of this subspecies, as well as other species coexist. Thus, our study is essential to provide the first comfirmed assessment of the North Chinese leopard wild population on density, distribution, and habitat use, for further conservation strategy and policy.

Camera trap surveys are non‐invasive methods that are widely used in wildlife ecology and conservation studies (Cutler & Swann 1999), especially for felids (Karanth & Nichols 1998; Wang et al. 2017). As the cost of camera equipment have decreased, camera trap surveys are used to study not only wildlife abundance but also their distributions (Wang et al. 2016) and behavior (Yang et al. 2018b, 2019). Precise estimation of population size is crucial for wildlife management and conservation (Stephens et al. 2015), and camera trap surveys suffer, as any wildlife survey method, from some limitations such as imperfect detection (Stephens et al. 2015). However, many spatially explicit models have recently been developed to account for such sampling errors (Mackenzie et al. 2002; Efford et al. 2009; Royle et al. 2009), allowing direct estimates of the density and resource use of cryptic or low‐density species by using small data sets from camera trapping surveys (Carter et al. 2015).

The objectives of this study were as follows: (i) to conduct the first camera trap survey aimed at estimating the North Chinese leopard density in Shaanxi Province in central China; (ii) to estimate population size and distribution; and (iii) to examine how prey species and human disturbance affect leopard habitat use.

MATERIALS AND METHODS

Study area

We conducted our study in Shaanxi Province, central China (35°30′−36°01′N, 108°30′−108°49′E) (Fig. 1). The area is in the hinterland of the Loess Plateau, and the Ziwuling Mountain to the west constitutes the border between Shaanxi and Gansu provinces. The average annual temperature range is 3.42–11.9 °C, the average annual sunshine duration is 2159.4 h, and the average annual precipitation range is 588.7–677.4 mm, with the rainfall mainly concentrated between June and September (Liu 2004; Zhang 2014). The study area is situated in a rugged, mountainous landscape with an altitude ranging from 1100 to 1750 m a.s.l. (Zhang 2014). The vegetation type is mainly temperate deciduous broad‐leaved forest (Zhang 2014) and mixed coniferous forest at high elevations. Two nature reserves are in the study area: the Ziwuling Nature Reserve (ZNR) in the north and the Qiaoshan Nature Reserve (QNR) in the south. The leopard prey species include the Siberian roe deer (Capreolus pygargus), the wild boar (Sus scrofa), as well as some small mammals such as the Asian badger (Meles leucurus) and the Tolai hare (Lepus tolai). Other predators, including the leopard cat (Prionailurus bengalensis), the red fox (Vulpes vulpes), and the yellow‐throated marten (Martes flavigula), coexist with North Chinese leopards in our study area (Zhao et al. 2020).

Figure 1

Map of the area sampled by camera traps in Shaanxi, China, during 2016–2017, showing the camera sites relative to oil wells, major roads, and nature reserves. The inset shows the location of the study site within China.

Because the ZNR is a national nature reserve, human activity is very rare, and forestry workers patrol frequently. In contrast, the QNR is a provincial park, where human activity and oil operations are common. Unlike in other areas in China (Alexander et al. 2015; Wang et al. 2017), livestock grazing is not very common in our study area.

Data collection and field methods

Based on Karanth and Nichols (1998), we divided the study area into 4 × 4 km cells and selected 1–2 camera trapping sites in each cell, except for those with no forest (farmland and villages), to ensure that multiple camera trapping sites were in each leopard territory (Fig. 1). Beginning in 2016, we established 52 camera trapping sites in the study area, covering 784 km2 (3 cells had 2 camera trapping sites) (Fig. 1). Cameras (Ltl Acorn 6210M, Shenzhen, China) were mounted on trees at a height of approximately 0.4–0.8 m off the ground, deployed along ridges, forest roads, and trails commonly used by leopards, and set to be active for 24 h a day, with a 1‐min delay between consecutive videos. The cameras were programmed to record the time and date when triggered. Each camera worked all year round and was visited 3–5 times to download videos and check batteries. A total of 78% (39 camera sites) of stations had paired cameras (Wang et al. 2017).

The cameras recorded not only leopards but also wild prey species, domestic livestock (cattle, horse, and goat), and human activity. We identified individual leopards based on their unique spot patterns (Karanth & Nichols 1998; Wang et al. 2017), and determined the sex based on the presence or absence of testicles. Consecutive videos of the same species within 0.5 h of each other were not included in the data analysis to avoid inflated counts caused by repeated detections of the same event (O'Brien et al. 2003). We also calculated the relative abundance index (RAI) for each species at each camera trapping site as the number of detections/100 trap days. Due to the difference in management between the northern study area and southern study area, we used the Mann–Whitney u statistics to detect significant differences in the RAI of each species between the 2 areas.

Estimation of density

Two approaches of spatially explicit capture–recapture (SECR), the maximum likelihood method (Efford 2004; Borchers & Efford 2008) and the Bayesian Markov chain Monte Carlo (MCMC) method (Royle & Young 2008; Royle et al. 2009), were used to estimate the density of the North Chinese leopard. SECR has been suggested to be advantageous over other methods for estimating animal density, especially for small data sets typical of studies of large and elusive carnivores (Royle et al. 2009; Sollmann et al. 2011). It is a hierarchical model that contains explicit models of the spatial locations of individuals and their movements to account for imperfect detections (Efford 2004; Royle & Young 2008). To meet the assumption of a closed population and to minimize changes in the activity centers of individuals within the trapping period (Karanth & Nichols 1998), leopard density was estimated in two 100‐day windows (from July to October 2016 and from February to May 2017) consisting of 10 occasions of 10 days of consecutive trappings. The leopard cubs were omitted from the density analysis due to a high mortality rate (Athreya et al. 2013; Wang et al. 2017). The spatial detection history was constructed based on where animals were photographed on a particular occasion (Wang et al. 2017). In both the maximum likelihood and Bayesian frameworks, we included the sex of leopards as a covariate in the models because of sex‐specific differences in the encounter rate and home range size of elusive carnivores (Goodrich et al. 2010; Sollmann et al. 2011; Efford & Mowat 2014).

A maximum likelihood‐based SECR model was fitted to estimate density using the secr package (Efford 2015) in the R environment (version 3.3.2) (Team 2017). Camera traps were treated as proximity detectors that allowed for repeated detections of each individual at a particular trap on any occasion (Efford et al. 2009; Wang et al. 2017). We fitted the density models in the secr package using full likelihood with a half‐normal hazard function (Efford 2015; Xiao et al. 2016). The Akaike information criterion (AIC) and AIC weights were used to rank candidate models, and the models with a ∆AIC < 2 were treated as competing models (Burnham & Anderson 2003). We also used the secr package to conduct a closure test (Otis et al. 1978), as well as calculate the root pooled spatial variance (RPSV) to determine the buffer width (we used 4 times the RPSV as the buffer width) (Efford 2004). The buffer area was defined in ArcGIS 10.1 as a fine mesh of equally spaced grids (here, 1 × 1 km), representing home range centers of all individuals in the survey area, from which we excluded any non‐forest cells and cells with centers of more than the buffer width (Hebblewhite et al. 2011).

The SCRbayes package was used for Bayesian estimation of density in an SECR framework (Royle et al. 2015). Data augmentation of this method allows the maximum potential population size N within the state space (Royle & Young 2008; Royle et al. 2009), which we set as 300 individuals in this study. We ran a half‐normal model for 100 000 MCMC iterations, with a burn‐in period of 20 000 and a thinning rate of 20. For comparison, the state space area was the same as that in the maximum likelihood‐based SECR model. The model's goodness‐of‐fit was assessed by the Bayesian P‐value, with 0.05 < P < 0.95 showing model adequacy (Royle et al. 2013).

Seasonal occupancy models

We assessed seasonal ecological correlates and human factors that influenced the distribution and abundance of the North Chinese leopards across our study area using single‐season occupancy models for both seasons (growing season: July 2016–October 2016; non‐growing season: December 2016–March 2017) (Mackenzie et al. 2002). Occupancy models account for imperfect detections and use presence–absence camera‐trapping data from repeat surveys to estimate the probabilities of occupancy (ψ) and detection (P). Given the territorial behavior of leopards (Mizutani & Jewell 1998) and the short gap between seasons (only 73 days), we chose to run a single‐season occupancy model for each season rather than a multi‐season occupancy model due to the meaninglessness of estimating the local extinction and colonization parameters in such a short sampling period. Given the assumption of occupancy (Mackenzie et al. 2002), detection at a site should be independent from detections by adjacent camera traps. However, the average distance between adjacent camera trapping sites was 2.61 (SE = 0.83) − 2.63 (SE = 0.85) km, spatial autocorrelation between camera trapping sites may affect the result of the occupancy model. Thus, to solve the possible violations of spatial autocorrelation, we employed a hierarchical spatial occupancy model that used a probit mixture framework and a reduced‐dimensional spatial process to improve algorithm convergence (Johnson et al. 2013).

To meet the assumption of demographic closure demanded by occupancy models (Mackenzie et al. 2002), two 100‐day windows, each composed of 10 occasions of 10 days of consecutive trapping during each season, were used. Given the lack of research on the habitat use of the North Chinese leopards, we explored a list of covariates (Table S1, Supporting Information) that may influence their habitat use or behavior, based on previous studies on leopards (Ngoprasert et al. 2007; Balme et al. 2007; Simcharoen et al. 2008; Sugimoto et al. 2016; Wang et al. 2016, 2017). Thirteen variables and 4 variables were considered predictors of leopard probability of habitat use and detection, respectively (Table S1, Supporting Information). The forest type was determined from field sampling. The elevation and TPI were derived from the Shuttle Radar Topography Mission (SRTM) 30‐m digital elevation model, and the TPI was calculated for each camera trapping site using a circular neighborhood with a 1‐km radius (De Reu et al. 2013). We used the RAI as the abundance of each prey species. Because leopards consumed a more diverse range of prey species (Hayward et al. 2006; Sugimoto et al. 2016), and due to the high relative abundance of small mammals in our study area (badger and hare) (Tables 1 and 5), we also included small mammals in our occupancy models. The RAIs of human activity and livestock grazing were also calculated from the camera trapping data as measures of human disturbance. Due to the differences in management between the ZNR and QNR, we designated the study area as a categorical variable in the habitat use model. The interactions between study area and prey species were included in our model as fixed explanatory variables. Oil operations are very common in the QNR, and many oil wells and roads for oil carriers are distributed throughout QNR, we calculated the number of oil wells within a 2‐km radius of each camera trapping site and the road density in each grid (m/km2 in each 4 × 4 km cell). We also derived the distance between each camera trap and the nearest road in ArcGIS 10.1. Road type was recorded during camera trap deployment. We standardized all the continuous variables to z‐scores to improve model convergence. We tested for collinearity among all the continuous variables with the variance inflation factor (VIF). When the VIF > 3 for a given model, one of the covariates was excluded from the model.

Table 1

Summary of wildlife, human activity, and grazing captured by the camera traps, showing the number of independent detections (N), RAI (mean ± SE), and number and proportion of camera traps where the species were captured in the northern study area (ZNR) and southern study area (QS) in central China

Common name	Northern study area (ZNR)					Southern study area (QS)
	N	% of all captures	RAI ± SE	No. of trap sites	% of all trap sites	N	% of all captures	RAI ± SE	No. of trap sites	% of all trap sites
North Chinese leopard	139	3.47	1.65 ± 0.37	21	84	120	2.23	1.46 ± 0.52	17	68
Small‐size mammals	1718	42.93	21.05 ± 3.02	24	96	3115	57.89	39.65 ± 7.51	25	100
Wild boar	902	22.54	11.01 ± 1.71	24	96	625	11.61	7.81 ± 1.02	25	100
Roe deer	974	24.34	11.71 ± 1.53	25	100	1044	18.66	13.23 ± 1.61	25	100
Human activity	229	7.72	2.81 ± 1.29	18	72	395	7.34	4.64 ± 1.79	17	68
Livestock grazing	40	1.00	0.47 ± 0.26	3	12	86	1.60	0.95 ± 0.38	8	32
Total	4002			25		5381			25

Bold face indicates significant differences (Mann–Whitney u test, P < 0.05).

Table 5

Seasonal parameter estimates and 95% credible intervals (CIs) from spatial occupancy models for the North Chinese leopard in central China

	Covariate	Mean	SD	95% CI	Z score
Growing season
Detection	(Intercept)	−1.70	0.19	(−2.00, −1.36)	0.10
	Road type (main road)	0.92	0.27	(0.46, 1.35)	1.90
	Road type (valley)	0.32	0.30	(−0.16, 0.83)	0.55
	Road type (ridge)	1.59	0.30	(1.07, 2.07)	0.40
Habitat use	Intercept	0.47	1.17	(−1,38, 2.32)	−0.35
	Elevation	1.86	1.09	(0.08, 3.59)	1.78
	TPI	−3.58	1.34	(−5.58, −1.42)	−1.08
	Small‐size mammals	2.06	1.22	(0.13, 4.03)	0.39
	Wild boar	0.30	1.32	(−1.81, 2.41)	0.40
	ZNR	2.04	1.47	(−0.28, 4.47)	−1.37
	ZNR × Wild boar	−3.83	1.76	(−6.53, −0.85)	−0.29
Non‐growing season
Detection	(Intercept)	−1.21	0.28	(−1.66, −0.73)	0.56
	Forest type (broad‐leaved forest)	0.37	0.26	(−0.05, 0.79)	−0.33
	Road type (main road)	1.04	0.31	(0.54, 1.54)	−0.39
	Road type (valley)	−0.1	0.27	(−0.54, 0.35)	−0.01
	Road type (ridge)	−0.26	0.24	(−0.65, 0.12)	0.46
Habitat use	(Intercept)	3.39	−1.21	(1.44, 5.34)	−0.35
	Wild boar	3.66	−1.7	(0.83, 6.34)	−0.78
	Livestock	1.33	−0.91	(0.06, 2.57)	0.27
	Forest type (broad‐leaved forest)	−2.08	−1.13	(−3.87, −0.25)	−0.09
	ZNR	0.77	−1.06	(−0.96, 2.42)	0.50
	ZNR × Wild boar	−1.11	−2.04	(−4.44, 2.25)	−0.44

Estimates of coefficients are reported for standardized covariates, scaled to the mean and standard deviation (SE). Bold face indicates that covariates had a significant association with leopard habitat use and detection because their 95% CIs did not overlap with zero. |Z| < 1.96 indicates model convergence. ZNR, Ziwuling Nature Reserve.

Leopards’ detection histories for each camera trapping site were developed based on the records of the cameras, where “1” indicated that the leopard was detected at a specific camera‐trap station on a specific occasion, and “0” represented no detection. We ran models in 2 phases: (i) models without spatial autocorrelation to select covariates and (ii) models with spatial autocorrelation. The occupancy models for the two seasons were modeled separately. The first phase was performed under a maximum likelihood framework and conducted with the R package unmarked version 0.12‐2 (Fiske et al. 2017). First, we modeled the detection probability (P) by using all combinations of covariates (Table S1, Supporting Information) while holding ψ constant, and the top models with a ∆AIC < 2 were considered to contribute equally and used to model probability of habitat use (Lele et al. 2013)in relation to the site covariates (Tan et al. 2017). Then, we ran each combination of site‐specific covariates while the detection probabilities were modeled following models selected previously. The top models with a ∆AIC < 2 were considered to contribute equally. We selected only covariates with summed model weights > 0.5 (Kalies et al. 2012) to model spatial autocorrelation in the next phase. We assessed the goodness‐of‐fit of the global model to evaluate the probability that the model would be correct (P > 0.5) and the accuracy of estimation determined by c‐hat (MacKenzie & Bailey 2004).

In the next phase, the R package stocc version 1.30 (Johnson et al. 2013) was used to model spatial autocorrelation with restricted spatial regression (RSR). The posterior predictive loss criterion (PPLC) was used to compare models without spatial autocorrelation parameters (top‐ranked models from the first phase) to models with spatial autocorrelation (Bayesian RSR) (Gelfand & Ghosh 1998). Due to the lack of home range information for the North Chinese leopard, we set the distance threshold for detecting spatial structure in neighboring sample locations to 4.33 km for the growing season and 4.7 km for the non‐growing season based on the result from leopard density estimation (see Table 3). The moran.cut cut‐off parameter was set to 4.7 for the growing season and 4.1 for the non‐growing season (equal to 10% of the number of camera‐trap sites) (Hughes & Haran 2013). For both the spatial and non‐spatial models, we set flat prior distributions for p and ψ and a gamma (0.5, 0.00005) distribution for the spatial component (Johnson et al. 2013). We ran the Gibbs sampler for 10000 iterations, with a burn‐in of 1000 iterations, to estimate the parameter mean, standard deviation (SD), and 95% Bayesian credible interval (CI). Covariates with a 95% CI that did not overlap with 0 were considered to have a significant association with leopard detection and habitat use. We used Geweke diagnostic statistics (Geweke 1992) to assess model convergence (|Z| < 1.96).

Table 3

Population size and density of the North Chinese leopard in central China from 2016 to 2017 from spatially explicit capture‐recapture models (maximum likelihood method and Bayesian method with an MCMC algorithm)

Method	Parameter	2016			2017
		Mean	SE	95% CI	Mean	SE	95% CI
Maximum likelihood	Density	1.70	0.48	0.98–2.93	1.24	0.40	0.67–2.28
	N	75	21	44–130	55	17	30–101
	λ 0	0.23	0.08	0.12–0.44	0.11	0.03	0.06–0.18
	σ males	2.96	0.51	2.12–4.13	10.86	1.78	7.89–14.95
	σ females	1.76	0.25	1.33–2.33	2.52	0.37	1.89–3.36
Bayesian with MCMC algorithm	Density	2.00	0.53	1.11–3.2	2.40	0.67	1.3–3.9
	N	88	23	50–142	110	30	59–174
	λ 0	0.18	0.08	0.15–0.46	0.08	0.02	0.05–0.12
	σ males	3.02	0.51	2.20–4.15	8.90	1.51	6.57–12.36
	σ females	2.00	0.28	1.53–2.64	2.30	0.31	1.79–3.01

Density is calculated as animals/100 km2; N is the population size of the leopard; λ 0 is the expected encounter rate; σ is the spatial scale parameter; MCMC indicates Markov chain Monte Carlo; SE is the standard error; and 95% CI is the 95% credible interval.

RESULTS

Abundance

From June 2016 to May 2017, 50 camera traps (two camera traps were damaged by livestock or interfered with by humans) continuously operated on 15,914 days (318.28 ± 9.38 trap days/camera). Twenty‐seven individual leopards (12 females, 7 males, 3 of unidentified, and 5 cubs) were captured for a total of 259 independent detections triggered at 38 camera trapping sites (Table 1).

The RAI of small mammals in the southern study area was significantly higher than that in the northern area (P = 0.043) (Table 1). The proportion of independent detections of wild boars in the northern study area was 1.94 times higher than that in the southern study area (Table 1). All camera traps recorded roe deer, but the RAI of roe deer in the south was higher than that in the north (Table 1). The RAIs of human activity and livestock grazing in the southern area were higher than those in the northern area, but this difference was not significant (P = 0.425 and P = 0.161, respectively) (Table 1).

Density of leopards

Two 100‐day sampling periods resulted in 54 and 82 independent detections of leopards in 2016 and 2017, respectively (Table 2). In 2016, 52 independent detections were used to estimate the density of leopards (excluding 2 videos of 1 cub and 1 individual of unidentified), and the detections included 15 individual adult leopards (10 females and 5 males). In 2017, 72 independent detections were used to estimate the density of leopards (excluding 6 videos of cubs and 4 videos of individuals of unidentified), and 9 females and 7 males were recorded by camera traps (Table S2, Supporting Information).

Table 2

List of North Chinese leopards recorded by the camera traps in 2016 and 2017, showing the independent detections (N), relative abundance index (RAI), and number of camera traps where the individuals were captured in the camera trapping study area

Common name	North Chinese leopard
	2016	2017
N	54	82
Work days	4538	4563
RAI	1.19	1.80
No. of camera traps	20	26
% of all camera traps	40%	53.1%

N is the number of independent events.

Because 13 adult individuals were captured in both years, indicating no seasonal shift, and each individual was recorded on average 3.47 (SE = 0.50) times at 1.87 (SE = 0.26) different sites in 2016, which was lower than the 4.50 (SE = 1.19) times at 2.5 (SE = 0.51) different sites in 2017, we calculated the RPSV (7600 m) for 2017, which resulted in a 30‐km buffer width. The closure test calculated with the secr package supported the assumption of a closed population for both 2016 (z = –0.61, P = 0.27) and 2017 (z = –0.70, P = 0.24).

Model selection according to AIC and ∆AIC values indicated that the spatial scale parameters (σ) were influenced by the sex in both 2016 and 2017 (Table S2, Supporting Information). The σ was 1.51–1.68 times and 3.87–4.30 times larger for male leopards than for female leopards in 2016 and in 2017, respectively (Table 3). The maximum likelihood SECR models estimated the mean leopard density to be 1.70 (SE = 0.48) adult individuals/100 km2 in 2016 and 1.24 (SE = 0.40) adult individuals/100 km2 in 2017, with an expected leopard population size of 75 (SE = 21) and 55 (SE = 17), respectively (Table 3). The leopard densities from the Bayesian SECR were 2.00 (SE = 0.53) adult individuals/100 km2 in 2016 and 2.40 (SE = 0.67) adult individuals/100 km2 in 2017, with a population size of 88 (SE = 23) in 2016 and 110 (SE = 30) in 2017 (Table 3). Most home range centers of captured individuals were in the study area covered by camera traps (Fig. 2). The Bayesian P‐values (0.57 in 2016 and 0.58 in 2017) implied model adequacy.

Figure 2

The mean posterior density and predicted home range center of individuals captured in 2016 (left) and 2017 (right) from the Bayesian SECR model.

Estimation of seasonal occupancy models

A total of 47 camera traps continuously worked on 4342 days during the growing season, and 41 camera traps sustained operation on 4100 days during the non‐growing season. We excluded 5 and 11 camera traps for the growing season and non‐growing season, respectively, due to interference from humans or livestock, making sure that only one camera trapping site was included per grid.

The mean RAI of prey species during the growing season was significantly higher than that during the non‐growing season (P < 0.05). Both human activity and livestock grazing were detected in fewer than 50% of sites over the 2 seasons (Table 4), and neither of them differed in frequency between seasons (P = 0.54 and 0.18 for human activity and livestock grazing, respectively) (Table 4).

Table 4

Summary of wildlife, human activity and grazing captured by the camera traps, showing the number of independent detections (N), RAI (mean ± SE), and number and proportion of camera traps where the species were captured during the growing season and non‐growing season in central China

Common name	Growing season					Non‐growing season
	N	% of all captures	RAI	No. of trap sites	% of all trap sites	N	% of all captures	RAI	No. of trap sites	% of all trap sites
North China leopard	51	1.61	1.11 ± 0.33	18	38	71	5.71	1.73 ± 0.43	25	61
Small‐size mammals	1371	43.37	31.23 ± 5.04	47	100	644	51.77	15.71 ± 2.84	38	93
Wild boar	645	20.40	14.40 ± 1.89	45	96	212	17.04	5.17 ± 1.19	33	80
Roe deer	865	27.36	20.26 ± 2.37	46	98	169	13.59	4.12 ± 0.76	29	71
Human activity	197	6.23	4.33 ± 1.59	20	43	109	8.76	2.66 ± 0.97	16	39
Livestock grazing	32	1.01	0.68 ± 0.31	6	13	39	3.14	0.95 ± 0.36	10	24
Total	3161	100		47		1244	100		41

Bold face indicates significant differences (Mann–Whitney u test, P < 0.05).

We removed road density (VIF = 6.34) for the growing season and road density and TPI (VIF = 6.34 and 3.96, respectively) for the non‐growing season. Under a maximum likelihood framework, the top‐ranked detection models (∆AIC < 2), namely, ψ (.) P (effort + road type) and ψ (.) P (road type) for the growing season and ψ (.) P (road type) and ψ (.) P (forest type + road type) for the non‐growing season, were used in the seasonal habitat use analyses. As results, there were 4 and 1 top models account for habitat use variables (∆AIC < 2), respectively for the growing season and non‐growing season (AIC < 2, Table S4, Supporting Information). According to Kalies et al. (2012) and Tan et al. (2017), covariates with summed model weight > 0.5 were used to assess spatial autocorrelation (Table S5, Supporting Information).

The result of the RSR model indicated that the random spatial effect should be considered in parameter estimation (PPLC: 62.607 vs. 62.696 for the growing season; 93.744 vs. 93.872 for the non‐growing season, respectively). During the growing season, leopards significantly preferred sites at a high elevation and with a flat area (TPI: mean = –3.58 (SE = 0.20)) (Table 5), they preferred the ZWL Nature Reserve, and their occurrence showed a positive association with small‐sized mammals and wild boar relative abundance (Table 5). The interaction between the ZWL Nature Reserve and wild boars had a significant negative effect on leopard occurrence (Table 5). During the non‐growing season, the increased relative abundance of wild boars and livestock grazing had a significant positive effect on the occurrence of leopards (Table 5). Furthermore, leopard habitat use strongly decreased in the broad‐leaved forest (Table 5). The z‐scores of all covariates ranged between –1.96 and 1.96 and showed good convergence (Table 5). The average probability of habitat use during the growing season (ψ = 0.64 ± 0.06) was significantly lower than the probability during the non‐growing season (ψ = 0.74 ± 0.06) (P < 0.01). The map based on the RSR models showed a lower probability of habitat use concentrated in the QS nature reserve (south study area) (Fig. 3).

Figure 3

Predicted probability of habitat use (ψ) and standard errors (SEs) for the North Chinese leopard during the growing season and non‐growing season in Shaanxi Province in central China estimated by a restricted spatial regression (RSR) model.

The leopard detection probabilities across the study area were 0.19 (SE = 0.02) and 0.22 (SE = 0.02) during the growing season and non‐growing season, respectively. The detection probabilities of leopards for both seasons were influenced by the road type (Table 5).

DISCUSSION

Population size and density of leopards

Robust estimates of the population size and density of large predators are essential for guiding conservation decisions. We established a monitoring program based on camera traps to study the population size and density of North Chinese leopards in Shaanxi Province, central China. Our estimated densities of the North Chinese leopards ranging from 1.70–2.00 leopards/100 km2 in 2016 and 1.24–2.40 leopards/100 km2 in 2017 based on different estimation methods (maximum likelihood and Bayesian methods) (Table 3) were lower than the leopard densities previously estimated in southern Asia (2.9–13.17 leopards/100 km2) (Kalle et al. 2011; Gray & Prum 2012; Borah et al. 2014; Thapa et al. 2014), but much higher than the densities of the Amur leopard estimated in north‐eastern China (0.30–0.42 leopards/100 km2) (Wang et al. 2017). The population size results based on SECR models estimated 75–88 individuals in 2016 and 55–110 individuals in 2017 (Table 3). The population size and density of the North Chinese leopard obtained in our study are higher than those estimated by a previous study (Alice et al. 2015). In this study area, the higher relative abundance of prey species and lower human disturbance may contribute to the high density of North Chinese leopards, as suggested in previous studies on large felids (Karanth et al. 2004; Ngoprasert et al. 2007; Butler et al. 2013; Steinmetz et al. 2013; Wang et al. 2016, 2017, 2018).

It is difficult to determine whether the difference in population estimates is related to differences in the methods used, but our result may represent a real increase in leopard number in our study area, which is also suggested by the higher number of leopard cubs we recorded in 2017 compared to 2016 (Tables 2 and 4). In both years, the population size and density of leopards based on the Bayesian MCMC method were larger than those based on the maximum likelihood method, especially in 2017 (approximately 2 times larger than in 2016) (Table 3). However, the 95% CIs overlapped, which indicated that the results of the 2 methods were not significantly different. Compared to the maximum likelihood method, the Bayesian MCMC method is preferred for estimation when capture–recapture data sets are small (Royle et al. 2009; Sollmann et al. 2011). Given that female leopards were captured 1.80 (SE = 0.29) times in 2016 and 1.78 (SE = 0.36) times in 2017, and that male leopards were captured 2.00 (SE = 0.55) times in 2016 and 3.43 (SE = 1.00) times in 2017, we employed 2 methods of SECR to estimate the density of leopards with more acceptable levels of precision, which could then be used to guide conservation and management actions (Gerber et al. 2014). According to Mizutani and Jewell (1998), the home ranges of female leopards do not overlap, but the home ranges of male leopards overlap with those of females. We found that the high‐density area in 2017 was more uniform and larger than that in 2016 (Fig. 2). We speculated that the recruitment of leopard individuals (F13, M03, and M09; Table S2, Supporting Information) contributed to the more uniform and larger area in 2017. Our results also showed that the area of the female leopards’ home range centers was relatively stable. Male leopards first choose an area occupied by females, and when other transient male leopards compete with resident male leopards, the home range changes (Mizutani & Jewell 1998). However, we cannot confirm that M03 and M09 occupied a new home range, and we will continue to study these individuals and confirm this speculation.

We found that the relative abundance of prey species (small mammals, roe deer, and wild boar) of leopards during the non‐growing season was significantly lower than during the growing season, while the opposite was true for the relative abundance of leopards (Table 4). The increasing leopards during the non‐growing season suggested that the population of leopards is growing, which was proven by the density estimate. When the temperature is below a certain threshold during non‐growing seasons, small mammals hibernate or sleep, especially badgers (Zalewski et al. 2007). In our study, the number of independent detections of badgers (746) during the growing season was approximately 21 times higher than during the non‐growing season (35) (Table 4). We speculated that the decreasing number of badgers was the reason why the relative abundance of small mammals was significantly lower during the non‐growing season compared to the growing season. The RAI of roe deer during the growing season was approximately 5 times greater than during the non‐growing season (Table 4). During the non‐growing season, roe deer tend to select low‐density, mixed broad‐leaved shrub forest and birch tree‐dominated forest (Jiang et al. 1996) and prefer valleys, but only 9 camera trap sites were located in such areas. Therefore, the low probability of detection of roe deer was likely the main reason for the lower relative abundance of roe deer during the non‐growing season. Similar to roe deer, the RAI of the wild boars during the non‐growing season was approximately 3 times lower than during the growing season. Wild boars prefer grassy areas, and the low temperature and snow cover, by reducing the activity of wild boar, also contributed to the low RAI of this species (Honda 2009).

Prey richness and availability are the most important factors sustaining leopards (Wang et al. 2017). The wild boars and livestock (goat) grazing exhibited a significant positive relationship with leopard presence at fine spatial scales during the non‐growing season, while these relationships did not hold during the growing season (Table 5). Earlier studies suggested that leopards primarily feed on small‐ to medium‐sized species (10–40 kg) (Hayward et al. 2006); thus, specializing on roe deer or goats instead of wild boars is a more appropriate strategy for leopards due to the suitable size of roe deer and goats. However, we did not find a positive relationship between leopards and roe deer (Table 5), perhaps because the roe deer is distributed throughout our study area and has a high abundance in China and, as a result, has not become the limiting factor for the distribution of North Chinese leopards (Table 5 and Table S1, Supporting Information). Adult wild boars, particularly males, are probably too dangerous for leopards to capture (Andheria et al. 2007; Sugimoto et al. 2016; Yang et al. 2018a). The similar habitat selections of leopard and wild boars contribute to the positive relationship between these two species (Wang et al. 2017, 2018) (Table 5). According to the positive relationships between leopards and small animals, the relatively high diversity, and predominant biomass contribution of small animals in leopard diets in other regions (Sugimoto et al. 2016), we speculated that the small animals may be the main prey species of leopards in this area.

The occupancy models suggested that leopards select habitats at high elevations and with a low TPI during the growing season, but we did not find evidence of such selection during the non‐growing season (Table 5). The results of the occupancy models also indicated that leopards tended to select dirt roads and ridge (Table 5). Leopards were previously reported to select high‐elevation areas and ridge trails (Wang et al. 2017), but we are the first to find that leopards tend to select dirt roads. We speculated that the spatial niche breadth of the North Chinese leopard, which is the apex predator in our study area, shifted because the leopard does not coexist with any other large predators. Leopards might select dirt roads and areas with a low TPI because they provide corridors facilitating travel and marking by leopards. The other reason for the selection of dirt roads could be that human disturbance is not extensive in our study area, especially in the ZWL Nature Reserve located in the northern area (Tables 1 and 5; Figs 1 and 3). Although the oil operations (oil wells and road density) do not seem to have a negative effect on leopards based on the occupancy models, the result from the RSR model indicated that the probability of habitat use was relatively low near high‐density oil wells and negatively associated with road density.

Conservation implications and recommendations

The most important step in protecting the North Chinese leopards is determining how many leopards are currently distributed in China and where they are distributed. Although Kitchener et al. (2017) revised the taxonomy of leopards and combined the North Chinese leopard with the Amur leopard due to the obscure biogeographical barrier between them, given the lack of more accurate molecular evidence, this study still considers the North Chinese leopard to be an independent subspecies. The protection of the North Chinese leopards requires a combination of scientific field surveys and molecular research to perform detailed investigations of the distribution area of the North Chinese leopards in China. Such research requires distributional data and DNA samples of the 4 subspecies of leopards throughout China to explore these subspecies and their population status. Therefore, the importance of the North Chinese leopard needs to be recognized, and we recommend that the government increases protection efforts and funding for leopard conservation. We also recommend that researchers increase the scope of long‐term field surveys, and strengthen interdisciplinary researches on population, habitat, prey species, reproduction, behaviors, and genetic analysis, to better understand and make policy to protect leopards.

CONCLUSION

In this study, we provide the first comprehensive evidence documenting the status of North Chinese leopards in central China as well as their habitat selection during different seasons. The population of North Chinese leopards is increasing in our study area, and the density of leopards is higher in this area than in other areas of China. Prey species are very important in determining the habitat use of leopards, but human disturbances, especially oil wells, seem to have negative impacts on the habitat use of leopards. Protection of the North Chinese leopard will require joint efforts by the government and researchers to prevent the species from becoming extinct.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

Table S1 Variables considered for occupancy models to predict habitat use of North Chinese leopard

Table S2 Summary of North Chinese leopard by camera traps in 2016 and 2017, showing ID number, gender, the number of camera traps, capture frequencies and capture history in 10 sampling occasions

Table S3 Model selection of the maximum likelihood‐based SECR models for the combined‐sex estimation of the North Chinese leopard population

Table S4 Top likelihood‐based occupancy (habitat use) models (ΔAIC ≤ 2) for the North Chinese leopard, with p (effort + road type) and p (road type) for the growing season and p (road type) and p (forest type + road type) for the non‐growing season, in central China

Table S5 Summary of summed model weights for parameters derived from the Top likelihood‐based occupancy (habitat use) models, see Table S3

Click here for additional data file.