Prioritizing areas for conservation outside the existing protected area network in Bhutan: the use of multi-species, multi-scale habitat suitability models

Understanding the environmental and anthropogenic factors influencing habitat selection of multiple species is a foundation for quantifying human impacts on biodiversity and developing effective conservation measures. To determine the effect of multiple scales of environmental/topographic and anthropogenic variables and landscape patterns on habitat suitability of terrestrial mammals in Bhutan, assess the effectiveness of the current protected area network, identify areas of high species richness outside of the existing protected area, and evaluate the potential effectiveness of indicator and umbrella species for conservation planning. We modelled multi-scale habitat selection of sixteen species of terrestrial mammals across Bhutan using data from a nation-wide camera trap survey. We used the predicted species distribution maps to assess the multi-species conservation effectiveness of the existing protected area network. We performed simulations to identify high priority areas for multiple species based on their habitat suitability, proximity to existing protected areas and overall connectivity within the predicted distribution of species. We used correlation analysis among predicted occurrence maps and multivariate cluster analysis to identify potential indicator species. We evaluated the potential utility of each species as umbrella species by assessing how well optimal protected areas for that species would protect suitable habitat for all 16 species simultaneously. Protected areas and forest cover were strongly associated with habitat use of most modelled species. Additionally, topographical features, like terrain roughness and slope position, contributed to habitat selection of multiple species, but often in different ways. Environmental and topographical variables were mostly selected at medium to broad scales. Anthropogenic variables (agriculture and built-up areas) were negatively associated with habitat suitability of most species at both fine and broad scales. Conservation effectiveness assessment of existing protected areas found protected areas in south-central Bhutan have high effectiveness in terms of both mean and total richness protected. Similarly, biological corridors in the south-central region offered high mean richness protection. Our simulation of optimal areas for additional protection found areas abutting protected areas in southern Bhutan offered high relative species richness protection. Our umbrella species analysis found muntjac, wild pig, serow, sambar and Asian golden cat are the most effective umbrella species for broader biodiversity protection. Our indicator species analysis found tiger, gaur, dhole, clouded leopard, Asian black bear and common leopard as effective indicator species. This study highlights the need to protect optimally located species-rich areas outside the current protected areas. This kind of multi-species habitat assessment provides important information to optimize future conservation and development plans at national and regional scales.


Objectives
To determine the effect of multiple scales of environmental/topographic and anthropogenic variables and landscape patterns on habitat suitability of terrestrial mammals in Bhutan, assess the effectiveness of the current protected area network, identify areas of high species richness outside of the existing protected area, and evaluate the potential effectiveness of indicator and umbrella species for conservation planning. Methods We modelled multi-scale habitat selection of sixteen species of terrestrial mammals across Bhutan using data from a nation-wide camera trap survey. We used the predicted species distribution maps to assess the multi-species conservation effectiveness of the existing protected area network. We performed simulations to identify high priority areas for multiple species based on their habitat suitability, proximity to existing protected areas and overall connectivity within the predicted distribution of species. We used correlation analysis among predicted occurrence maps and multivariate cluster analysis to identify potential indicator species. We evaluated the potential utility of each species as umbrella species by assessing how well optimal protected areas for that species would protect suitable habitat for all 16 species simultaneously. Results Protected areas and forest cover were strongly associated with habitat use of most modelled species. Additionally, topographical features, like terrain roughness and slope position, contributed to habitat selection of multiple species, but often in

Abstract
Context Understanding the environmental and anthropogenic factors influencing habitat selection of multiple species is a foundation for quantifying human impacts on biodiversity and developing effective conservation measures.

Introduction
Humans have impacted biodiversity for millennia and are linked to the decline of megafauna and prehistoric extinctions (Bartlett et al. 2016). However, the recent explosion in human populations, coupled with increasing per capita resource consumption, has accelerated these impacts, with more than 50% of global land cover now altered by human activities (McGill et al. 2015). These rapid land-use changes have driven large declines in wildlife populations worldwide (WWF 2016). As a result, it is estimated that current extinction rates are as much as 1000 times the background rate in the fossil record (Johnson et al. 2017).
Global initiatives to address the extinction crisis and improve biodiversity conservation have produced ambitious visions. Amongst them, the Aichi Biodiversity Targets 2020 aims not only to reduce species loss but also to safeguard habitat and acknowledge ecosystem services, while also improving the management and sustainable use of natural resources (Tittensor et al. 2014). Despite such efforts, the rate of global biodiversity loss has likely accelerated substantially in recent decades (Tittensor et al. 2014). Therefore, more effective conservation programs that optimize protection of critical habitat are essential to stem biodiversity loss and enhance species recovery (Venter et al. 2014).
Animals have varied requirements for habitat extents, patterns and characteristics (Torres-Contreras et al. 1997;Kelt et al. 1999;Finlayson et al. 2008) and their habitat use and distribution are influenced by biotic and abiotic factors, such as food availability, predation, guild interaction and competition (Falkenberg and Clarke 1998). Thus, knowledge of factors influencing species habitat use is a prerequisite to habitat protection, conservation and management. Further, species respond to environmental, topographical and anthropogenic features at different scales (Wiens 1989;Collingham et al. 2000;McGarigal et al. 2016). Specifically, scale plays an important role in resource and space use of organisms in heterogeneous landscapes (Kotliar and Wiens 1990;Holling 1992). Different species select habitat at different scales (Roland and Taylor 1997), and a given species may select different resources at different scales (Grand et al. 2004;Wasserman et al. 2012). Therefore, it is critical to incorporate the effects of scale in species-habitat relationships and more importantly in planning and designing reserves (Thompson and McGarigal 2002).
To stem species loss and habitat conversion, global protected area (PA) designations have increased in extent and number in the last few decades . However, PAs have often been designated without sufficient strategic planning, and as a result, existing PAs are often inadequate or inefficient in their protection of critical biodiversity and their habitats (Jenkins et al. 2015;Chundawat et al. 2016;Johansson et al. 2016). Assessment of the biodiversity effectiveness of existing PAs is often lacking (Chape et al. 2005). It is largely unknown whether PAs adequately represent biodiversity and ensure sustainable protection (Margules et al. 2002). Pragmatically, it is impossible to protect all wilderness areas; hence, it is important to identify critical biodiversity areas and prioritize them for intensive conservation (Waldron et al. 2013).
When prioritizing land use and protected area establishment, past studies have frequently intersected connectivity models with protected areas and prioritized areas for protection based on the strength, centrality or connectedness of areas outside the existing connectivity network Kaszta et al. 2020;Ahmadi et al. 2020;Ashrafzadeh et al. 2020). These approaches, however, don't use formal optimization to select or prioritize new areas for protection and often have been limited to one ) or a few (Cushman et al. 2013) species. Formal optimization approaches, such as MARXAN (Ball et al. 2009), have the advantage of being objective, stochastic and replicable, which can improve the efficacy and robustness of predictions. However, past optimization approaches have been limited for use in connectivity analyses as they typically are based on habitat suitability and do not formally include connectivity or spread algorithms within them. In this study, therefore, we wrote a new simulation model to select areas of the highest priority for multiple species based on their habitat quality, the proximity of their habitat to existing protected areas, and its overall connectivity within the full predicted distribution of the species across Bhutan.
The Eastern Himalaya is a global biodiversity hotspot, with one of the world's highest rates of species endemism but is under pressure from a high rate of habitat loss driven by anthropogenic disturbance (Myers et al. 2000;Dorji et al. 2018). Bhutan, within the Eastern Himalayan biodiversity hotspot (Myers et al. 2000), currently has over 70% of land area under forest cover and more than 50% of the total land area receives official protection, supporting a broad spectrum of species and ecosystems (DoFPS 2011). However, as a developing country, Bhutan is beset with challenges to biodiversity conservation, including road development, urban expansion, human-wildlife conflict, wildlife poaching and depletion of natural resources. Hence, validating the effectiveness of PAs in the face of human land-use change and climate warming is imperative. Additionally, Bhutan has a constitutional mandate to protect 60% of its total land as forest cover. Although the current protected area and biological corridor network in Bhutan (hereafter PAN) is extensive (constituting 51% of the national area, 19,581 km 2 ), the efficacy of the protected area network in protecting wildlife, particularly megafauna, is largely unknown. Further, rapid human infrastructure development calls for identification and urgent protection of key biodiversity areas outside the PAN for the sustainability of conservation efforts. Lastly, there is a general interest in the conservation and scientific community regarding the use of indicator and umbrella species in conservation planning (Dufrêne and Legendre 1997;Roberge and Angelstam 2004). In this paper, we define indicator species as those whose presence indicates (by high correlation with) the presence or absence of other species, and umbrella species as those whose protection would simultaneously protect optimal habitat for many other species.
To address these knowledge gaps and management needs, this study aims to: (1) investigate the multiple scale habitat relationships of 16 terrestrial mammal species of conservation importance and predict their distribution across Bhutan, (2) evaluate the effectiveness of the existing PAN based on the amount of total species richness protected, (3) identify areas of conservation importance outside of the existing PAN that could be chosen to optimally expand the protected area network, and (4) evaluate the efficacy of each of the 16 focal species as conservation umbrella and indicator species. Our models are replicable, enabling incorporation into long-term monitoring of critical biodiversity areas by governments and conservation agencies.

Study area
This study was conducted in Bhutan as a part of the nation-wide tiger survey in 2014-15. Bhutan, with an area of 38,394 km 2 , is a small country in the eastern Himalayas between India and China (Fig. 1). About 70% of the total geographical area is covered by forest and another ~ 10% is a natural shrub and alpine scrub cover and ~ 4% under agriculture and built-up area (FRMD 2017). Vegetation varies with elevation; lower foothills are characterized by broadleaved forest, the temperate zone by conifer-broadleaved mixed forest and higher elevations by alpine and sub-alpine 1 3 scrub and shrub. The altitude ranges from 100 m in the south to above 7500 m in the north. The wet season (monsoon) occurs between July and September, with annual precipitation ranging between 300 and 6000 mm (NCHM 2017). The topography is rugged, with steep terrain, narrow river valleys and deep gorges. Bhutan has ~ 120 species of mammals, ~ 700 birds and ~ 5000 plants (Gyeltshen et al. 2020).

Field survey
The nation-wide camera-trap survey was conducted between March 2014 and March 2015 with a break during the intervening rainy season. For logistical reasons, the country was divided into two blocks: north (n = 681) and south (n = 448), (Table S1). Five different camera models were deployed (Roconyx, HCO Scoutguard, Cuddeback, U-Way, Bushnell). All stations had paired cameras, and none was baited. The deployment followed a 25km 2 grid based on the size of a female tiger home range (Karanth et al. 2002) and excluded settlements and areas above 4500 m (DoFPS 2015). The minimum distance between stations was 2 km. The monitoring team checked camera stations for battery drain, cleared bushes obstructing the camera lens and downloaded memory cards every month.

Predictor variables
Resource use by a species is influenced by natural and anthropogenic factors. While topographical factors impose physical constraints and drive major patterns of climate and vegetation, human land use often perturbs ecological systems, altering their structure, composition and suitability for many species (Magioli et al. 2019;Macdonald et al. 2020). We selected 20 predictor variables reflecting environmental, topographical and anthropogenic conditions that we hypothesized would affect habitat suitability of the 16 studied species (Table 1). These variables, as described below, were further used to calculate six landscape metrics, which quantify the composition and configuration of environmental patterns, (McGarigal et al. 2012) (Table 1).
To predict effects of topography on habitat selection, the Shuttle Radar Topography Mission (STRM) digital elevation model (Jarvis et al. 2008) was transformed into topographical roughness (which measures the complexity of the terrain), slope position (measuring the relative position of a location from valley bottom to ridge top) and compound topographic index (CTI) using the Geomorphometry and gradient metrics toolbox (Evans et al. 2014) in ArcGIS 10.6.1 (ESRI 2018). CTI measures the flow accumulation (wetness index) such that broad valleys with large drainage areas have high CTI and steep mountain ridges have low CTI values (Beven and Kirkby 1979;Macdonald et al. 2018). Forest cover (Hansen et al., 2013) was classified into three classes: non-forest (0-20% tree cover), open-forest (21-40% tree cover) and closed forest (> 41% tree cover) to enable an analysis of the spatial pattern of forest density and canopy closure on species occurrence patterns. Land cover (FRMD 2017) was classified into eight categories: agriculture, bare, broadleaved forest, built-up, conifer forest, meadow-shrub, water bodies and snow-glacier areas. We also included the spatial layers of settlement, road, river densities and the extent of protected areas (Table 1).
To identify spatial extents at which habitat components and anthropogenic factors most strongly affect species presence, we calculated each predictor variable for each species at seven spatial scales: 500 m, 1 km, 2 km, 4 km, 8 km, 16 km and 32 km. This scale range was chosen to bracket the range of expected optimal scales for each species' response to each variable (e.g. Mateo Sánchez et al. 2014;McGarigal et al. 2016;Macdonald et al. 2018Macdonald et al. , 2019 and is consistent with the range of scales in other multi-scale studies of wide-ranging species in Asia (Macdonald et al. , 2019. We categorize these scales as fine-scale (500 m, 1 km, 2 km), medium-scale (4 km, 8 km) and broad-scale (16 km, 32 km).
For settlements, rivers, roads and protected areas we ran a focal mean (FM) analysis in ArcGIS 10.6.1 (ESRI 2018) by applying a circular moving window of a radius equal to each of the seven investigated scales. To investigate the effects of spatial configuration and composition of the landscape mosaic at different scales on habitat use, we calculated landscape metrics on the reclassified land cover layer described above and the three classes of forest cover (McGarigal et al. 2012). These metrics included: PLANDpercentage of the landscape of a particular habitat class; CWED-contrast-weighted edge density which is a measure of edge effect of similar or contrasting habitats; GYRATE_AM-correlation length which measures the expected value of the distance from a random location in a cover type to the edge of that cover type moving in a random direction and is an effective measure of patch extensiveness; and PD -patch density. These metrics were chosen based on past research showing them to be among the metrics that most strongly and consistently affect wildlife habitat selection (e.g. Grand et al. 2004;Chambers et al. 2016). The metrics were calculated for all forest cover and land cover classes within a circular moving window of a radius corresponding to the chosen seven scales (Table 1)  All initial raster maps were resampled to 90 m resolution. The value for each predictor variable at each focal scale was then extracted at each camera station. All 308 predictor variables were standardized to mean 0 and standard deviation 1 before analysis. After removing variables which lacked sufficient information due to small coverage (e.g., unvarying in > 95% of stations), the final number of variables used in the multi-scale modelling was 291.

Multi-scale habitat suitability models
Multi-scale habitat suitability models (sensu McGarigal et al. 2016) were used to identify key environmental and anthropogenic variables that were predictive of the occurrence of each species, and the optimal scale of species-habitat relationship for each variable (Macdonald et al. , 2019. The species detection data at each camera site were collated into binary form, with 1 for detection and 0 for non-detection. We then employed a three-step modelling approach to optimize scale and assess species-habitat relationships, as suggested by (McGarigal et al. 2016). All analyses were performed in R (R Core Team 2019).
In the first step, we fit univariate generalized linear models (GLM) with a binomial link function for each species as a response variable and each variable at each scale as the predictor variable. The best scale for each variable for each species was selected based on the lowest AICc score (Akaike Information Criterion corrected for small sample size) (Burnham and Anderson 2002).
In the second step of the modeling process, we screened predictor variables for multi-collinearity (Dormann et al. 2013). We used Pearson's correlation coefficient threshold of |r|≥ 0.6 to screen out correlated predictors. We ran univariate GLM for each correlated pair and retained the one with lower AIC value. We also assessed spatial autocorrelation using Moran's I statistic in the spdep R package (Bivand and Wong 2018). When significant spatial autocorrelation was found, we added an auto-covariate term with a distance threshold of the mean distance between camera traps in the final multivariate model.
In the third and final step, we ran multivariate GLM with a binomial link on all remaining scaleoptimized predictor variables. We then used the 'dredge' function in R package MuMIn (Barton 2019) to test all possible combinations of predictor variables and ranked models using ∆AICc (the difference between AICc of subsequent models with the top model). The final models were considered competitive if they were within ∆AICc of 2. We employed model averaging based on AIC model weights to produce final coefficient estimates (Burnham and Anderson 2002), which we then used to map habitat suitability for each species.
We evaluated the performance of the model-averaged predictions using the area under the receiver operating characteristic curve (AUC). The mean habitat suitability of each study species across Bhutan was calculated as the mean of all pixels with suitability value > 0 using the Zonal Statistics tool in ArcGIS 10.6.1 (ESRI 2018). Twelve variables had three metrics (PLAND/CWED, GYRATE_AM and PD) and each metric had seven scales resulting in 12 × 21 = 252 predictor variables, and eight variables had one metric (FM) at seven spatial scales resulting in 8 × 7 = 56 variables. The total variables initially tested were 308

Conservation effectiveness of existing PAN
To evaluate conservation effectiveness, we calculated the mean predicted relative species richness within each protected area and biological corridor and the proportion of the total sum of richness across Bhutan (per cent richness) in each protected area and biological corridor. To calculate species richness, we summed the predicted habitat suitability layers of the 16 species. We then calculated the mean richness value for each protected area and biological corridor and the percentage of the sum of species richness represented by each protected area and corridor. The cumulative layer of the sum of predicted occurrence probability for all 16 species was used to identify hotspots.

Simulation to identify additional conservation areas
We wrote a new model in Arc Macro Language for ArcInfo Workstation Grid to simulate optimal locations for protected area expansion, given its high efficiency for spatial analysis and functionality of nested and looped spatial functions. We designed the model to optimize the selection of areas outside of the current protected area network that had simultaneously high local habitat quality and high connectivity to existing protected areas for multiple focal species. We ran the model for each of the 16 focal species and the sum of predicted probability of all focal species (17 simulation scenarios). The routine consists of five components and the simulation works through two nested iterative loops (Fig. 2). The basic method is to select additional cells for new protected areas for each species, or all species jointly (depending on scenario m), based on two criteria: (1) habitat suitability for that species or collective total sum suitability across species, and (2) cost distance from existing protected areas across the suitability surface. This results in new cells being stochastically selected based on the combination of their quality as habitat and their cost proximity to cells previously selected, which ensures that new protected areas are both high-quality habitat and well connected to current and additionally selected protected area cells.
There are several steps in the iterative simulation process. First, we create a difference raster by subtracting the inverted and rescaled cost distance (between 0 and 1) from protected area cells across the resistance surface (Fig. 2d) and a uniform random raster (Fig. 2e). This creates a stochastic process to select cells of low relative cost distance. This is done by selecting all cells with value > 0.75 from the difference raster. A resistance surface, as used in cost distance modelling, reflects the cost of resistance to movement of each cell for the process being modelled (Spear et al. 2015). In this case, we are modelling connectivity for a focal species or for a community of focal species (as the scenario may be), and the resistance layer was defined as the inverse of the predicted occurrence probability layer (or the sum of predicted occurrence probability layers, for the multispecies scenario), rescaled between 1 and 100. In each time step of the iterative model, the algorithm described in the first step above selects all cells that are 75% greater than the random value generated for that cell and results in a small proportion (~ 0.25% per time step) of the most optimal cells being selected. The choice of a value of 0.75 for this threshold is subjective but was calibrated such that it produced in each time step a selection of cells that was sufficiently large to seed a meaningful optimization but small enough to avoid seeding many suboptimal pixels.
In the next iterative step, the pixels that were stochastically selected from the most well-connected cells in the landscape are then added to the raster of protected areas (Fig. 2c) and the iterative process just described is repeated 50 times, creating additional optimal stochastically selected protected area cells. Once this loop is completed the entire process is repeated in another exterior loop (Fig. 2) 50 times. The choice of 50 iterations for both loops was intended to provide enough Monte Carlo runs to produce a precise estimate while minimizing processing time.
At the end of the total simulation, the protected area rasters produced are summed, creating values between 0 (never selected as an optimal new protected area) and 50 (selected as an optimal protected area in each exterior iteration). This provides a quantitative measure of how optimal each cell is for inclusion in a new protected area based jointly on habitat quality and connectivity for each focal species.
We computed these optimal protected area selections for all 16 species individually, and for all 16 species jointly. The prioritization for all 16 species jointly was done in two ways: (1) based on the sum of their predicted probability maps, (2) based on a weighted sum of predicted probability maps. The weighting reflects the relative conservation priority of each of the 16 species based on their local and global conservation status (Table S2). Specifically, in the first case, we summed the predicted probability of occurrence maps for all 16 species without any transformation, while in the second we weighted the sum such that the species with higher conservation importance (Table S2) received more priority. As noted above, these predicted probability surfaces (for each species scenario) and the sum of predicted probability surfaces (for the two multi-species scenarios), were transformed through inversion and rescaling between 1 and 100 to produce resistance layers for use in the spatial simulation.
For each of these 18 scenarios, we produced the raster of prioritization (0-50) for cells to include in new protected areas. We present all of these rasters in the supplemental material ( Figure S22). We include those for several focal species of particular interest and one of the multispecies optimizations as figures in the Results. We also evaluated the relative similarity and overlap of optimal cells for expansion of the protected area network between scenarios (16 species and 2 multispecies) by computing the pixel-wise correlation between values (0-50) in the output prioritization rasters.

Fig. 2
Schematic of the simulation model developed to stochastically optimize new protected areas based on habitat quality, proximity to existing protected areas and connectivity across the population. There are six inputs/outputs (a-f) and five iterative steps (1-5). The inputs/outputs are: a habitat suitability map, b resistance map produced from (a), c existing protected areas, d cost distance from existing protected areas across the resistance map (b), e uniform random raster for stochastic simulation, f selection of new optimal cells for addi-tional protected areas as the difference of (d) and (e). The iterative steps are (1) conversion of habitat suitability to resistance (outside of all iterative loops), (2) calculation of cost distance from protected areas on the resistance surface (inside iterative loops), (3) stochastic selection of new protected areas (inside iterative loops), (4) iterative loop 1, 50 repeats of the sequential additive selection of new protected area cells, (5) iterative loop 2, repetition of the full selection 50 times

Umbrella and indicator species
We assessed the effectiveness of each species as umbrella species, in which conservation of optimal habitat for that species would simultaneously conserve habitat for the others. The umbrella species analysis was done by (1) computing the average and (2) sum of habitat suitability for each species within the optimal protected area cell selection raster (from the simulation analysis described in the previous section) for all species sequentially, which quantifies how well the optimal protected area network for each species would protect habitat for the other species. Furthermore, we evaluated the efficacy of the two multi-species prioritization strategies in providing optimal protection for each species, by calculating the mean and sum of habitat suitability for each species within the optimal protected area polygons developed for each of the multispecies selection strategies. We also evaluated the relative ability of each species to serve as an indicator species for the sampled mammalian community using cluster analysis and correlation analysis. Cluster analysis generated groups of species based on the similarity in predicted occurrence probabilities. The pairwise correlations calculated between the predicted occurrence probability rasters produced from the multi-scale habitat modelling among all pairs of species identified which species had the highest average correlation with all other species (Cushman et al. 2010).

Results
We retrieved images from 849 camera stations (of 1129 initially installed; Table S1). Fifty-four mammal species were recorded over 73,259 trap days. The overall mean trap nights per camera was 86.7 (21.3). The number of species detected per camera station ranged from 0 to 17. The naïve occupancy (stations with species detection divided by the total number of camera stations) was highly variable among the study species with highest observed for muntjac Muntiacus muntjak (0.58) and lowest for Asian elephant Elephas maximus (0.10).

Habitat suitability models
About 39% of the predictor variables were selected at broad scales while small and medium scales were selected in equal proportions (~ 31% and 30% respectively) across all species (Table S4). Multi-scale single species habitat suitability models showed that forest-related variables were strong predictors of Asiatic golden cat Catopuma timminckii, muntjac, wild pig Sus scrofa, yellow-throated marten Martes flavigula, clouded leopard Neofelis nebulosa, serow Capricornis thar, tiger Panthera tigris and dhole Cuon alpinus habitat suitability. Surprisingly, masked palm civet Paguma larvata, Asian elephant, leopard cat Prionailurus bengalensis, common leopard Panthera pardus and sambar deer Rusa unicolor displayed a negative relationship with forest-related variables ( Table 2). The density of rivers in a local landscape had significant relationships with tiger, serow and wild pig while having significant negative relationships with smaller carnivores like dhole and marbled cat Pardofelis marmorata.
Anthropogenic variables were generally negatively associated with most modelled species. Forest obligate species like serow, Asian elephant, marbled cat, tiger, dhole and Asiatic black bear Ursus thibetanus displayed strong negative response to anthropogenic factors. However, wild pig and sambar deer showed positive associations with anthropogenic variables, suggesting that they may not avoid and may even select areas of relatively higher human impact.
The response of mammals to topographical predictors was highly variable. For example, roughness, CTI and slope position were positively associated with habitat use of Asian golden cat, wild pig, yellow-throated marten, marbled cat, common leopard and dhole but negatively associated with habitat use of Asian elephant, Asiatic black bear, serow, gaur Bos gaurus, leopard cat and sambar deer (Table 2). Muntjac had the highest mean habitat suitability across Bhutan (0.53) while clouded leopard had the lowest (0.13, Table 2).

Effectiveness of existing PAN
Some protected areas, notably Jigme Singye Wangchuck National Park (JSWNP) and Royal Manas National Park (RMNP) had relatively high effectiveness based on both mean and total relative

Simulation of areas of high conservation importance
We present the simulated optimal protected area selection rasters for each species in the supplementary material (Figures S20). As exemplars, given space limitations, we present only two here, in addition to the weighted multiple species selection raster (Fig. 4).

Conservation areas for dhole and Asian black bear
Dhole (or also known as Asiatic wild dog) had a high proportion of potential additional optimal conservation areas adjoining the protected areas JKSNR, JSWNP, PNP and BWS, and biological corridors connecting JSWNP and JDNP, JSWNP and WCNP and PWS and RMNP (Fig. 5a). Similarly, for Asiatic black bear, a large proportion of optimal additional conservation areas lie adjoining JKSNR, JSWNP, PNP and BWS and biological corridors connecting JSWNP and JDNP, RMNP and PNP and RMNP and JSWNP, JSWNP and WCNP and JWS and SWS (Fig. 5b).

Relative similarity and overlap of optimal cells between species
The full pair-wise correlation matrix (Fig. 6a) shows that, for several species, their optimal additional protected areas intersect relatively poorly with those of most other species and for the multi-species scenarios (e.g. wild pig, masked palm civet). It also shows AUC Area under the receiver operating characteristic curve several that appear to be relatively strongly associated with optimal new conservation areas for multiple other species. For example, there is a high correlation between tiger-dhole, tiger-marbled cat, tiger-clouded leopard, leopard cat-gaur, common leopard-dhole, common leopard-marbled cat.

Efficacy of multi-species prioritization strategies in providing optimal protection for each species
There were large differences in the performance of the two multi-species selection scenarios in terms of their overall efficacy of providing optimal selection for each species, and across species (Fig. 6b).
The unweighted sum of habitat suitability scenario had generally low performance, as measured by the correlation between the optimally selected areas for protection and predicted suitability across species (Fig. 7a). In particular, the unweighted analysis performed very poorly in protecting the species that had the largest overall umbrella performance (marbled cat, tiger, clouded leopard, common leopard, guar). It had a higher performance for civet, muntjac, black bear, which are notably among the weaker umbrella species and also relatively lower conservation priority.
In contrast, the conservation-priority-weighted multi-species scenario performed much better overall, based on average correlation across individual species (0.53 vs 0.26), and performed particularly well relative to the most effective umbrella species (Fig. 7b, marbled cat, tiger, clouded leopard, common leopard and gaur) and the species with the highest conservation importance (Table S2).

Umbrella species effectiveness
Most ungulate species emerged as effective umbrella species based on the mean percentage of total suitability (Fig. 8a). Muntjac was the most effective umbrella followed by wild pig, serow and sambar. Amongst the carnivores, Asiatic golden cat was relatively strong umbrella species. On the contrary, gaur was the least effective umbrella species probably due to small and restricted geographic range which may not offer wider habitat protection like the more ubiquitous muntjac or wild pig. Indicator analysis based on cluster dendrogram yielded five distinct clusters (Fig. 9). Mean relative richness (Mean) = average of habitat suitability pixels of each species Per cent habitat suitability (%) = sum of pixels with values < 0 in each PA or BC / total sum of suitable pixels across all PAs and BCs * 100

Indicator species effectiveness
The multivariate hierarchical cluster analysis identified five species groups based on similarity in their spatial patterns of habitat suitability (correlation of habitat suitability predictions, Fig. 9 & Fig. 8b). Cluster 1 consisted of three carnivores (marbled cat, clouded leopard and yellow-throated marten) and one omnivore (Asian black bear). Cluster 2 had two large ungulates (Asian elephant and gaur) and one omnivore (Masked palm civet). Cluster 3 was the largest group with five species: three of which were ungulates (muntjac, wild pig and sambar) and the other two were felids (leopard and leopard cat). Tiger was the singleton species in cluster 4. Finally, cluster 5 consisted of two medium-sized carnivores (Asian golden cat and dhole) and one ungulate (Himalayan serow). Asiatic black bear and yellow-throated marten in cluster 1 were paired more similarly together than marbled cat and clouded leopard and occur widely in Bhutan. Therefore, we would suggest one of these species (Asiatic black bear or yellow-throated marten) as a good choice of indicator for species in that cluster. In cluster 2, Asian elephant and gaur were clustered more closely than masked palm civet. The range of the former two is smaller than the latter but they represent better indicator species for lowland forest ecosystems and megafauna that depend on them. In cluster 3, besides common leopard, the four species (muntjac, wild pig, sambar and leopard cat) were closely clustered and are fairly ubiquitous. Thus, they represent a strong indicator of each other while common leopard has a relatively weaker association with the full cluster. Therefore, we suggest common leopard as a focal species, but sambar as an indicator species for these lowland forest associated species. Tiger, in the fourth cluster, is the sole indicator species. In the fifth cluster, all three species join near each other and hence they are equally good indicator species. However, from a conservation point of view, dholes are a better indicator, given their conservation interest, ecosystem role as an apex predator, and their association with a wide range of habitat.
Muntjac had the highest mean habitat suitability correlation with all other species (0.408) serving as a top overall indicator species across all species groups. Leopard cat (0.383) and masked palm civet (0.363) also displayed high mean correlation thus serving as important carnivore and omnivore indicators, respectively.

Discussion
Our analyses highlight the importance of multiscale assessment in identifying predictors influencing species distribution. In general, habitat suitability was positively influenced by environmental variables while anthropogenic variables had an overall negative influence. We assessed the effectiveness of the existing protected area network (PAN) and using simulations, identified potential conservation areas outside the existing PAN. Forests outside PAN across the south-central region have a high potential for terrestrial mammal conservation. Our indicator and umbrella species assessments suggest that multiple common and widely distributed species could offer broader protection than a single apex predator.

Species-habitat relationships
Multi-scale species-habitat modeling revealed that our study species, in general, were positively associated with protected areas and forest cover, and negatively related to anthropogenic variables. Most of these relationships were strongest at relatively broad spatial scales. The selection of variables at a broad scale suggests the importance of preserving large landscapes under legal protection. The broad-scale negative association with anthropogenic variables, Optimization of priority conservation areas (red polygon) for dhole a and Asian black bear b identified through simulation based on habitat quality, proximity to existing protected area and connectivity across the population Fig. 6 a Pixel-wise correlations between values in the output protected area selection raster (value 0-50) between each pair of species and the two multi-species scenarios. The height of the bar is proportional to the strength of the correlation. This shows the strength of the match between the predicted optimal protected areas between each pair of species. A high correlation means that a pair of species has high overlap between areas of a high frequency of selection in the simulation of optimal new protected areas. A high correlation means that areas that are optimal for protection for one species are optimal for protection of the other as well. b Frequency of different strengths of pixel-wise correlation between the optimal protected area selection raster for each species with the 15 other species. This figure shows the results of Fig. 6a as histograms, ranked from highest frequency of high correlation to lowest. A species with many high correlations of optimal additional protected area with other species optimal protected areas means that protecting optimal areas for that species would also protect optimal areas for many other species Fig. 7 a Mean pixel-wise correlation between the optimal protected area rasters for each species and all 15 of the other species. This gives a single value of the relative effectiveness of protecting optimal additional areas for each species in terms of its ability to protect optimal areas for all other species simultaneously. A species with a high average correlation between the optimal protected area rasters of all other species is a good umbrella species whose protection provides effective protection of many other species as well. b Pixel-wise correlation between the value of the optimal protected area selection raster for each species and the two multi-species selection scenarios. This figure shows the effectiveness of the two multi-species conservation scenarios (weighted and unweighted) in terms of their ability to protect the optimal unprotected areas for each species conversely, indicates that human land uses can have effects that extend beyond the immediate footprint of these activities, up to 16 km in our models. Our findings suggest that terrestrial wildlife species in Bhutan often require extensive habitat and are often sensitive to human disturbance at relatively broad scales (Karanth and Nichols 2017). In terms of management, this highlights the importance of protecting large extents of natural ecosystems with low human impacts for species that are sensitive to anthropogenic disturbance and which have large home range requirements.
Forest-related variables had a strong positive effect on 50% of the modelled species (n = 8, Table 2), highlighting the importance of protecting forest which provides habitat for native wildlife in Bhutan. A similar finding was reported in a more diverse composition of 45 mammals showing a strong positive affinity for forest cover (Penjor et al. 2020). Topographical complexity seems to be positively related to habitat use of the majority of species we studied in Bhutan, as observed through the positive association between habitat use and topographical variables represented by slope and terrain roughness. The generally rugged topography of Bhutan is inaccessible to human settlement and unsuitable for agriculture farming. However, many areas that are potentially of the highest inherent ecological stability for wildlife are also facing the greatest impact by humans. These rough landscapes likely offer refugia from human disturbance for large mobile species such as tigers (Krishnamurthy et al. 2016;Reddy et al. 2019) as well as smaller meso-carnivores like yellow-throated marten (Vergara et al. 2016). Approximately 50% of the modelled species were negatively associated with agriculture, suggesting that a landscape perforated or fragmented by agriculture, in general, is less suitable for the occurrence of the wildlife species (Table 2). Although agriculture may offer herbivores access to high-quality forage, it is often at the cost of lethal conflict (Woodroffe et al. 2005). This is often even more serious for carnivores, which face a severe risk of fatal conflict in agricultural-village mosaic lands (Rostro-Garcia et al. 2016). While fragmentation and intermediate extents of an agricultural mosaic can sometimes increase habitat diversity, leading to higher densities and diversity of wildlife species (Cushman and McGarigal 2003), increasing fragmentation of natural ecosystems by agriculture often brings animals close to humans and thereby escalates conflict (Sukumar 2003). Past studies demonstrated that increasing extent and heterogeneity of agricultural land either contributed to elevated conflict between large felids and human (Rostro-García et al. 2016) or decreased occupancy of large ungulates in Bhutan (Penjor et al. 2020).
Biodiversity conservation in agricultural landscapes is costly and has produced mixed results (Green et al. 2005), raising the appeal of landsparing, where existing farmland is intensively managed and improved to 'spare' other wildlife habitats (Green et al. 2005;Feber et al., 2015). However, this tradeoff is context-dependent, and in places where remnant habitat in highly productive agricultural landscapes is essential for biodiversity, and/or the protected area is small in size compared to the species home range, it may be necessary to adopt a combination of land-sparing and land-sharing approaches (Johansson et al. 2016). Multi-species, quantitative prioritization of areas of highest biodiversity value, as presented here, provide an objective basis to evaluate the relative Fig. 9 Cluster dendrogram used for identifying indicator species. The five clusters (numbered 1:5) are based on the predictions from the generalized linear model fitted to binary data of 16 terrestrial mammals in Bhutan merits of land-sparing and land-sharing conservation strategies.

Effectiveness of the existing PAN in protecting species richness
Bhutan has the constitutional mandate of maintaining 60% of the total geographical land under forest and identifying and protecting key biodiversity areas outside the existing PAN. However, there is little scientific guidance currently available to guide the optimal selection of areas for designation as new protected or conservation areas. Our analyses show that the south-central region encompassing JSWNP and RMNP had high effectiveness in protecting mean and total relative species richness (Fig. 3). PWS and JWS were effective in protecting mean richness per unit area but less effective in protecting the total proportion of species richness. JSWNP and RMNP, together with Manas Tiger Reserve in India, form one of the largest contiguous protected landscape in Asia, and one of the important global tiger conservation landscapes (WWF 2015). The ecosystems within this landscape range from tropical grassland to snow-capped mountain, supporting a wide spectrum of wildlife. Further, JSWNP and RMNP have been recently accredited as CA|TS sites (conservation assured tiger standards). CA|TS certifies and credits important tiger conservation sites with high standards in protecting tigers, their prey and habitat (Pasha et al. 2018). Thus, our findings provide evidence that these PAs are not only important for large species but also medium-sized and yet smaller terrestrial mammals (Dorji et al. 2019) and these additional conservation areas could bolster the protected landscape from anthropogenic pressure and perhaps also foster the maintenance of the ecological integrity of this important wildlife conservation landscape in the eastern Himalayas.
Biological corridors also displayed high variability in terms of mean relative species richness and proportion of total cumulative relative species richness. A similar finding was reported by Dorji et al. (2019), who found that mammal diversity in biological corridors and non-protected areas were lower than in protected areas, confirming that PAs are critical to large mammal conservation in Bhutan. But we believe that this does not necessarily mean that the current BCs do not serve their purpose. We found that several designated BCs provided important habitat for species and biodiversity overall.

Optimal protection of species-rich areas outside PAN
By implementing a newly developed model to optimize the selection of new protected areas to maximize protection of quality habitat that is functionally connected to existing protected areas, we found several areas of conservation importance abutting the existing protected areas and biological corridors. Based on the spatial optimization to identify areas for protection outside of the current protected area network, we found that the conservation-priority weighted sum of habitat suitability scenario for multi-species presented an overall better picture of optimal protection than the unweighted sum of habitat suitability scenario. This indicates that prioritizing species weighting based on conservation status improves the solution for all species collectively, and not just for those of highest conservation risk. Such weighted prioritization of species is vital to provide a better picture of conservation in an increasingly changing landscape (Saetersdal and Gjerde 2011). Several optimal areas adjoining the existing PAN were identified through simulation. For example, subtropical forests were identified abutting RMNP, PWS and JWS and BCs 3, 5 and 8. These areas are minimally inhabited, or not at all in some cases, and could offer important refuge to a wide spectrum of species (indicated by the weighted sum of multispecies richness). The weighted multispecies scenario performed well for species which could be considered as conservation focal species or umbrella species, such as tiger, common leopard, clouded leopard, gaur and elephant. These species are declining rapidly globally and are increasingly at risk of extinction even though a large portion of conservation funding is directed to securing their future. Our analysis further strengthens the evidence that conservation of the habitat of these species would help provide overall protection for a wide range of species under their umbrella. It is relatively easy for conservation managers to garner support for species which are inherently respected and popular in history, religion and culture (Macdonald et al. , 2017. Tiger and elephant, for example, are widely represented in Buddhist and Hindu mythology and although these species come into conflict with humans, they are better tolerated than other species such as wild pig or Asiatic black bear (Sukumar 2003;Kusi et al. 2020). This makes them more effective as an ambassador (Macdonald et al. 2017) or flagship species for conservation.
Based on the hotspots of relative species richness, we identified several focal areas outside of the existing PAN (Fig. 3 P1-P7). These are (i) subtropical and cool broadleaved forest in Samtse (Fig. 3 P1) adjoining JKSNR, (ii) cool broadleaved and temperate forest in Zhemgang (Fig. 33 P3) contiguous to PNP and BC4 (connecting PNP, JSWNP and RMNP), (iii) subtropical broadleaved forest in Pemagatshel (Fig. 3 P4) adjoining BC5 which connects RMNP to JWS, (iv) subtropical and cool broadleaved forest in Chhukha (Fig. 3 P2) forming contiguous forest corridor with Buxa Tiger Reserve in India, (v) subtropical and cool broadleaved forest adjoining JWS (Fig. 3 P5 & P6) and (vi) cool broadleaved forest adjoining BWS (Fig. 3 P7) to the east. Further, we identified potential corridors (Fig. 3 C1-C6) which connects the potential areas (P1-P7) and to the existing PAN. These focal areas could be high priority candidates for conservation outside protected area network because of their high predicted species richness and could serve as habitat for species while also bolstering the current PAs against the growing anthropogenic impacts.

Identification of umbrella and indicator species
Muntjac, wild pig, Himalayan serow, sambar and Asian golden cat emerged as the most effective umbrella species for optimal protected area designation. It might be surprising that an ungulate was identified as the most effective umbrella species in a landscape containing powerful and charismatic predators like tiger and common leopard. However, our results, based on a simulation study, suggest that trophic level and body size alone are not necessarily the best indicators of umbrella species effectiveness (Cushman and Landguth 2012). Common species may often serve as a better umbrella due to their adaptability and likelihood of surviving dynamic and hostile environments (Berger 1997). It would be sensible to design conservation programs to protect areas using more than one umbrella species. Our results suggest the selection of four ungulates and one medium-sized carnivore as umbrella species, which have a broader habitat niche, more extensive occupied habitat and the higher probability of long-term survival than the rare apex carnivores (Caro 2003). Lastly, for an organism at the higher trophic level to survive (e.g., the apex predators), healthy prey populations are essential. Designating ungulate species, which are generally much more widespread and abundant than apex predators, as umbrella species ensures that top carnivores like tiger, leopard and dhole are provided with sufficient prey populations, which enhances their reproductive potential and ultimately regulating food chain, providing check and balance to the ecological food web.
Our indicator analysis revealed findings that will benefit management as well as enhance our understanding of the choice of indicator species. The five clusters are the five groups of species that are most internally similar in terms of their patterns of habitat quality across Bhutan. Species in cluster 1, marbled cat, clouded leopard, Asiatic black bear and yellowthroated marten, are associated with highly forested upland habitats. One of these species could be chosen as an indicator to monitor the condition of this group of species and other species associated with the same forested upland environment. Cluster 2, in contrast, consists of three species, including two herbivorous megafaunas, elephant and gaur. Given the high conservation importance of elephant and guar, it is likely that future monitoring will be employed to track the populations and distributions of both species. Their close clustering, however, shows they inhabit very similar low land, flat forested habitats. The third cluster contains muntjac, wild pig, common leopard, leopard cat and sambar. This species group is associated with low elevation river valley ecosystems and lower slopes of foothill landscapes, often in mosaics of agricultural and forest lands. It is likely that common leopard will be individually the focus of conservation given its high charisma (Macdonald et al. 2017). Our results suggest that sambar would be a good choice as an additional indicator of the habitat condition of valley mosaic habitats and the species associated with them. Tiger is the sole species in cluster 4, suggesting that its habitat characteristics are most different from all of the other species. It, therefore, is a poor indicator of other species; but given its exceptional high conservation profile, it likely will be the focus of conservation planning in its own right. It is important, however, to note that based on our results conservation planning optimized for tiger is not likely to be optimal for most other species. This emphasizes the importance of optimizing multispecies conservation planning. The fifth cluster consists of golden cat, dhole and serow. These species are associated with upland ecosystems with relatively low human footprints and high forest cover.
Conservation planning based on indicator and umbrella species is controversial (Carignan and Villard 2002). On one hand, focus on umbrella/indicator species can make conservation planning simpler, monitoring more tractable and communication with public stakeholders easier, which helps managers leverage better community support (Macdonald et al. 2012). On the other hand, many of the species identified as the best indicators and umbrella species in our analysis are conflict-prone species, which can reduce public support for conservation planning built around them.
Conservation should strive to a strike balance between species protection and human wellbeing. Bhutan has a largely agrarian population, subsisting on small-scale farming and livestock rearing. Crop and livestock depredation by wild animals results in significant loss of annual household income in Bhutan (Wang and Macdonald 2006;Sangay and Vernes, 2008) leading to increasingly aversive sentiment towards wildlife (Woodroffe et al. 2005). Therefore, alternatives that promote coexistence should be integrated into conservation management plans. Community participation in conservation is increasingly seen as an essential approach to wildlife conservation and management (Andrade and Rhodes 2012). A joint management approach involving communities could be initiated in areas of high conservation importance outside protected areas, offering stewardship and regulating resource extraction by providing alternatives such as ecotourism.

Scope and limitations
Although the current PAN in Bhutan protects a wide array of ecosystems, it omits some areas of critical biodiversity value that are increasingly at risk due to expanding infrastructure, hydropower and other developments. While our study provides a comprehensive multi-species and multi-scale approach for identifying important habitat and conservation areas, future work should improve and enhance the robustness of the models by including human dimensions of conservation, such as the threat to wildlife arising from direct persecution, human-wildlife conflict and ancillary pressure on wildlife habitats from expanding human settlement. Besides, our assessment of BCs is based on their habitat suitability and not their functional performance in facilitating movement. Future work could use telemetry (e.g. Cushman et al. 2016) or landscape genetics (e.g. Cushman et al. 2006) to demonstrate the effectiveness of the BC network for facilitating movement and gene flow (Zeller et al. 2018). We also acknowledge that the current habitat suitability analyses do not account for detection heterogeneity arising from the non-detection of species (MacKenzie et al. 2002). Non-detection does not imply absence, but it could mean that the species was not present during our survey or the survey effort was inadequate to detect it. It is difficult to distinguish the observation process from the latent ecological process when detection probability is unaccounted. Models that cannot distinguish between state process (or ecological process, e.g., site occupancy) and observation process (i.e., variability in sampling process) often confuse between the true state of the animal and its absence due to one-time survey of a site and hence results in an underestimation of occupancy. Confusion between these processes risks spurious results, especially when considering species that are both cryptic and of high conservation concern. Temporal replication, therefore to estimate detection probability, improves the estimates of species diversity for a site (Iknayan et al. 2014). Hence, we caution readers to interpret our habitat suitability only as a probability of relative use (apparent) and not strictly as true occupancy. Finally, the low AUC for some species (< 0.75; Table 2) cautions us that for these species the prediction accuracy is only slightly better than by chance. However, AUC cannot be considered as a sole measure of a model's predictive power because it is in itself inappropriate for presence/absence data (Lobo et al. 2008).

Conclusions
We urge policymakers and conservation managers to broaden the scope of conservation by exploring and prioritizing alternative conservation areas based on multi-species optimization that combines both habitat suitability and network connectivity, as we present here. This study highlights the need to protect optimally located species-rich areas outside the current PAN (Tyrrell et al. 2020) and to reduce direct humanwildlife conflict through intensive anti-poaching patrols and law enforcement and community participation. Taking proactive measures such as optimally identifying and expanding protection would not only help fulfil global biodiversity targets (e.g., Aichi Targets) but also set new standards for science-based conservation and replenish the motivation for conservation amid escalating wilderness loss to human pressure.