Modeling Distribution and Habitat Suitability for the Snow Leopard in Bhutan

The snow leopard (Panthera uncia) is one of the world's most elusive felids. In Bhutan, which is one of the 12 countries where the species still persists, reliable information on its distribution and habitat suitability is lacking, thus impeding effective conservation planning for the species. To fill this knowledge gap, we created a country-wide species distribution model using “presence-only” data from 420 snow leopard occurrences (345 from a sign survey and 77 from a camera-trapping survey) and 12 environmental covariates consisting of biophysical and anthropogenic factors. We analyzed the data in an ensemble model framework which combines the outputs from several species distribution models. To assess the adequacy of Bhutan's network of protected areas and their potential contribution toward the conservation of the species, we overlaid the output of the ensemble model on the spatial layers of protected areas and biological corridors. The ensemble model identified 7,206 km2 of Bhutan as suitable for the snow leopard: 3,647 km2 as highly suitable, 2,681 km2 as moderately suitable, and 878 km2 as marginally suitable. Forty percent of the total suitable habitat consisted of protected areas and a further 8% of biological corridors. These suitable habitats were characterized by a mean livestock density of 1.3 individuals per hectare, and a mean slope of 25°; they closely match the distribution of the snow leopard's main wild prey, the bharal (Pseudois nayaur). Our study shows that Bhutan's northern protected areas are a centre for snow leopard conservation both at the national and regional scale.


INTRODUCTION
Understanding species distribution and habitat suitability are crucial for species management (Aryal et al., 2016;Bai et al., 2018). Managing species is increasingly important due to the growing anthropogenic pressure on natural habitats and the impending threats from climate change (Forest et al., 2012;Aryal et al., 2016;Watts et al., 2019). The latter is expected to be most severe on species that inhabit mountain ecosystems (Forest et al., 2012;Aryal et al., 2016); a prime example of such a species is the snow leopard (Panthera uncia) (Ghoshal et al., 2019). The snow leopard occurs in 12 countries: Afghanistan, Bhutan, China, India, Kazakhstan, the Kyrgyz Republic, Mongolia, Nepal, Pakistan, Russia, Tajikistan, and Uzbekistan . The main threats to its persistence are habitat loss and fragmentation, natural prey depletion, poaching, and retaliatory killing (Watts et al., 2019;Li et al., 2020). Although the population estimates are yet to be updated, the global snow leopard population is currently estimated at 2,710-3,386 mature individuals (McCarthy et al., 2017).
Snow leopards occur at low densities, which in concert with their elusive nature renders distribution and habitat suitability studies difficult (Ghoshal et al., 2019;Kalashnikova et al., 2019;Watts et al., 2019). Even though studies on snow leopard distribution exist, most of them were conducted in high-quality habitats (Suryawanshi et al., 2019) and are restricted to a few countries such as Pakistan (Hammed et al., 2020), India (Watts et al., 2019;Singh et al., 2020), Nepal (Aryal et al., 2016;Shrestha and Kindlmann, 2020), China (Bai et al., 2018;Li et al., 2020), Southern Russia (Kalashnikova et al., 2019), and Kazakhstan (Holt et al., 2018), whereas in the case of Bhutan, only little is known. Such imbalances in the knowledge of the snow leopard's distribution can result in an incomplete understanding of its distribution nationally and globally. From a conservation perspective, understanding the species' distribution as well as the protection status and connectivity of suitable habitats is important, because information on whether a suitable habitat falls inside or outside a protected area provides insight into habitat protection and helps in identifying areas that require conservation focus Thinley et al., 2021).
Based on the knowledge of the snow leopard's ecology, its distribution or habitat suitability are expected to be influenced by several biophysical and anthropogenic factors such as wild prey availability and livestock density, habitat type, elevation, slope, and aspect (Bagchi and Mishra, 2006;Aryal et al., 2014Aryal et al., , 2016Sharma et al., 2015;Alexander et al., 2016a). In the Himalayas, snow leopards are known to occur between 3,000 and 5,800 meters above sea level (a.s.l.) (Fox and Chundawat, 2016;Chetri et al., 2017) and on predominantly south-facing slopes of 24 to 40 • (Jackson and Hunter, 1996;Fox and Chundawat, 2016). Their main wild prey in the Himalayas is the bharal (Pseudois nayaur), commonly known as the blue sheep (Alexander et al., 2016b;Aryal et al., 2016;Holt et al., 2018;Filla et al., 2020). Bharal distribution has a positive effect on snow leopard distribution in Nepal (Aryal et al., 2016). Although wild prey is a key determinant of snow leopard distribution, dependence on livestock cannot be ignored, because snow leopard signs (e.g., scats and pugmarks) in Bhutan can be found in alpine meadows where livestock and bharal graze together . Snow leopards are known to depredate livestock such as yaks, horses, goats, and sheep (Chetri et al., 2017;Khatoon et al., 2017;Jamtsho and Katel, 2019;Lham et al., 2021), suggesting livestock can supplement snow leopard diet in areas with low density of natural prey.
In addition to biophysical factors, bioclimatic factors can also affect snow leopard distribution (Aryal et al., 2016;Bai et al., 2018;Hammed et al., 2020;Shrestha and Kindlmann, 2020). Hence, accounting for them can improve the accuracy of species distribution models (SDMs) (Aryal et al., 2016). Some of these effects can be indirect. Increasing temperature and precipitation, for example, can lead to an upward shift of the treeline, which, in turn, will reduce the areas of alpine meadows, the foraging grounds of bharals, and increase interspecific competition between livestock and bharal (Shrestha and Wegge, 2008;Forest et al., 2012;Shen, 2020;Shrestha and Kindlmann, 2020;Yang et al., 2021). Such a situation may affect the snow leopard's distribution and its habitat.
Species distribution models (SDMs) approximate speciesspecific habitat suitability and can be effective tools for conservation management and planning (Hammed et al., 2020;Thinley et al., 2021). They are often referred to as habitat suitability models (HSMs) or ecological niche models (ENMs) (Barbet-Massin et al., 2012;Atzeni et al., 2020;Hammed et al., 2020;Thinley et al., 2021). Based on the relationship between species occurrences and environmental factors, SDMs essentially assign suitability values to an area to predict its suitability for a target species (Bai et al., 2018;Watts et al., 2019;Thinley et al., 2021). Despite the availability of a wide range of SDMs, choosing an appropriate model remains a challenge because models differ in assumptions, data requirements, and predictive accuracy (Thuiller et al., 2009;Grenouillet et al., 2011). A commonly used SDM is the maximum entropy model (MaxEnt) (Phillips et al., 2006). Although MaxEnt is widely used due to its simplicity, user-friendliness, and low needs for computational power, it also has been criticized for model simplification and model overfit (Kaky et al., 2020). The accuracy and predictive uncertainty from individual models such as MaxEnt can be improved by using an ensemble of models to estimate average predictions from multiple models (Thuiller et al., 2009;Grenouillet et al., 2011;Ahmad et al., 2020;Kaky et al., 2020).
In this study, we investigate snow leopard distribution in Bhutan, identify suitable habitats, and assess the adequacy of protected areas in securing protection of suitable habitats. To achieve our goals, we built an ensemble of SDMs to predict the snow leopard's distribution and understand the effects of the underlying biophysical, bioclimatic, and anthropogenic factors. We also spatially overlaid the layer consisting of protected areas and biological corridors on to the habitat suitability layer to assess the adequacy of the snow leopard's protection in terms of areas encompassed by protected areas and connectivity in terms of areas occurring within the biological corridors. We expect that snow leopards occur mainly at high elevations inside the network of protected areas and biological corridors. We further expect that the distribution of the snow leopard follows those of the bharal and livestock. Increase in temperature and precipitation are further expected to negatively affect the habitat suitability for snow leopards.

Study Area
Bhutan is a landlocked country in the eastern Himalayas, sharing borders with China to the north and India to the south, west, and east. It has a geographic area of 38,394 km 2 [Forest Resources Management Division (FRMD), 2016], divided into 20 dzongkhags or districts (Figure 1). About half of the country (51.44%) is designated as a protected areas network, which includes five national parks, four wildlife sanctuaries, one strict nature reserve, and eight biological corridors (Lham et al., 2018). Bhutan is predominantly a mountainous country with rugged terrain and deep river valleys (Tshering et al., 2020;Lham et al., 2021). Owing to an elevation gradient from 97 m a.s.l. in the south to 7,750 m a.s.l. in the north (Tshering et al., 2020), Bhutan has high floral and faunal diversities (Tempa et al., 2019). The latter comprises 129 mammalian species, including top predators such as the tiger (Panthera tigris), the common leopard (Panthera pardus), the dhole (Cuon alpinus), and the snow leopard  (Figure 1) . Camera trapping studies conducted after this survey revealed snow leopard occurrences also in the Bumdeling Wildlife Sanctuary (BWS; personal communication, Tshering Dendup) and the Jigme Singye Wangchuck National Park (JSWNP; Letro et al., 2021). People residing in the northern region of these areas are mainly pastoralists (Leki et al., 2018) raising mostly yaks, cattle, horses, and sheep. Yaks are led to alpine meadows (>3,500 m a.s.l.) in the summer and to lower elevations (around 3,000 m a.s.l.) in the winter (Lham et al., 2021).

Snow Leopard Occurrence Data
Snow leopard occurrence data were collected during the national snow leopard sign survey from August 2014 to July 2015, and the camera trap survey from October to December 2015 Thinley et al., 2016). The sign survey and the camera trap survey were conducted in 395 and 227 survey grids, respectively. Each grid measured 4 × 4 km 2 . Prior to the sign survey, 100 members of forest staff were trained to reliably identify and differentiate signs of snow leopards from signs of other sympatric carnivores such as the Tibetan wolf (Canis lupus chanco) (Lham et al., 2021). Snow leopard scats were identified based on scat size, shape, and the presence of scrapes or pugmarks (Lham et al., 2021). Before the camera trap survey, a team of 150 personnel were trained to setup, monitor, and retrieve data from camera traps. Within each surveyed grid, search effort focused on areas that had the highest probability of snow leopard presence, such as ridgeline saddles, bases of cliffs, confluences of ridgelines near streams or rivers, and well-defined trails (Jackson and Hunter, 1996;Jackson et al., 2005;Thinley et al., 2015). The snow leopard signs included sighting, track, scat, and scrape (Devkota et al., 2013;Lovari et al., 2013). The sign search was maximized during dry periods when signs were most durable and easiest to detect (Thinley et al., 2015). A total of 227 camera trap locations were selected based on the results from the snow leopard sign survey, FIGURE 2 | Snow leopard occurrence records (black dots) in Bhutan, overlaid on protected areas (green), biological corridors (blue), and dzongkhags (districts, black outlines). and two cameras facing animal trails but separated by at least 1 km (Bai et al., 2018) were placed in each grid. During camera monitoring depleted batteries were changed, memory cards were emptied or replaced, misaligned cameras were repositioned, and lost or non-functional traps were replaced .
A total of 420 snow leopard occurrence records were obtained, 345 from the sign survey and 77 from the camera trap survey (Figure 2). Due to lack of adequate snow leopard images, the 77 occurrence records could not be identified to the individual level. We spatially filtered the occurrence records to maximally one occurrence per 1 × 1 km 2 grid (Boria et al., 2014;de Oliveira et al., 2014;Aryal et al., 2016;Holt et al., 2018;Gaspard et al., 2019;Watts et al., 2019). This grid size corresponds to the resolution of the bioclimatic data from the World Climate Database (www.worldclim.org, version 2.1; Fick and Hijmans, 2017). Consequently, we used 233 snow leopard occurrence records in our ensemble model. The data preparation and analyses were performed in R version 4.0.4 (R Core Team, 2020).

Biophysical, Bioclimatic, and Anthropogenic Factors
To model and predict snow leopard distribution and habitat suitability, we selected 12 environmental covariates for modeling: bharal distribution, habitat type, elevation, slope, aspect, livestock density, annual mean temperature (labeled as BIO1 in the World Climate Database, https://www.worldclim.org/ data/bioclim.html), mean monthly temperature range (BIO2), maximum temperature of the warmest month (BIO5), annual precipitation (BIO12), precipitation of the driest month (BIO14), and precipitation of the driest quarter (BIO17). The six bioclimatic factors were selected because they influenced snow leopard distribution in Nepal, Pakistan, and China (Aryal et al., 2016;Bai et al., 2018;Hammed et al., 2020;Shrestha and Kindlmann, 2020). These data were derived from the World Climate Database version 2.1, which comprises climate data for the period 1970-2000 at a resolution of 1 × 1 km 2 (Fick and Hijmans, 2017). We assigned to each snow leopard occurrence the bioclimatic values that corresponded to the occurrence location. Bharal distribution was modeled using occurrences recorded during the national snow leopard's prey sign and camera trap survey from 2014 to 2016 Thinley et al., 2016). We used modeled bharal distribution as a function of elevation, livestock density, and annual mean temperature (see Supplementary Figures 1, 2). We used a habitat layer composed of five habitat types: forest, meadow, rocky outcrop, scrub, and other (see Supplementary Figure 3; Lham et al., 2021). We then created a buffer with a 500 m radius around each occurrence record to align the spatial resolution of the bioclimatic data, 1 × 1 km 2 . We assigned a habitat type to each occurrence record that corresponded to the dominant habitat type inside each buffer. We computed the mean elevation, slope, and aspect for each occurrence buffer from the digital elevation raster data using the function get_elev_raster from the R package elevatr (Hollister et al., 2021). We calculated livestock density per gewog (the subdistricts of a dzongkhag) using the livestock abundance data from the 2015 Annual Livestock Statistics [Department of Livestock (DoL), 2015]. We assigned to each snow leopard occurrence the livestock density of the gewog in which the occurrence was located. Lastly, to remove collinear factors, we calculated Pearson's correlation coefficient (r), and considered any pairwise combination of factors collinear if |r| > 0.5 (Lham et al., 2021). After removing collinear factors, we used bharal distribution, livestock density, and slope to model snow leopard distribution. The six bioclimatic factors and elevation were highly collinear with bharal distribution and therefore not used in the modeling. The effects of bioclimatic factors, specifically annual mean temperature (BIO1) and elevation, was expressed through bharal distribution, which was modeled as a function of annual mean temperature (BIO1), elevation, and livestock density. Protected areas and biological corridors layers were obtained from the Forest Resources Management Division (FRMD), Department of Forests and Parks Services, Bhutan.

Species Distribution Modeling
We analyzed our data using multiple SDMs including MaxEnt, following an ensemble approach. To predict snow leopard distribution from the spatially filtered occurrence records, we formed an ensemble of SDMs using the R package BIOMOD2, version 3.4.6 (Thuiller et al., 2020). BIOMOD2 is a platform that can simultaneously run 10 types of SDMs: generalized linear model (GLM), generalized additive model (GAM), generalized boosting model (GBM), classification tree analysis (CTA), artificial neural network (ANN), surface range envelop (SRE), factorial discriminant analysis (FDA), multivariate adaptive regression spline (MARS), random forest (RF), and maximum entropy (MaxEnt) (Thuiller et al., 2020). In the ensemble model, we tested the 10 SDMs to predict snow leopard distribution. We randomly generated three times as many pseudo-absence data records as the occurrence records throughout Bhutan and used both datasets as input to the SDMs. We used 80% of the occurrence records to calibrate the models and 20% to test the models' predictive performance (Ahmad et al., 2020;Thuiller et al., 2020). This data split was repeated 10 times, resulting in 100 calibrated models. We then converted the cellwise predictions of the 100 models into binary values using a threshold to maximize sensitivity and specificity. Sensitivity and specificity measure the percentage of presence and absence correctly predicted (Grenouillet et al., 2011;Ahmad et al., 2020). We used true skills statistics (TSS) and the area under the receiver operating characteristic (ROC) curve to assess the effectiveness of each modeling technique (Thuiller et al., 2009;Ahmad et al., 2020). To form the ensemble model, we selected SDMs using TSS and ROC threshold values > 0.7 (Ahmad et al., 2020). To predict the snow leopard distribution in Bhutan, we then averaged the selected model outputs.
The resulting snow leopard distribution was imported to ArcGIS version 10.8.1. [Environmental Systems Research Institute (ESRI), 2020] and reclassified into four categories of occurrence probability: highly probable, moderately probable, marginally probable, or not probable. These probabilities correspond to the species' suitable habitat and realized niche, a subset of the fundamental niche (Phillips et al., 2006). We identified and characterized the surface as highly suitable (values ranging from 1.0 to 0.8), moderately suitable (0.7-0.4), marginally suitable (0.3-0.1), and unsuitable (0). Following Thinley et al. (2020), we overlaid protected area and biological corridor layers on the snow leopard suitability layer to assess how much of the area deemed as suitable was encompassed within the protected areas and how well-suitable habitats were connected. Finally, we calculated the importance of each variable using the function get_variables_importance (Thuiller et al., 2009). This function uses a random procedure that is independent of the modeling technique and enables direct variable comparison across models (Thuiller et al., 2009(Thuiller et al., , 2020. The procedure calculates a correlation coefficient between the fitted and predicted values whereby the covariate under evaluation is permuted randomly.

RESULTS
Except for three models (ANN, GAM, and MARS), all SDMs converged. Amongst the converged models, GBM and RF had the highest mean TSS and ROC values (0.81 and 0.93, respectively). CTA and MaxEnt had the lowest TSS value (0.79 each), and CTA had the lowest ROC value (0.89). The remaining models had mean TSS and ROC values > 0.80 and > 0.91, respectively, indicating good model performance. Therefore, GLM, GBM, CTA, FDA, RF, and MaxEnt were used to form the ensemble model.
Bharal distribution contributed the most to snow leopard habitat suitability (93%), followed by livestock density (4%) and slope (3%), all effects being positive. The average suitable habitat was characterized by the presence of bharal, had a mean livestock density of 1.3 individuals per hectare, and a mean slope of 25 • .
The ensemble model had mean TSS and ROC values of 0.81 and 0.95, respectively. The ensemble model predicted that suitable snow leopard habitat covered 7,206 km 2 (18%) of Bhutan, encompassing 3,647 km 2 (9%) of highly suitable, 2,681 km 2 (7%) of moderately suitable, and 878 km 2 (2%) of marginally suitable habitats (Figure 3; Table 1; the predicted mean distributions of the single SDMs are presented in the Supplementary Figure 4). When protected areas and biological corridors layers were overlaid on the snow leopard suitability layer, protected areas encompassed 6,403 km 2 (40%) of suitable and 9,995 km 2 (60%) of unsuitable habitats, whereas biological corridors encompassed 248 km 2 (8%) of suitable and 3,086 km 2 (92%) of unsuitable habitats. In total, only 555 km 2 (3%) of suitable habitats were predicted to be outside the network of protected areas and biological corridors.
Amongst the protected areas, the majority of suitable snow leopard habitat occured in JDNP (2,606 km 2 ) followed by WCNP (2,345 km 2 ), BWS (648 km 2 ), Paro (460 km 2 ), and JKSNR (242 km 2 ) ( Table 2). The two largest highly suitable snow leopard FIGURE 3 | Extent of unsuitable and marginally-to-highly suitable habitats for snow leopards in Bhutan, as predicted by an ensemble of Species Distribution Models. habitats occurred in JDNP and WCNP. JKSNR had the smallest area of highly suitable and unsuitable habitat amongst the protected areas, and PTFD had the largest area of unsuitable habitat ( Table 2).

DISCUSSION
We modeled the distribution of snow leopards in Bhutan using an ensemble model composed of six SDMs. Based on the ensemble model, suitable snow leopard habitat sums up to one fifth of the country's total area. Further, as expected, we found that the snow leopard is distributed mostly in the northern high-altitude region of Bhutan, and the majority of suitable areas fall inside the network of protected areas and biological corridors. Moreover, the suitability of an area for the snow leopard is largely explained by the presence of its main prey, the bharal, confirming previous studies (Aryal et al., 2014;Lyngdoh et al., 2014;Sharma et al., 2015;Alexander et al., 2016a,b;Suryawanshi et al., 2017;Holt et al., 2018;Bhattacharya et al., 2020;Singh et al., 2020;Lham et al., 2021). While bharal is by far the most important prey species of the snow leopard in Bhutan (Lham et al., 2021) and hence determines its distribution, as shown by our analyses, other prey species, such as Himalayan tahr (Hemitragus jemlahicus) in Nepal , may play a similar role in other regions, and their distribution should be considered.
Overall, Bhutan has a suitable snow leopard habitat of 7,206 km 2 . This area is about 3,000 km 2 larger than previously suggested by Forest et al. (2012). This discrepancy in suitable habitat may be due to the differences in factors used to identify suitable snow leopard habitat. Forest et al. (2012) identified suitable snow leopard habitat based on habitat types and topographical features such as ruggedness and elevation, while they did not account for the availability of bharal and livestock. The difference in the size of snow leopard suitable habitats from the two studies shows that accounting for ecological factors such as prey distribution may be crucial for a more accurate habitat modeling of a wild carnivore. Within Bhutan, the two largest protected areas, Jigme Dorji National Park (JDNP) and Wangchuck Centennial National Park (WCNP), encompass most of the highly suitable habitat. The high suitability of JDNP is explained by the presence of high bharal density (8.5-9.3 individuals per km 2 ; Leki et al., 2018). Although WCNP has a lower bharal density [1.8-2.4 individuals per km 2 ; Wangchuck Centennial National Park and Word Wildlife Fund (WCNP and WWF), 2016] than JDNP, it still has a relatively large area of highly suitable snow leopard habitat. This finding implies that WCNP may have the potential for future bharal population expansion, which would result in larger suitable snow leopard habitat. Paro Territorial Forest Division (PTFD) has the third largest highly suitable snow leopard habitat, which is about half the size of highly suitable habitat in JDNP. PTFD has bharal densities similar to WCNP (Lham et al., 2021), and given its size, PTFD is an important snow leopard conservation area. Bumdeling Wildlife Sanctuary (BWS) and Jigme Khesar Strict Nature Reserve (JKSNR) have the smallest areas of highly suitable habitat. For BWS, this finding is explained by the lower elevation of the sanctuary; bharals inhabit sites at ≥3,500 m a.s.l. (Lovari and Mishra, 2016). For JKSNR, the result is partly explained by JKSNR's small area, although JKSNR has high bharal densities similar to JDNP (Lham et al., 2021). In addition, aside from seasonal yak herding, JKSNR lacks human activities. This lack makes JKSNR not only an ideal snow leopard conservation area but also provides potential for expanding bharal populations. Overall, our results show that we need to understand the spatiotemporal trends in bharal distribution to be able to conserve suitable snow leopard habitats in Bhutan (Leki et al., 2018).
Suitable snow leopard habitat within Bhutan is contiguous from JKSNR in the west to WCNP in the central north. BWS in the north-east, however, is disconnected from WCNP. The gap likely consists of areas that are below 3,500 m a.s.l. or without bharal in the eastern part of WCNP [Wangchuck Centennial National Park and Word Wildlife Fund (WCNP and WWF), 2016]. This gap might explain why no snow leopards were caught on camera in BWS during the national snow leopard survey , and only one individual was detected in a subsequent survey (Tshering Dendup, personal communication). The most western biological corridor connecting JKSNR and JDNP encompasses only a portion of the suitable snow leopard habitat. Instead, the majority of suitable snow leopard habitat between JKSNR and JDNP lies inside the PTFD. Therefore, to ensure greater protection and landscape connectivity of the suitable habitats, we recommend enclosing this snow leopard habitat in the network of protected areas by either extending the biological corridor to the north or readjusting the boundaries of JKSNR or JDNP. Furthermore, three small, disconnected patches between 89 • 30 ′ 30 ′′ N and 91 • 31 ′ 30 ′′ N were predicted to be suitable. These patches are surrounded by unsuitable habitats and have no records of bharal. Nonetheless, a camera trap survey caught a snow leopard in the mid-patch, which falls inside Jigme Singye Wangchuck National Park (JSWNP) (Letro et al., 2021). JSWNP is connected to the nearest snow leopard suitable habitat, WCNP, in the north via a biological corridor, BC8 (Letro et al., 2021). The presence of high elevation (>3,500 m a.s.l.) ridgelines and the presence of ungulates such as goral (Naemorhedus goral), musk deer (Moschus leucogaster), and Himalayan serow (Capricornis thar) may have caused snow leopard movement from WCNP in the north to JSWNP in the south-central. Corridor BC8 connecting WCNP to JSWNP appears effective in ensuring snow leopard movement; yet, it is important to assess if the suitable habitats in JSWNP can sustain snow leopard sub-populations.
For the persistence of the snow leopard, it is important to ensure landscape connectivity at the regional level. Our results highlight the need for transboundary snow leopard conservation with China to the north and India to the west to maintain landscape connectivity. To the west of Bhutan, JKSNR serves as an important link connecting Bhutan to India through the Kanchenjunga landscape conservation development initiative (KLCDI). KLCDI is a regional conservation strategy between Bhutan, India, and Nepal (Oli et al., 2013). To the north and to the east of Bhutan, PTFD, JDNP, WCNP, and BWS are directly connected to the Tibetan Autonomous region of China. The Tibetan Autonomous region of China is an important snow leopard habitat in China. A well-connected habitat across Bhutan, China, and India will foster effective snow leopard movement and contribute to maintaining a viable snow leopard sub-population.
People in northern Bhutan are pastoralists rearing livestock (Lham et al., 2021). In the study area, livestock numbers are increasing (Thinley et al., 2014;Leki et al., 2018). This increase in livestock is an issue because livestock and bharal share resources in the alpine meadows, and livestock is gradually displacing bharal (Mishra et al., 2004;Bagchi and Mishra, 2006;Shrestha and Wegge, 2008;Bhattacharya et al., 2020). A reduced bharal availability may lead snow leopards to attack livestock, likely causing conflict with people (Lyngdoh et al., 2014;Suryawanshi et al., 2017;Holt et al., 2018), and to move to unsuitable habitats in search of alternative food (Letro et al., 2021). For the persistence of the snow leopard, it is therefore important to understand resource competition between livestock and bharal and to identify what must be done to ensure that a bharal population persists.
Among the three factors explaining the snow leopard distribution, slope explained the least. In this study, suitable snow leopard habitats were predicted to have an average slope of 25 • . The selection of 25 • slopes is due to the high presence of the bharal and the low human activities (Shrestha and Wegge, 2008). Further, the selection of slope by the bharal may depend on the season and the presence or absence of livestock in the area (Shrestha and Wegge, 2008;Filla et al., 2020). In India, suitable snow leopard habitat had an average slope of 24 • (Fox and Chundawat, 2016). In Kazakhstan and Nepal, slope was an important determinant of snow leopard prey such as argali (Ovis ammon) and urial (Ovis orientalis) (Chetri and Pokharel, 2005;Shrestha and Wegge, 2008;Fox and Chundawat, 2016;Holt et al., 2018;Filla et al., 2020). This indirect association of slope and prey further reiterates the importance of prey presence in a habitat suitable to snow leopards.
Although elevation and annual mean temperature (Aryal et al., 2016;Holt et al., 2018;Watts et al., 2019;Filla et al., 2020;Singh et al., 2020) were expected to influence snow leopard habitat suitability, our model did not pick up on these factors. Yet, they were selected as important determinants of bharal distribution, which, in turn, was selected as an important determinant of snow leopard distribution. Therefore, effects of elevation and annual mean temperature on snow leopard distribution may have been mediated through their effects on prey distribution.
In Bhutan, the snow leopard is a protected species [Forest and Nature Conservation Act (FNCA), 1995], and our study indicates that the majority of the suitable snow leopard habitats occur inside a network of protected areas. Moreover, these suitable habitats are mostly well-connected. Our study highlights the importance of maintaining regional connectivity with Nepal to the west through the KLCDI and the need to initiate transboundary conservation efforts with the Tibetan Autonomous region of China. Moreover, because bharal is of utmost importance to the long-term persistence of the snow leopard in Bhutan (Lyngdoh et al., 2014;Lham et al., 2021), we need to understand better the spatiotemporal trends in bharal distribution (Leki et al., 2018) and population dynamics. Due to the increasing presence of livestock in bharal habitat, it is also important to study resource competition between bharal and livestock. Understanding how suitable snow leopard habitat is influenced by these variables will allow park managers to design site-based actions to ensure the conservation of the snow leopard, and conserving the snow leopard will contribute to maintaining a healthy mountain ecosystem (Aryal et al., 2016;Filla et al., 2020).

DATA AVAILABILITY STATEMENT
Most data used in the analysis will be made available by the corresponding author. Due to the sensitive nature of some information, such as camera-trap coordinates and data coming from the national snow leopard camera trap survey of Bhutan, some data will not be available on public domain.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because there is no direct involvement of animals in this study.

AUTHOR CONTRIBUTIONS
DL: conceptualization, study design, data curation, formal analysis, and writing the original draft. GC and AO: analysis conceptualization, review, and editing. SS: review and editing. PT: conceptualization, study design, and review. NW and SW: conceptualization and study design. All authors contributed substantially to the article and approved it for submission.

FUNDING
The national snow leopard sign and camera trap surveys was funded by the WWF Bhutan (grant agreement no. CA 42), the Bhutan Trust Fund for Environmental Conservation (grant MB2015-16/05/0166), the Global Environment Facility (GEF-5 cycle), the World Bank through the Strengthening Regional Wildlife Conservation project (IDA-49830), Nature and Biodiversity Conservation Union (NABU), and the Royal Government of Bhutan. Further, we thank the Swiss Government and the University of Zurich (UZH) for financial support of DL's research at the UZH through a Swiss Government Excellence Scholarship for Foreign Students (2016/1113) and a Forschungskredit CanDoc (FK-20-092), respectively.

ACKNOWLEDGMENTS
We thank Mr. Chencho Norbu and Mr. Phento Tshering, the former directors of the Department of Forests and Parks Services, Ministry of Agriculture and Forests, Bhutan, for permission to use the Department's data (WCD/ICS/SLCP/2015-2016/4055 dated 07/04/2016). We also thank the forest staff involved in the national snow leopard sign and camera trap surveys.