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 (Li et al., 2020). 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 (Li et al., 2020; 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., 2014, 2016; Sharma 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 (Thinley et al., 2016). 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 species-specific 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.

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² (18%) of Bhutan, encompassing 3,647 km² (9%) of highly suitable, 2,681 km² (7%) of moderately suitable, and 878 km² (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² (40%) of suitable and 9,995 km² (60%) of unsuitable habitats, whereas biological corridors encompassed 248 km² (8%) of suitable and 3,086 km² (92%) of unsuitable habitats. In total, only 555 km² (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²) followed by WCNP (2,345 km²), BWS (648 km²), Paro (460 km²), and JKSNR (242 km²) (Table 2). The two largest highly suitable snow leopard 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).

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 (Shrestha et al., 2018), 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². This area is about 3,000 km² 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²; Leki et al., 2018). Although WCNP has a lower bharal density [1.8–2.4 individuals per km²; 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 (Lham et al., 2016), 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).

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.

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

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.

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.

SW was CEO of Natural Resources Development Ltd. while contributing to the manuscript. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.