Multi‐scale habitat modelling identifies spatial conservation priorities for mainland clouded leopards (Neofelis nebulosa)

Deforestation is rapidly altering Southeast Asian landscapes, resulting in some of the highest rates of habitat loss worldwide. Among the many species facing declines in this region, clouded leopards rank notably for their ambassadorial potential and capacity to act as powerful levers for broader forest conservation programmes. Thus, identifying core habitat and conservation opportunities are critical for curbing further Neofelis declines and extending umbrella protection for diverse forest biota similarly threatened by widespread habitat loss. Furthermore, a recent comprehensive habitat assessment of Sunda clouded leopards (N. diardi) highlights the lack of such information for the mainland species (N. nebulosa) and facilitates a comparative assessment.


| INTRODUC TI ON
Deforestation, fire, and land conversion (e.g., to large-scale oil palm and Acacia monocultures) are rapidly altering South and Southeast Asian landscapes (Cushman, Macdonald, Landguth, Malhi, & Macdonald, 2017;Tacconi, 2003), resulting in some of the highest rates of habitat loss worldwide (Gaveau et al., 2016;Miettinen, Shi, & Liew, 2011). Among the many species facing declines in this region, clouded leopards rank notably for their charisma , umbrella capacity , and ambassadorial potential . Species of conservation concern in their own right, clouded leopards are also powerful levers of conservation action for broader forest conservation programmes, including umbrella protection for diverse forest biota similarly threatened by widespread habitat loss (Collins, Milner-Gulland, Macdonald, & Macdonald, 2011).
The mainland clouded leopard, Neofelis nebulosa, ranges from the Nepali Himalayas in the west to southern China in the north and east, and extends south into Peninsular Malaysia; its sister species, the Sunda clouded leopard (N. diardi), occurs on Borneo and Sumatra (Buckley-Beason et al., 2006;Can et al., ;Kitchener, Beaumont, & Richardson, 2006). Both species are listed as Vulnerable on the IUCN Red List of Threatened Species (Grassman et al., 2016;Hearn et al., 2015), and India recently added N. nebulosa to its Recovery Programme for Critically Endangered Species (National Board for Wildlife, 2018). It is estimated that fewer than 10,000 adults remain of either N. nebulosa (Grassman et al., 2016) or N. diardi (Hearn et al., 2015). The two species diverged ~1.9 mya (Wilting et al., 2011), and while the degree of differentiation between the two species suggests they are as genetically distinct as lions and tigers (Buckley-Beason et al., 2006), similarities and differences between their ecological niches remain largely unknown. Understanding the habitat preferences of these reclusive and vulnerable felids is critical to inform conservation planning to mitigate further losses.
A recent assessment of habitat associations for N. diardi revealed that both recent forest loss and large-scale plantations strongly and negatively influenced detection rates (Macdonald, Bothwell, et al., 2018). Conversely, N. diardi were positively associated with forests, higher elevations and ridgelines. That study also found significant differences in poaching among regions, with substantially greater poaching pressure in Sumatra compared to Borneo; N. diardi detections decreased rapidly when even a few poachers were observed (Macdonald, Bothwell, et al., 2018). The recent comprehensive assessment of habitat use for the Sunda species now makes possible the comparison with N. nebulosa presented here.
To compare ecological niches between the two clouded leopard species, we developed a multi-scale habitat selection model that encompasses the full range of N. nebulosa from Nepal to Malaysia, following the same approach previously applied to N. diardi on Borneo and Sumatra (Macdonald, Bothwell, et al., 2018). It has long been recognized that species-habitat relationships are scale-dependent (Levin, 1992;Wiens, 1976Wiens, , 1989, and scale optimization has been shown to substantially increase model predictive power when compared with non-optimized, single-scale models (Timm, McGarigal, Cushman, & Ganey, 2016;Wan et al., 2017). Yet, a recent review found that <5% of current habitat modelling papers used robust approaches to optimize multivariate scale relationships (McGarigal, Wan, Zeller, Timm, & Cushman, 2016). Scalar relationships are a critical component of habitat use and provide valuable information for conservation management planning and reserve design; therefore, it is critical to incorporate multi-scale approaches into the statistical modelling framework when describing species-environment relationships (Thompson & McGarigal, 2002). Using N. nebulosa as an example, we demonstrate a two-step, multi-scale modelling framework, recommended as the most robust method currently available for multi-scale optimization . We first utilize a univariate approach to identify the optimal spatial scale for each environmental predictor, and second, combine scale-optimized predictor variables in a multivariate model to describe N. nebulosa's association with its environment.
Here, we aimed to identify N. nebulosa's primary habitat requirements, limiting factors and sources of threat. Specifically, we sought to (a) identify key environmental and anthropogenic variables influencing N. nebulosa habitat use, (b) determine the spatial scale at which each variable most strongly influences clouded leopard detection, and (c) draw comparisons between N. nebulosa and findings previously reported for N. diardi. The comparison between their habitat associations, which we make here for the first time, offers valuable insights into the ecology of these elusive species. Furthermore, these results provide critical information to assist in conservation management of N. nebulosa and the associated forest biodiversity for which it is an ambassador .  (Penjor et al., 2018), nine sites in Peninsular Malaysia were contributed by Tan et al. (2017), and E. Ash (unpublished data) contributed surveys from five locations in Thailand. Camera stations were primarily located in national parks, reserves, and other protected areas. In Bhutan, camera traps were widely deployed, irrespective of land protection status, resulting in 10/21 sampling locations in protected areas (Penjor et al., 2018). Given the challenge of detecting clouded leopards in dense forests (Wilting et al., 2011), cameras were primarily deployed along ridgelines, streams, and forest breaks (e.g., man-made trails, abandoned logging roads), where detection rates have been shown to be higher Macdonald, Bothwell, et al., 2018). Imperfect detection can bias estimates of true occupancy probability; however, our aim was rather to assess the overall strength and direction of relationships with environmental factors to provide guidance for N. nebulosa habitat conservation. For the purpose of assessing species-habitat relationships in rare or difficult-to-detect species, Banks-Leite et al. (2014) found that accounting for imperfect detection (e.g., with occupancy modelling) did not provide significantly more accurate results. Furthermore, unadjusted capture frequencies tend to be highly correlated with adjusted estimates (Kelly, 2008).

| Data collection
Thus, while we acknowledge the potential for variance in probability of detection among camera stations, relative strength and direction of relationships should be preserved. Paired camera stations were situated ~40 cm above ground and spaced 1-2 km apart. This camera density was originally designed for estimating spatial capture-recapture density. Data were spatially rarefied to one station/1.0 km 2 to negate pseudo-replication for subsequent modelling, resulting in removal of 56 camera stations. We further tested for and found very low levels of spatial autocorrelation (Moran's ̄I = 0.02; Table 1; 'spdep' r package; Bivand & Piras, 2015), thereby validating the assumption of independence among samples.
TA B L E 1 Collection information, including nine countries, 45 sampling locations, camera traps/site (n), total N. nebulosa detections, detection rate, mean number of trap nights (camera effort), and Moran's I spatial autocorrelation and significance

| Covariates
From previous regional studies, we identified a suite of predictor variables related to N. nebulosa habitat use. Closed-canopy forests are a primary requirement for this semi-arboreal species Sollmann, Linkie, Haidir, & Macdonald, 2014;Tan et al., 2017). Conversely, deforestation and subsequent land conversion to large-scale palm and Acacia plantations have been identified as major threats for both the mainland (Tacconi, 2003) and Sunda Hearn et al., 2017;Hearn et al., 2019;Macdonald, Bothwell, et al., 2018; species. Previous studies of N. nebulosa in Bhutan (Penjor et al., 2018) and of N. diardi in Borneo and Sumatra (Hearn et al., 2016;Macdonald, Bothwell, et al., 2018;Sollmann et al., 2014) also found positive associations with ridgelines and slope, and negative associations increased with density of human settlements and land use intensity.
Based on careful biological consideration of these relationships, we identified 13 predictor variables that we hypothesized are driving Neofelis habitat use. We then transformed these into 46 more biologically informative variables using class-and landscape-level spatial statistics (Table S1). To facilitate comparison with the recently published N. diardi habitat model (Macdonald, Bothwell, et al., 2018), we used the same covariates here, plus mean annual precipitation and temperature (MAP and MAT) to account for substantial climatic variability across the large study region. by similarity or contrast among habitat types. We generated hypothesized edge density weightings such that moving between very similar habitat types incurred a minimum weighting of 0, whereas maximally contrasting habitats incurred a weighting of 1 (Table S2).
For example, we hypothesized that N. nebulosa would perceive the difference between closed forest and urban areas as high contrast (CWED = 1), whereas moving from shrubland/grassland to mosaic cropland would be a less abrupt transition (CWED = 0.5). Because CWED is a metric related to habitat fragmentation, we hypothesized that clouded leopard detection would exhibit a negative association with higher CWED values.
To assess scalar relationships between N. nebulosa and its en-

| Data analysis
In Step 1 of the multi-scale modelling approach, we identified the optimal scale and functional form (linear or quadratic) for all speciesenvironment relationships via univariate generalized linear mixedeffects models (GLMMs; 'lme4' package (Bates, Maechler, Bolker, & Walker, 2015)) in r 3.3.2 (R Core Team, 2016). We included camera effort (i.e., total active nights per camera trap) as a fixed effect and sampling location nested within country as random effects. We selected best-supported scales and functional forms for each variable based on Akaike's information criterion, adjusted for small sample size (AICc; Burnham & Anderson, 2002). We note that the 13 origi- Thus, testing a limited number of a priori 'best-guess' hypotheses at best provides limited scope and insight into complex relationships, and at worst is easily corrupted by investigator bias. Data-driven approaches that identify best-supported relationships by more fully interrogating the true complexity of the underlying hypothesis space can provide a more objective understanding of species-habitat relationships.
We then applied four filtering steps to reduce the number of variables included in the multivariate model (Table S3). (a) Twelve variables were removed that occurred at <10% of camera stations.
Lack of representation of these habitat types at sampling locations resulted in insufficient data to reliably test their influence on clouded leopard detections. (b) AIC model selection identifies the best model relative to a chosen set of models, but it does not provide information on strength of relationship or effect size.
Therefore, we also evaluated model performance of the top univariate GLMM selected for each variable and further required that optimal models identified via AIC also demonstrated strong evidence of relationship with the response variable (Macdonald, Bothwell, et al., 2018;Šímová et al., 2011). It is possible for uni- contrast to other countries also impacted by the war (e.g., Cambodia, Laos) but lacking international market access (Dudley, Ginsberg, Plumptre, Hart, & Campos, 2002
Additionally, mean annual precipitation (MAP), which influences vegetation density and biomass, was selected at the 32-km radius scale. Protected area correlation length (i.e., extent) and mean slope position were best-supported at the 8-km radius scale. Fine-scale relationships included mean compound topographic index (CTI; 500-m radius) and percentage mosaic habitat (1-km radius). Furthermore, quadratic relationships were best-supported for all variables except mean MAP and shrubland/grassland extent.

| Multi-scale model selection and validation
AIC model selection produced four top models with ΔAICc ≤2 (Table   S4); model averaging then identified a final model with six variables (Table 2). Percentage closed forest was the strongest predictor of N. nebulosa detections (16-km radius focal landscape). When closed forest habitat increased from 65 to 100% of the landscape, we observed ~25% increase in detections (Figure 1a). MAP was also strongly, positively correlated with detections (32-km radius).
Regions receiving <170 cm MAP were associated with low detection frequencies (<0.12); however, as precipitation increased above 170 cm, we observed a steady increase in detections (Figure 1b).
Broad extent of protected areas was also positively associated with increasing detections (8-km radius). Variables negatively associated with N. nebulosa detection frequency included percentage mosaic habitat (1-km radius) and mean CTI (500-m radius; Figure 1d).
Using MESS analysis, we identified geographic locations with environments falling outside the range occurring at sampling locations used for model training, thus representing extrapolation ( Figure   S2a). Due to the cost of setting up and maintaining camera grids, and the challenge of detecting the rare and notoriously illusive N. nebulosa, limited resources necessitated placing cameras where there was some chance of detecting the species. While sampling locations were biased towards protected areas, locations were specifically chosen to maximize variation encountered by N. nebulosa across its range throughout Southeast Asia. Overall, multivariate environmental space at training data locations was highly similar to conditions encountered throughout the core of N. nebulosa's range ( Figure S2a).
Non-analog environments were identified along the species' range margins, where training data were more limited. Univariate similarity assessments revealed that lack of closed forest and high CTI were the strongest drivers of dissimilarity in central and south-western China, India, Bangladesh, and central Thailand. Very high precipitation along the western coastal lowlands of Myanmar and lower precipitation in central and south-western China also contribute to non-analog environments ( Figure S2b). Model inference should be considered cautiously in these regions.

| Regional variation
The projected model revealed that 9.44% of the Southeast Asian landscape is highly suitable for N. nebulosa (Figure 2; representing top 10% of pixels), yet this varies substantially among countries ( Figure S1). Neofelis nebulosa's core distribution coincides with Laos, Malaysia, and Myanmar, with 52.5%, 35.4%, and 30.7% of their landscapes harbouring highly suitable habitat, respectively. These countries also boast the largest mean patch sizes, greatest contiguous patch extensiveness and highest patch densities-highlighting their importance for clouded leopard conservation (  Figure S5).
We summarized multivariate habitat configuration and extent relationships among nations via PCA ( Figure S3a-c). The first three axes collectively explain 88.5% of total variance among metrics (Table S6). PC1 is strongly aligned with percentage and extensiveness of high-quality habitat (≥90th percentile). Laos and Myanmar, followed by Malaysia, Bhutan and India, are characterized by high values ( Figure S3a); conversely, Bangladesh, China, Nepal, Cambodia and Thailand exhibit low percentage and extensiveness of high-quality habitat. Along PC2, China and Myanmar have the greatest extensiveness of medium-to high-quality habitat (≥50th percentile).
Interestingly, China is an outlier, with a low-to high-quality habitat extent, but given its vast size, the highest area of medium-to highquality habitat.
The cluster dendrogram shows strong partitioning into five hierarchical groups (Figure 3) the best ability to discriminate among clusters (Table S7), suggesting this metric can be a simple way to compare habitat patterns among countries. A bar plot of correlation length shows clear differentiation among countries and is highly consistent with the dendrogram clusters identified above (Figure 4). Additionally, we quantified the summed pixel values of predicted habitat suitability for each nation to provide a simple, intuitive and non-threshold-based comparison.
This metric similarly ranks Myanmar and Laos well above the other nations in total predicted habitat suitability ( Figure S4).
To assist with prioritization of habitat patches for conservation management, we identified the top 28 patches >1,000 km 2 , ranked by descending high-quality habitat area ( Figure 5; Table S8). The largest patch (1) extends throughout Laos and western Vietnam, and is 1.8 × larger than the next largest patch. We also identified Myanmar as a critical hub of N. nebulosa habitat; harbouring four of the top six patches (2, 3, 4, 6), these together cover 1.5 × the area of Patch 1. The 5th largest patch also identifies important habitat in Malaysia for future conservation focus. Assessing overlap between the top patches and protected areas ( Figure 6,  Our projected model provides crucial information to assist conservation management, including the identification of highly suitable core habitat and medium-quality habitat likely critical to clouded leopard meta-population viability through its provisioning of essential connectivity corridors for dispersal and mating among core populations. This is also the first empirically based comparative assessment of environmental niche space of two keystone, Southeast Asian carnivores.

| Niche comparisons
Using the same modelling approach as Macdonald, Bothwell, et al. (2018) allowed us to compare habitat models between the mainland clouded leopard and its allopatrically distributed sister species, the Sunda clouded leopard. We acknowledge that lack of data for prey availability and intra-guild competitive dynamics limit the utility of habitat models to accurately describe these species' realized niches.
However, given that clouded leopards are threatened by rapid habitat loss, these models provide a timely assessment and represent the best models currently available for these species. Cross-model comparison identified many of the same variables and spatial scales as being important for both species, supporting model reproducibility and suggesting substantial niche conservatism at the genus level.
F I G U R E 2 Predicted N. nebulosa habitat suitability across Southeast Asia, showing 50th, 75th, 90th and 97.5th percentiles of habitat suitability TA B L E 3 Comparative analysis of landscape metrics for high-quality (≥90th percentile) and medium-to high-quality (≥50th percentile) thresholds of the projected model. ≥97.5th and 75th percentile analyses are presented in Closed-canopy forest was the strongest predictor of clouded leopard detection; few N. nebulosa were observed when the landscape was <65% closed forest. Furthermore, detections were most strongly correlated with forested landscapes at the 16-km radius scale, consistent with previous estimates of N. diardi's relatively broad-scale movements (Hearn et al., 2013;. Both species are strongly associated with landscapes dominated by extensive forest cover at broad spatial scales. Both the mainland and Sunda models also showed strong negative associations with mosaic cropland habitat and large-scale plantations, respectively. Deforestation and agricultural intensification are clearly dire threats to clouded leopards. Neofelis nebulosa detections were also strongly, positively associated with increasing precipitation. Our model suggests a threshold response; low detections were observed below 170 cm MAP and steadily increased above 170 cm. This association is particularly informative for predicting future climate change risks to clouded leopards. Regions that experience increasing aridity will have reduced capacity to support high-density forests, and by association, clouded leopards. Although percentage availability of highly suitable habitat is comparable among mainland Southeast Asia (9.44%), Borneo (10.04%) and Sumatra (8.98%) (Macdonald, Bothwell, et al., 2018) Table 2). While the difference in optimal scale may reflect less extensive habitat available to N. diardi, it may also reflect differences in prey availability, which is known to strongly influence habitat selection in clouded leopards (Mohamad et al., 2015).

Comparative analysis revealed another distinction between
the Neofelis species with respect to shrubland/grassland habitat.  (Bothwell et al., 2017;Cushman et al., 2018;Eckert, Samis, & Lougheed, 2008). Divergence in optimal CTI may also reflect greater deforestation intensity in Borneo and Sumatra. Vast

| Future research needs
Clouded leopards co-occur with a varying guild of felids and prey species across their range. The habitat selection model developed here provides a timely assessment to support critical habitat conservation initiatives for a species facing rapid habitat loss; however, both intraguild competitive dynamics and hunting opportunities are influential drivers of habitat selection for these carnivores (Mohamad et al., 2015). Detectability may be affected by larger felid trail use, as observed in other mammal communities (Harmsen, Foster, Silver, Ostro, & Doncaster, 2010). The distribution of and variation in prey species, including a wide variety of birds, squirrels, monkeys, deer and wild pigs, may contribute to a narrower realized niche than observed for the current habitat model. For example, we suspect the positive association of N. nebulosa with shrublands and grasslands may be driven by lower poaching intensity and greater prey opportunities in this habitat type on the mainland in contrast to N. diardi in the Sunda Islands. We are currently exploring both intra-guild competitive dynamics and associated community biodiversity patterns, which we expect will provide additional insights into clouded leopard habitat selection.
Additionally, poaching poses a critical threat to both species.

| Scale and conservation management planning
Multi-scale modelling provides an efficient means to optimize conservation management planning to the scales at which environmental variables impact species occurrence. For clouded leopards, this suggests the need for broad-scale, landscape-level reserve network design, with targets for forest cover >65% and a minimum 800 km 2 of highly suitable habitat. Our scale optimization suggests a minimum 8 km buffer around protected areas should be managed for low coverage of plantations and mosaic cropland (Table 2). To assist in prioritizing continuing conservation efforts, we identified the top 28 high-quality habitat patches ( Figure 6; Table S7) and highlight key gaps in protection and connectivity ( Figure 6).
Three key recommendations emerged based on overlap be- broad spatial scales (>800 km 2 ), highlight clouded leopard's value as an ambassador for conservation of broader forest biodiversity , and support its capacity to act as an umbrella for the protection of co-occurring Southeast Asian species similarly threatened by widespread deforestation (Macdonald, Burnham, Hinks, & Wrangham, 2012;Macdonald et al., 2017;Roberge & Angelstam, 2004).

ACK N OWLED G EM ENTS
The majority of the team and data were part of the core WildCRU effort supported by a Robertson Foundation Grant to DWM, for which we are deeply grateful. WildCRU personnel were assisted by partner organizations and staff at every sampling location; we warmly acknowledge their assiduous and tireless work. We ac- National Protected Area crew members and generous in-kind support. We also thank WWF Cambodia, FFI, and the Ministry of Environment for supporting the surveys and providing personnel.

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.

DATA ACCE SS I B I LIT Y
All input GIS layers, the habitat suitability model and the high-quality habitat patch map are archived on Dryad. Given the extremely sensitive nature of clouded leopard occurrence data with respect to illegal wildlife trade, locations of camera traps will not be made public so as to avoid further endangering the species. However, we welcome correspondence with scholars and conservationists regarding collaborative use of the data to advance science and conservation of clouded leopards, their associated communities and habitats.