Identifying important conservation areas for the clouded leopard Neofelis nebulosa in a mountainous landscape: Inference from spatial modeling techniques

Abstract The survival of large carnivores is increasingly precarious due to extensive human development that causes the habitat loss and fragmentation. Habitat selection is influenced by anthropogenic as well as environmental factors, and understanding these relationships is important for conservation management. We assessed the environmental and anthropogenic variables that influence site use of clouded leopard Neofelis nebulosa in Bhutan, estimated their population density, and used the results to predict the species’ site use across Bhutan. We used a large camera‐trap dataset from the national tiger survey to estimate for clouded leopards, for the first time in Bhutan, (1) population density using spatially explicit capture–recapture models and (2) site‐use probability using occupancy models accounting for spatial autocorrelation. Population density was estimated at D^Bayesian=0.40 (0.10 SD) and D^maximum−likelihood=0.30 (0.12 SE) per 100 km2. Clouded leopard site use was positively associated with forest cover and distance to river while negatively associated with elevation. Mean site‐use probability (from the Bayesian spatial model) was ψ^spatial=0.448 (0.076 SD). When spatial autocorrelation was ignored, the probability of site use was overestimated, ψ^nonspatial=0.826 (0.066 SD). Predictive mapping allowed us to identify important conservation areas and priority habitats to sustain the future of these elusive, ambassador felids and associated guilds. Multiple sites in the south, many of them outside of protected areas, were identified as habitats suitable for this species, adding evidence to conservation planning for clouded leopards in continental South Asia.

important element of the landscape approach to carnivore conservation (Chapron et al., 2014). For example, proportion of area occupied (PAO) and extent of occurrence (EOO) are used extensively to calculate species distributions in the IUCN Red List (IUCN, 2001), thus enabling the identification of important areas for conservation management and protection. In addition, estimating reliable population density is critical in animal ecology as inferences on the population dynamics of species that naturally occur in low densities (such as most felids) are central to their effective conservation and management (Karanth, Nichols, Kumar, & Hines, 2006;Kéry, Gardner, Stoeckle, Weber, & Royle, 2010;Obbard, Howe, & Kyle, 2010).
Density and habitat-use information are also important in assessing risk under climate, environmental, and land-use scenarios (Dormann et al., 2007).
The clouded leopard Neofelis nebulosa is the smallest of the large felids (Austin, Tewes, Grassman, & Silvy, 2007; Grassman et al., 2016). The species is not well known, and the rapidity of environmental change in Bhutan adds urgency to acquiring ecological data to inform its conservation (Penjor, 2016).
Illegal persecution and trade in skins and parts threaten clouded leopard survival (Nowell & Jackson, 1996;Nijman & Shepherd, 2015;Min et al., 2018). Rapid development, increasing human population, deforestation and habitat fragmentation, and ineffective protected area management add to the pressing threats to the clouded leopard across its range (Choudhury, 1993;D'Cruze & Macdonald, 2015;DeFries, Hansen, Newton, & Hansen, 2005;Rabinowitz, 1988). Against this background, a reliable estimate of the population is an important datum for the understanding of forest carnivore guilds in Southeast Asia and also important to underpin conservation planning (Murphy & Noon, 1992).
More broadly, development is rapid in Bhutan and forest loss is accelerating faster than expected (DoFPS, 2011).
We quantified the site use of clouded leopards using spatial occupancy models that account for spatial autocorrelation (SAC hereafter; Johnson, Conn, Hooten, Ray, & Pond, 2013) and for imperfect detection (MacKenzie et al., 2002(MacKenzie et al., , 2006. The spatial occupancy model partitions out the spatial component from the environmental F I G U R E 1 Clouded leopard Neofelis nebulosa captured in one of the camera traps in the study area effects. This approach improves inference when observations are spatially autocorrelated and helps to reduce bias and overestimated precision when observations are not conditionally independent given the covariates (Johnson et al., 2013). Although MacKenzie et al., (2009) state that using appropriate covariates would account for SAC, we decided to include a spatial component in our model because of nonindependence between camera stations (the same individual clouded leopards were "captured" at multiple stations) and correlated variables (Broms, Johnson, Altwegg, & Conquest, 2014;Dormann et al., 2007;Johnson et al., 2013). Density is estimated using both maximum-likelihood (Efford et al., 2009) and Bayesian (Royle, Chandler, Sollmann, & Gardner, 2014) models, enabling us to compare the estimates from the two (Noss et al., 2012). These spatial capture-recapture (SCR) models incorporate spatial information on detection and allow for heterogeneity based on sex, age, and response behavior of animals . These improve reliability of the density estimates when compared to nonspatial models (Efford et al., 2009).
This study used data from the national tiger survey, the first comprehensive nationwide camera-trapping survey conducted to assess the tiger population in Bhutan. We aimed to (1) provide the first population density estimate of the clouded leopard in Bhutan, (2) identify the environmental and anthropogenic variables that influence site use, and (3) predict the site-use intensity across Bhutan and identify important areas for protection. We hypothesized that clouded leopards would prefer forested habitat away from anthropogenic features such as settlements and roads (Austin et al., 2007;Sunarto et al., 2012).

| Study area
Bhutan (27°24′14″N, 90°24′2″E) is a small landlocked country bordered by the Tibetan Autonomous Region (China) in the north and India in the east, west, and south ( Figure 2). The area of the country is 38,394 km 2 with a total population of 779,666 (NSB, 2017). The climate is characterized by four distinct seasons: winter (December-February), spring (March-May), summer (June-August), and autumn (September-November). Rainfall is persistent and heavy during the monsoon (July-September), largely fed by moisture-laden winds from the Bay of Bengal. Precipitation varies between 300 mm (north) and 5,000 mm (south) per year, and the temperature ranges from subzero in the north to above 35°C in the humid subtropical south.
As a result, the vegetation type varies from lush subtropical broadleaved forest to dry alpine meadow as the altitude rises from <100 to >5,000 m a.s.l.

| Field survey
The survey was organized in two major blocks: south and north.
It was conducted between March 2014 and March 2015. The survey was not conducted at settlements, agriculture land, and areas above 4,500 m because we judged these areas as nonhabitat and there was no previous record of clouded leopard reported at these places. After eliminations of these areas, a grid of 1,522 grid cells and 25-km 2 grid space was laid out across the country ( Figure 2; for details of survey area refer to Table 3). In total, 1,129 F I G U R E 2 Location of Bhutan in continental Southeast Asia and camera locations (n north = 681; n south = 448) camera stations were deployed along animal tracks or locations where field signs suggested maximum detection probability. Five different camera models (Bushnell™, CuddeBack™, HCO-ScountGuard™, Reconyx™, and U-Way™) were used for the survey. In each station, two cameras were installed, 45 cm above the ground and at least 5 m apart, facing each other.
A minimum distance of 2 km between each station was maintained whenever possible in the context of impassable rivers, steep slopes, and settlements. The camera traps were deployed for four months in the south (March-June 2014) and 5 months in the north (October 2014-March 2015. Monthly checks were carried out to retrieve data, change batteries, and clear obstacles in front of the camera lenses. We used images from stations that captured both flanks of clouded leopard to estimate their density. At some stations, the clouded leopard was captured at only one camera trap. Such single images cannot be identified as individuals and so were excluded from the density estimation but not the occupancy analysis (Foster & Harmsen, 2012). We used data from the south block only to estimate clouded leopard population density because the sample size was adequate (n south = 448) and the cameras were deployed on a single season (March-June 2014). The number of captures in the north block was few (9 of 681 stations), and while this is a revealing result in itself, the numbers were too low for density estimation ( Figure S2). We used the entire sample to analyze site-use probability (n total = 849).

| Site use
The photographic records of clouded leopard were converted into binary detection histories representing 1 for detection and 0 for nondetection (MacKenzie et al., 2002). To minimize the risk of violating the closure assumption when estimating site use, the capture period used for the analysis was the first 120 days of each camera station's history (Rota, Fletcher, Dorazio, & Betts, 2009).
The 120-day subset was collapsed into sampling periods to increase temporal independence among occasions and overall detection probability (Dillon & Kelly, 2007;Otis, Burnham, White, & Anderson, 1978;Tan et al., 2017). The optimum number of days per occasion was selected based on chi-square goodness-of-fit test for the multivariate global model (MacKenzie & Bailey, 2004). A 15-day period proved optimal and was used for further occupancy analyses (Table S3).
Occupancy ψ is defined as the probability that a species will occupy a random site at a given time period (MacKenzie et al., 2006).
We refer to ψ as the probability of site use (MacKenzie et al., 2006, p. 105), which is the probability of clouded leopard using a sampling unit (camera station; Mohamad et al., 2015). Occupancy models can accommodate covariates, and hence, detection and occupancy probabilities can be modeled as a function of a survey and site-specific covariates (MacKenzie et al., 2006). Site covariates were selected based on literature on clouded leopards (Mohamad et al., 2015;Tan et al., 2017). Site covariates for the whole of Bhutan were processed using QGIS 2.14 (QGIS Development Team, 2016) and ArcGIS 10.2 (ESRI, 2011). Each site covariate value was the mean of raster cells bound within concentric circles of 500 m radius around each camera-trap station. It was calculated using the "zonal statistics" tool in QGIS. This radius distance was chosen to represent the average site characteristics around each camera station. Elevation, slope, and aspect values were extracted from a 30-m resolution digital elevation model raster data (DEM; USGS, 2016). Vegetation data came from two different sources: (1) a 250-m resolution vegetation continuous field (VCF) which shows the percentage of pixels covered by vegetation above 5 m height (DiMiceli et al., 2011) and (2) a 30-m resolution global forest change (GFC) cover for the year 2014 (the camera traps were deployed in that year) (Hansen et al., 2013). With the GFC layer, users can define the percentage of tree cover to be considered as forest. GFC thresholds of 30%, 50%, 75%, and 90% tree cover were tested independently for their effects on occupancy via univariate modeling (Table S2). The different threshold levels were chosen to avoid subjective selection as we did not have supplementary data or field verification (Tan et al., 2017;Tropek et al., 2013). Distance to major rivers, settlements, roads, protected areas, and commercially logged forest was generated using the Euclidean distance tool in ArcGIS. These vector layers were rasterized to generate the distances in meters.
We ran single-season, single-species occupancy models with the R package "unmarked" (Fiske & Chandler, 2011) to estimate the maximum-likelihood probability of site use, accounting for imperfect detection (MacKenzie et al., 2006). Comparison between all possible models of covariates was made using the R package "AICcmodavg" (Mazerolle, 2015). We used AIC c corrected for small sample size (AIC c ) for model selection (Burnham & Anderson, 2004).
All multivariate models within ΔAIC c < 2 of the best-performing multivariate models were considered to be strongly supported by the data (Burnham & Anderson, 2004). The parameter estimate of each covariate was averaged using the R package "MuMIn" (Barton, 2016).
Detection probability p was modeled as a function of survey covariates viz., survey area and effort (number of active camera-trap days during each sampling occasion). Site-use probability ψ was modeled as a function of site covariates (aspect, elevation, distance to logged forest, GFC and VCF forest cover, distance to river, distance to road, distance to protected area, distance to settlement, and slope). A two-stage modeling approach was adopted to reduce the number of combinations of every possible covariate (Long, Donovan, MacKay, Zielinski, & Buzas, 2011). First, we modeled the detection p(.), while keeping ψ(.) constant by running possible combinations for detection covariates (effort and survey area) additively, that is, p(variable) ψ(.). Probabilities of detection might differ among survey areas, which could be due to various factors such as disturbance level in the survey area or topography (Tan et al., 2017). We then retained significant detection covariates and ran multivariate models for occupancy with site covariates.
All continuous site covariates were z-standardized x −x∕σ to a mean of 0 and unit standard deviation prior to analysis in order to facilitate model convergence and comparisons among covariates (Stanton, Thompson, & Kesler, 2015). A preliminary set of 13 covariates was tested for collinearity using Pearson's correlation in R 3.3.1 (R Core Team, 2016). Any pairwise combination with coefficient |r| ≥ 0.6 was considered correlated (Table S1). Multivariate collinearity among predictor variables is often found to hinder the analysis of ecological processes, thus producing biased estimates (Cade, 2015). Therefore, from the correlated pairs, the covariate that performed better (low AIC c ) in the univariate single-season occupancy models was retained for the multivariate model. The goodness-of-fit test for the most parameterized multivariate model was performed with 1,000 bootstraps (MacKenzie & Bailey, 2004).
To account for SAC in the occupancy models, we used R package "stocc" (Johnson, 2015) to perform the spatial occupancy modeling of site use (Johnson et al., 2013). The best model from the maximumlikelihood analysis was fitted using a Bayesian framework following Broms et al. (2014;Kéry & Schaub, 2012). This method employs restricted spatial regression (RSR) using the probit link formulation (Johnson et al., 2013). RSR improves computational efficiency, minimizes confounding between parameters, and corrects for spatial effects in the covariates (see Johnson et al., 2013;and Broms et al., 2014, for details). The final model describing the site-use probability was where β 0 refers to intercept, β i are the coefficient estimates of the covariates, i is the site surveyed, η i(τ) represents the spatial process for the ith site, and τ ~ Gamma (0.5, 0.0005) is a parameter that controls the spatial process. Following Johnson et al. (2013), the prior for the τ heavily weights large value (implying less spatial autocorrelation). As a result, any spatial effect seen is strong evidence for spatial autocorrelation.
The Moran cut used in the spatial model was 10% of the total number of camera stations (0.1 × 849 = 84.9; Hughes & Haran, 2010). The Moran statistic is a measure of SAC, analogously interpreted as the correlation coefficient for the correlation in site-use probabilities across sites (see Johnson et al., 2013; for details). The threshold of 8.7 km was used. This indicates that all sites within this distance are considered neighbors and this distance was chosen based on the movement parameter derived from spatial capture-recapture analysis (see 3.2 estimates in the Section 3; Table 4).
Each Bayesian model was run with 1 chain for 50,000 iterations.
We discarded the first 10,000 as burn-in and employed a thinning rate of 5.
The mean untransformed beta coefficients and 95% credible intervals (CRI) were used to examine the degree and direction of the covariate effect on the site-use probability. We considered covariates to have a strong influence on the clouded leopard's site use if their 95% CRI did not include 0.
We used coefficient estimates from the Bayesian model corrected for SAC to predict clouded leopard site use in forested pixels across Bhutan. Univariate modeling revealed that GFC (forest cover) was an important variable influencing clouded leopard site use. For each forested pixel (for the year 2014), we extracted values of important covariates using a 90-m resolution. We then predicted site use for these pixels with parameter estimates from the best multivariate model. These pixels were then imported into QGIS to map out clouded leopard site suitability across Bhutan.

| Density
Density was estimated using maximum-likelihood and Bayesian approaches both of which are spatially explicit capture-recapture (SCR) methods (Efford et al., 2009;Royle et al., 2014). We chose to apply both methods to the same dataset for comparison (Noss et al., 2012). Spatial methods address edge effects using temporal and spatial information in addition to detection data. Further, they eliminate ad hoc estimation of the effective survey area (Efford, Borchers, & Byrom, 2009;Noss et al., 2012). The R package "secr" (Efford, 2016) was used to estimate maximum-likelihood density (Efford et al., 2009), and the package "SPACECAP" (Gopalaswamy et al., 2014) was used to estimate density using a Bayesian modeling framework.
It is assumed that animals have circular home ranges with fixed centers . The half-normal detection function was used to model home range, and this is comprised of two parameters: g 0 , the encounter rate at the home range centre, and the scale parameter σ, which describes the decline of the encounter rate from home range centers (Tobler & Powell, 2013).
A spatial mask comprising of a 25-km buffer was created around the outermost camera-trap stations (Soisalo & Cavalcanti, 2006). This buffer distance corresponds to approximately 3× the scale parameter (σ) (8.7 km) and is large enough to contain potential home range centers of all clouded leopards caught on our survey grid. All unsuitable habitats, such as settlements, roads, and large rivers, were assigned 0 as opposed to suitable habitats (assigned 1). Equidistant points in the space mask represented potential activity centers. The spacing between points was entered as 500 m. This is the shortest distance allowable and enables finescale modeling of where actual home range centers might be based on the camera-trap data (Efford & Efford, 2016;Gopalaswamy et al., 2014). For the SPACECAP analysis, the input data were prepared as prescribed in the SPACECAP manual (Gopalaswamy et al., 2014). We executed 60,000 Markov chain Monte Carlo (MCMC) iterations, a burn-in of 10,000, and a thinning rate of 1. The data augmentation value was set to 190 (10 times the number of the total number of identified individuals). Geweke diagnostic scores were used to check convergence of chains (Gopalaswamy et al., 2012). We report posterior means with standard deviations and 95% highest posterior density intervals (HPDI). Due to ambiguity probit(ψ i ) =β 0 +β ele elevation i +β for forest cover i +β riv distance to river i +η i (τ) in photographs, sex-based identification was not possible; thus, we did not report the sex-specific density estimates in both approaches.
A surface density map was generated to depict posterior estimates of pixel densities. The final discrete habitat map consisted of equally spaced grid cells of 0.25 km 2 , for visualizing the fine-scale variation in density across southern Bhutan ( Figure S3).

| RE SULTS
We retrieved camera traps from 849 stations (of 1,129 stations), the rest were lost to animal damage, theft, and malfunction. All stations in the south were recovered (448 of 448), while only 401 of 681 stations were recovered from the northern block.

| Site use
Clouded leopards were detected in 114 of 849 stations totaling an effort of 62,739 trap days and giving the overall naïve occupancy estimate of 0.134. All the GFC forest cover variables were correlated with one another and with VCF. Within each correlated pair, the covariate performing better in the univariate occupancy models were GFC at 90% threshold (GFC90), elevation (ELE), slope (SLO), aspect (ASP), distance to logged forest (LOG), distance to major river (RIV), distance to protected area (PA), distance to road (ROA), and distance to settlement (SET). MacKenzie and Bailey (2004)  The probability of detection was influenced by survey area and effort (number of active camera-trap days per sampling occasion; Table 1). The effort was positively correlated with detection probability β (ŜE) = 0.085 (0.023) and varied across the surveyed areas (Table 3; Figure S4). The model with elevation, forest cover, and distance to river had the highest support from the likelihood-based analysis ( Table 2).
As for the Bayesian analyses, a slightly lower posterior predictive loss criterion was reported for the spatial model as compared to the nonspatial model (PPLC spatial = 316.806; PPLC nonspatial = 317.085). This indicates that the inclusion of the spatial random effect could be necessary. Similarly, for the spatial model, the posterior distribution of the spatial variance parameter (σ = 1∕ √ τ) was not very far from zero (95% CRI: 0.053-34,104.68), implying that additional spatial effect contributed to the variability of site use (Johnson et al., 2013) but not to a large extent. The spatial occupancy model resulted in lower site-use probability and narrower CRI compared to likelihood-based estimates (ψ nonspatial = 0.829 (0.066 SE) versus ψ spatial = 0.448 (0.076 SE); Table 3, Figure 4). There was strong evidence to suggest that site use was negatively associated with elevation and positively associated with forest cover and distance to rivers as the 95% CRIs did not include zero (Table 3). The mean predicted site-use probability in forested areas is 0.448 (0.076 SE). Distance to logged forests and protected areas were negatively associated with site-use probability, while site-use probability was positively associated with distance to road, slope, and aspect. However, these effects were at best weak, with high standard errors and 95% CIs overlapping zero (Table S4).
In estimating site-use probabilities, we found that the mean siteuse probabilities of protected and nonprotected areas were similar (ψ PA = 0.431(0.056 SE); ψ outsidePA = 0.451 (0.060 SE)). Sites inside and outside of protected areas in the southern region were predicted to have higher use probability than sites inside and outside protected areas in the north (Figure 4; Figure S1). Models strongly supported by the data (ΔAIC c < 2) are shown. Site covariates tested: elevation (ELE), global forest change of 90% tree cover threshold (GFC90), slope (SLO), aspect (ASP), distance to rivers (RIV), to logged forest (LOG), to road (ROA), to settlement (SET), and to protected area (PA). All models include different survey areas (Survey Area) and the number of active camera days per sampling occasion (Effort) as detection covariates, p(Survey Area + Effort).

TA B L E 2 Multivariate model selection results of clouded leopard site-use probability
TA B L E 3 Parameter estimates, standard errors, and 95% CRIs from the best-fitting model for clouded leopard site use

| Density
We used 321 images and 20 videos from south stations only (n stations = 448) for the year 2014 for SCR analysis. The total effort was 23,249 trap days. The mean spacing between the successful stations was 2,963 m (±1,238 m SD); however, the distance between all stations was not uniform due to terrain (range: 80-16,600 m).
Nineteen individuals were identified. The number of captures and recaptures ranged from one to three for identified individuals (mean ±SD = 1.37 ± 0.58; Figure S2).
Results from the Bayesian analysis are presented in Table 4. The maximum-likelihood density estimate for clouded leopards was 0.30 (±0.12 SE) individuals per 100 km 2 . We ran only three models for the maximum-likelihood analysis (Table S5). Due to computational limits and a small number of captures/recaptures ( Figure S2), complex models of more than one parameter and the time-based model were unable to run. The addition of parameters, however, did not improve the model, and we reported our results from the most parsimonious model, g 0 ~ 1 (Table 4). Capture probability at home range centre (g 0 ) was estimated at 0.30 (±0.10 SE) per 100 km 2 , and the scale parameter σ was 7.59 km (±1.73 km SE). The density estimates from the two approaches were similar (Table 4).
There were pockets of high-density regions in the south ( Figure   S3), most of which were overlapping with high predicted site use ( Figure 4).

| D ISCUSS I ON
Using bycatch from national tiger survey, we have estimated density and site-use probability of a cryptic species, the clouded leopard, following methodologies developed in previous studies (Borah et al., 2014;Mohamad et al., 2015 for N. nebulosa;Sollmann et al., 2014;Haidir et al., 2013 for N. diardi).
The naïve site-use probability of 13.4% was underestimated compared to the overall estimated probability of 44.8%, which was approximately three times the former. This confirmed the need to account for imperfect detection and doing so substantially improved the predictive ability of occupancy models. For a notoriously elusive species, failure to detect them when they are present at a site will induce a common source of survey bias (Linkie, Dinata, Nugroho, & Haidir, 2007). Thus, explicitly estimating detectability is important to estimate spatial distribution reliably. Probabilities of detection differed among survey areas which might be due to behavioral responses (e.g., to disturbance level in the survey areas or to the season; Tan et al., 2017).
Other factors that might influence clouded leopard abundance and/ or detectability include the availability of prey (Mohamad et al., 2015) and/or other large carnivores such as tigers and leopards (Tan et al., 2015). As the focus of this study is on the influence of geographical factors, the effects of interspecific interactions on the detection of  (Table 3, Figure 3a). We suspect that the lowest altitude in Bhutan might be close to the highest in more tropical landscapes; hence, our finding of clouded leopards in lower elevations may not be at odds with those from earlier studies. This also suggests that clouded leopards are less likely to be found at higher elevations in forested habitat at high elevations may be used by clouded leopards.
Elsewhere, studies have concluded that elevation was the main determinant of occupancy (Haidir et al., 2013 for N. diardi), while others concluded it was an indirect factor operating through an effect on prey abundance (Mohamad et al., 2015 for N. nebulosa). An earlier study in Bhutan suggested that elevation was not correlated with the abundance of putative prey of clouded leopards , although knowledge of this felid's diet is poor (Mohamad et al., 2015).
We found that a higher percentage of forest cover in the 500 m buffer area of each camera station favored clouded leopard site use ( Figure 3b). This supports our hypothesis and concurs with a previous study that demonstrated the clouded leopard preference for forest cover in Peninsular Malaysia (Tan et al., 2017). The use of forested habitat over grassland and secondary vegetation by clouded leopard has also been shown in previous studies (e.g., Austin et al., 2007).
Clouded leopards are suspected to be arboreal (Rabinowitz et al., 1987) and forest-obligate species; forest cover is thought to play a vital role in providing concealment while hunting (Nowell & Jackson, 1996). Our findings indicate that it is important to protect large tracts of forest both inside and outside protected areas. Contiguous forest cover over extensive landscapes may facilitate movement of clouded leopards, serving as corridors for dispersal. High forest cover has the potential to provide high-value habitat that could serve as a refuge and food resources within the species' home range (Galvez et al., 2013).
High turnover of site-use probability in corridors and areas beyond the protection offered by protected areas could reflect the importance of forest cover. Forest cover can affect the hunting success of carnivores, providing concealment (Balme, Hunter, & Slotow, 2007).
Forest protection is affected by infrastructure development. Recently, Bhutan has seen the extensive construction of road network to connect all remote villages to district headquarters. In the process, the forest (largely subtropical broadleaved forest) has undergone major loss (Reddy et al., 2016). We suspect such losses are malign indicators of the decrease in important wildlife habitat.
Clouded leopards avoided habitat close to rivers (Figure 3c).
Most settlements are concentrated near major rivers, probably contributing to their unattractiveness for clouded leopards (see Sunarto et al., 2012;and Tan et al., 2017). Our findings are at a grain too coarse to comment on the role of smaller rivers deeper in forests.
Failure to account for SAC resulted in overprediction, evident in the predictive map derived from the nonspatial model ( Figure 4).
Sites to the north and central Bhutan, with human settlement and higher elevation, were predicted to have high site use by the clouded leopard. However, based on field experiences of authors (U.P., S.W., and T.), those sites were unlikely to be used by clouded leopards.
These effects showed that the lack of spatial independence in detection rates resulted in overestimated site-use probability. The clouded leopards were moving between stations on a single occasion (as suggested by the high 8.7 km movement parameter), and hence, the detection rates of trap stations closer to one another are more similar.
The spatial model, when taking into account this similarity, produced a lower site-use probability than the nonspatial model. The large difference in occupancy estimates between spatial and nonspatial models suggests that the spatial structure is prominent; that is, we had tight concentrations of high detection rates versus low detection rates. Additionally, models accounting for SAC revealed residual spatial correlation that better captured the clouded leopard's presence. Hence, when planning space use and protection of forests, the use of the range map accounting for SAC would be more conservative and appropriate.
The density estimates from maximum-likelihood and Bayesian approaches were similar (D maximum−likelihood = 0.30∕100 km 2 and D Bayesian = 0.40∕100 km 2 ). These estimates are comparable to, though lower than, those from produced by SCR models in previous studies (Table S5). Although little is known of clouded leopard spatial organization, female clouded leopards probably have home ranges of about 16-40 km 2 with a core of about 5.4 km 2 (radii of core and home range = 1.3 km and 2.99 km, respectively; Austin et al., 2007;Hearn et al., 2016). The grid of 5 km × 5 km with mean spacing between stations of 2.96 km (±1.24 km SD) can certainly miss such cores, leaving "holes" in the grid (Noss et al., 2012). The camera traps would be placed further apart for a survey grid designed for tigers than for clouded leopards. This spacing would help in minimizing the violation of the assumption of occupancy analysis that the "detection of the species at a site is independent of detecting the species at any other sites." However, for density estimation, this spacing would introduce "holes" and likely yield conservative estimates. Our analyses also do not account for biases such as a possible difference in detecting males and females on and off trails (Balme, Hunter, & Slotow, 2009;Wegge, Pokheral, & Jnawali, 2004).
Such methodological biases are widely recognized and deserve further research (Foster & Harmsen, 2012; see also Sollmann et al., 2011;Broekhuis & Gopalaswamy, 2016). Our study provides the first density estimates of the clouded leopard in Bhutan. This information sets an important benchmark for clouded leopard monitoring programs and underscores the rarity of this species to prompt conservation planning.

| Methodological considerations
Although we successfully used the same dataset to estimate both density and occupancy (site use), there are certain assumptions underlying SCR and occupancy modeling that must be adhered to strictly in order to reduce bias in estimates. SCR requires adequate captures and recaptures at multiple sites, while occupancy, on the other hand, requires spatial independence between sites. Notwithstanding the risk of violating the underlying assumptions, we reconciled these requirements by explicitly modeling occupancy (site use) using the restricted spatial regression (RSR) method incorporating both detection probability and correcting for spatial autocorrelation (see Johnson et al., 2013 for details). The RSR model had two main effects: a decrease in the average probability of site use and narrowing of the width of confidence intervals of the covariates (Table 3). The posterior predictive loss criterion (PPLC) favored the spatial model over the nonspatial model, although the difference was small. The performance of PPLC for hierarchical models is still unknown and therefore should be interpreted with care. It is analogous to Akaike information criterion (AIC) for model selection but does not show whether the model fits the data well (Broms et al., 2014).

| Management implications
A suite of species-habitat relationships for clouded leopard was identified, some of which accords well with the results of previous studies, whereas others are new. For example, we found a negative relationship between site use and elevation, the opposite of previous studies (Tan et al., 2017 for N. nebulosa;Haidir et al., 2013 for N. diardi), probably due to differences in altitude ranges between studies. We recommend, therefore, that future studies on clouded leopards in Bhutan place camera traps at lower altitudes. Our use of bycatch information derived from the large-scale tiger survey illustrates the great potential for additional returns on investment in camera-trap studies (Mohamad et al., 2015;Sollmann et al., 2014).
The predictive mapping indicated that forests inside and outside protected areas were equally suitable for the clouded leopard ( Figure S1). Habitat suitability was high in the low-altitude areas of the southern region (Table 5). Clouded leopards have considerable potential as ambassador species, fostering the protection of sympatric carnivores (Asiatic golden cat Catopuma temminckii, marbled cat Pardofelis marmorata, and leopard cat Prionailurus bengalensis; Singh & Macdonald, 2017). Although growing human footprint is a cause for concern, we did not find evidence that it exerts pressure on the current distribution of clouded leopard. Nevertheless, our finding of an effect of distance to river may indicate an indirect effect of disturbance on habitat selection. Management should seek to minimize loss of forest cover and identify, at a finer scale, the limiting factors influencing site use in particular areas. The predictive spatial occupancy map, together with the density map, also revealed many isolated but seemingly optimal sites for the clouded leopard conservation beyond protected areas ( Figures S1 and S3).
We hope these maps will prove useful for future land-use and conservation planning. Bank, WWF-Bhutan, and Bhutan Foundation for the financial support to conduct the national tiger survey. The authors thank heartily all the national tiger survey field crews for painstakingly collecting the data.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R S' CO NTR I B UTI O N
U.P., D.W.M., and C.K.W.T. conceived the ideas and designed analyses; S.W. designed and coordinated the field survey; T., with support from field staff, collected data; U.P. analyzed the data; U.P. led the writing of the manuscript. D.W.M. and C.K.W.T. provided input on the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

DATA ACCE SS I B I LIT Y
All the data used in the analyses will be made publicly available with the supporting information. Few information (such as coordinates) due to sensitive nature (data from national tiger survey) will not be available on public domain. Bold values indicate that the 95% CI of the mean predicted site use is below and not overlapping the overall average of 0.448. Divisions are nonprotected areas. Refer to Supporting information Figure S1 for site location.