The DNP (75°27′E, 35°00′N) is a 3,584 km² alpine plateau east of Nanga Parbat Peak, Northern Areas, Pakistan. Elevations range from 3,500 to 5,200 m above sea level, with about 60% of the area between 4,000 and 4,500 m above sea level (Nawaz et al., 2008; Figure 1 ).

Wetland regions and several tiny lakes are found on the DNP plateau. Their genesis is owing to snow and ice melting from the surrounding high mountains, as well as the existence of a permafrost table, which forms a thermokarst landscape (Karlsson et al., 2012). Deosai’s and surrounding highland habitat is known for its unpredictable, unstable, highly seasonal, and harsh environmental conditions. The typical daily temperatures range from -20°C to 12°C, and yearly precipitation ranges from 510 to 750 mm, largely as snow (Nawaz et al., 2008).

The vegetation in DNP is primarily herbaceous perennials, with a high number of cushion-forming and tufted species. Grasslands and areas with persistent snow or ice cover make up the majority of the land cover in high-elevation regions. Instead, the valley are mainly constituted by bogs and ponds with associated vegetation such as sedges, grasses, and plants like Aconitum violaceum, Swertia perfoliata, and Saxifraga hircus (Bellemain et al., 2007).

The latest results on bear population in DNP estimated ca. 58 – 77 bear specimens (Abbas et al., 2022). In addition to the bear, other mammalian species live in the DNP, such as the Himalayan red fox (Vulpes vulpes montana), the Mongolian wolf (Canis lupus chanco), the long-tailed marmot (Marmota caudata), the Himalayan ibex (Capra ibex sibrica), and 17 other smaller species (Nawaz, 2008; Ali H. et al., 2021). Many permanent human settlements exist inside and along the boundaries of the DNP, such as Sadpara, Shilla, Dappa, and Karabos. In addition to these permanent communities, it is common to find transhumance activity coming from the valleys (Nawaz, 2007).

Aside from the issues of overgrazing and soil degradation that affect DNP and its surrounding areas, the most pressing concern is the regions with water bodies and wetlands that are nourished by glacier/snow and permafrost melting throughout the summer. These regions, as well as their surrounds, are critical for flora and animals to live and thrive.

We collected bear occurrences by using camera traps in DNP and neighbouring areas during the field surveys from June 2021 to May 2022. This method was widely used to capture images of several bear species in their habitat (e.g., Ahmad et al., 2024; Lacy et al., 2024; Thakur et al., 2023; Fida et al., 2024). We used a motion-sensing digital camera equipped with infrared flash (HC500/PC900 model, Reconyx brand, Holmen, Wisconsin, USA) to capture bear pictures placed in 89 camera stations (e.g., Ahmad et al., 2024) with a capture rate value of 65.5%. Each camera was placed on a ground-mounted steel pole (about 50 cm tall), with a minimum distance of 1 km (mean = 1.824 km; maximum = 5.84 km) between camera stations (e.g., Gubbi et al., 2023; Wolfson et al., 2023; Melo-Dias et al., 2024). We used this distance value following a similar study carried out in the south of Iran with similar ecosystem conditions (Almasieh et al., 2016). The camera traps were positioned such that vegetation in the front could not activate them; they were designed to take pictures continuously at 1-second intervals when triggered (Bischof et al., 2014). The camera traps were normally left operational at each point for 10–40 days (mean = 25 days per station, Ahmad et al., 2022). To reduce the effect of spatial autocorrelation, we thinned all the species occurrences to ensure one presence record per environmental raster cell (e.g., Freeman et al., 2022). Finally we obtained 52 bear presence records ( Figure 1 ) and used them in the bear’s spatial analyses. Although previous publications have reported bear presence records for this area, we prefer to use only our records for two main reasons: 1) they are based on a more recent and reliable camera trap survey (e.g., Fida et al., 2024); and 2) they cover the whole park without any spatial autocorrelation.

We selected the environmental predictors based on the bear ecology obtained from other studies carried out in similar areas (Nawaz et al., 2014; Almasieh et al., 2016, 2019; Mohammadi et al., 2021a). The elevation map was derived from a 30-meter Digital Elevation Model (DEM) created by the Shuttle Radar Topography Mission (SRTM, http://earthexplorer.usgs.gov). We generated the slope maps starting from the DEM map by using the Surface tool in ArcGIS v. 10.3 (http://www.esri.com/software/arcgis). We considered 23-class global land cover maps (GlobCover version 2.3; ESA, 2009) for our analysis, but, after an initial screening, we selected only those most preferred by the bears, such as: 1) closed to open grasslands (>15%); 2) bare areas (23.5%); and 3) permanent snow or ice (22.6%). We used a moving window with a radius of 5 km was in the Focal Statistics tool to create density maps for the cover types selected (Mohammadi et al., 2021a; Almasieh et al., 2023a). All the following distances were calculated in the Euclidean space between a point in the map and the corresponding predictor. Distance from DNP borders was added in consideration of the current level of bear protection. Distance from rivers was used for water resources (Ansari and Ghoddousi, 2018), and distance from roads and human settlements was used for human disturbances. Distances to areas used by livestock for grazing and livestock movement routes represented both human and livestock presence. The coordinate system utilized was WGS 1984 UTM zone 43, and all variables were used at a 30-m resolution ( Supplementary Materials 2 Figures S1–S9 ).

We used the Variance Inflation Factor (VIF) to evaluate environmental predictor multicollinearity. VIF > 3 was set as the threshold (Zuur et al., 2010), and variables exceeding the threshold were excluded using the USDM package (Naimi et al., 2014) in R. 6.0 (R Core Team, 2019). All the predictors used in our modelling analysis were included in the Table 1 .

We explored the bear’s habitat suitability by using SDMs in an ensemble forecasting approach, as implemented in the R (v. 4.0.2) package “biomod2: Ensemble Platform for Species Distribution Modeling” (Thuiller et al., 2019; the R code is available in the Supplementary Material 1 ). We combined predictions from different model types and fitted multiple suitability models because it has been shown that this methodology improves model accuracy (Araújo and New, 2007; Almasieh et al., 2022; Mohammadi et al., 2022; Di Febbraro et al., 2023). We combined five algorithms into two main categories, such as: 1) the regression-based models—Generalized Linear Model (GLM) and Multivariate Adaptive Regression Splines (MARS); 2) the machine learning models—Maximum Entropy Model (MaxEnt), Random Forest (RF), and Generalized Boosting Model (GBM).

A 3:1 training-testing division was used to allocate occurrence points. A total of 500 pseudo-absence points were randomly generated outside a 1 km radius around occurrences. We chose this value of pseudo-absence beyond a radius of 1 km from occurrence since the ecosystem around the DNP was identical to that of the park itself (Barbet-Massin et al., 2012). Each model was run with 20 replications to increase reliability (Barbet-Massin et al., 2012; Kaboodvandpour et al., 2021). For each variable, we calculated the mean contribution to the model with the best performance by using biomod2, which also shows the response curves. We converted the continuous habitat suitability map generated by the ensemble model into a binary map through the 10th percentile threshold (Ahmadi et al., 2020; Almasieh et al., 2023b). This threshold (the value above which the model classifies correctly 90% of the training locations) is a conservative value that excludes outliers from the training presence points and is commonly used in species distribution modelling studies, especially when we have areas with similar environments (Ten Caten et al., 2023). Patches overlapping with the main presence records were assumed to be the core of the bear’s habitat suitability, and the total coverage by DNP was calculated.

The initial round of SDMs utilised all 10 models; then, we selected only models with an AUC > 0.9 and a True Statistic Skill (TSS) > 0.75 (Eskildsen et al., 2013; Almasieh and Cheraghi, 2022). As a threshold-independent measure of predictive accuracy, AUC provides an overall summary of the model’s diagnostic accuracy (Fielding and Bell, 1997). Following Swets (1988), the evaluation criteria for the AUC were classified as fail (0.5–0.6), poor (0.6–0.7), fair (0.7–0.8), good (0.8–0.9), and excellent (0.9–1). The TSS is a threshold-dependent measure that combines sensitivity (true positive rate) and specificity (true negative rate) to assess model performance. It is expressed as Sensitivity + Specificity – 1 (Allouche et al., 2006). The evaluation criteria for TSS were classified as poor (<0.4), good (0.4–0.75), and excellent (>0.75) (Landis and Koch, 1977).

Enhancing landscape connectivity is one of the most recommended anthropogenic pressure-adaptation strategies for protecting biodiversity as the global changes (e.g., Velázquez et al., 2017; Kabir et al., 2017; Velázquez et al., 2022; Tang et al., 2023). In order to generate a landscape connectivity map, we used the ensemble modelling output to produce a resistance layer following the methods proposed by Wan et al. (2019). This approach is based on the notion that animals select dispersal movement paths based on the same features that they use to select habitat. If this assumption is correct, habitat suitability would be the inverse of landscape resistance (i.e., areas with high habitat suitability would have low resistance to dispersal movement). The ensemble layer was rescaled in 0-1 mode by using the Rescale by Function tool in ArcGIS. A negative exponential function calculates resistance as follows: R=1000^{(-1×Ensemble Suitability)}, where R is the pixel resistance value (Mateo-Sánchez et al., 2015). Resistance values were linearly rescaled from 1 to 10, where 1 and 10 were minimum and maximum resistance, respectively (Wan et al., 2019).

We well know that there are numerous methods for evaluating landscape connectivity, including least-cost path analysis, resistant kernel models, and circuit-theory-models (Phillips et al., 2021). Of these methods, circuit-theory based models have shown consistently high performance in estimating connectivity, and there has been a notable uptick in their use in connectivity modelling (Phillips et al., 2021). Additionally, omnidirectional methods like Circuitscape demonstrated accurate estimations of connection corridors, particularly in situations where the animal uses random exploration of the underground —as is the case with bears during dispersion migrations, for instance (Zeller et al., 2020).

Electrical circuit theory models (McRae et al., 2008) simulate habitat connectivity based on the random walk approach and electrical circuits, with the current that imitates the organisms moving between focal nodes (occurrences or patches) based on voltage (movement probability) and resistance (permeability) (McRae et al., 2008; Roever et al., 2013). This approach identifies multiple potential movements, unlike least-cost methods, which identify only one route (Urban and Keitt, 2001; Urban et al., 2009), making it useful for modelling gene flow (McRae and Beier, 2007).

We carried out electrical circuit modelling using Circuitscape v4 (McRae and Shah, 2009) through an all-to-one approach that maps current density occurrences, consuming less RAM (McRae et al., 2008). The resistance map and bear occurrence records were used as inputs. We determined connections to eight neighbours per cell. Lower resistance indicates a higher movement density.

We found that five models (GLM, MaxEnt, RF, MARS, and GBM) showed high model validation values with AUC > 0.9 and TSS > 0.75, respectively. RF obtained the best performance ( Table 2 ). We detected that elevation and distance from DNP, livestock, and rivers were the most important variables for predicting the potential distribution of bears ( Table 3 ).

Low habitat suitability values increased with distance from DNP borders (up to ~10 km). The bear’s probability of presence was consistent about 2 km from rivers. Up to around 20 kilometers from roadways, bear habitat suitability appeared to improve. Distance from livestock regions decreased gradually, unlike distance from livestock transportation routes, which decreased quickly before flattening at around 20 km ( Figure 2 ). At least, high probabilities of bear presences were found between 3,750 and 4,250 meters above sea level.

According to the ensemble suitability maps, the bear’s habitat suitability was highest in the central region of the research area (inside DNP), while there were also other suitable patches in nearby areas ( Figure 3 ). These locations are crucial for encouraging bear migration over their whole range ( Figure 3 ).

We found that the habitat patches covered 1,125 km² of the entire study area, and ca. 73% of them were located inside the DNP. The core of the bear’s potential distribution fell mainly within DNP, with a total area of about 1,020 km² (13% of the study area; 80% coverage by DNP, Figure 4 ). We detected some other little patches in areas close to the border of DNP ( Figure 4 ). Our results showed that ca. 82% of livestock grazing areas occurred within habitat patches, and livestock grazing areas created about 11% of habitat patches for the bear. Total road and river densities within habitat patches were about 151 and 341 m/km², respectively. There was robust connectivity among occurrence points of the bear in the central part of the study area within DNP ( Figure 5 ), mainly due to the high number of bear records and low resistance values, i.e., a high value of probability of presence. Connectivity to several locations outside DNP close to the national park boundaries was discontinued ( Figure 5 ).

Along the altitudinal gradients, we found that the bears preferred an elevation between 3,500 and 4,500 m a.s.l. (e.g., Dar et al., 2021). Here, the possibility of sighting the bear increases due to the presence of meadows, bare patches, and prolonged snow and ice. Instead, in other areas, the bear preferred to avoid high-elevation areas because food availability is lower than in forested habitats at low elevations (e.g., valley), increasing the human-bear interactions (Martin et al., 2012; Almasieh et al., 2019; Mohammadi et al., 2021a; Kumar et al., 2022). In the valleys, uncontrolled overgrazing and activities such as beekeeping in the highlands can increase the risk of conflict with local communities (Sharief et al., 2020; Kumar et al., 2022).

Our results confirmed that one of Pakistan’s best bear-protected areas is DNP, where bear populations are growing due to effective conservation initiatives, as reported by Abbas et al. (2022). A few patches of bear habitats located outside the boundary of the park support the survival of the bear population, promoting the dispersion of this species also in other areas of its range (Nellemann et al., 2007; Martin et al., 2010; Fernández et al., 2012; Ziółkowska et al., 2016; Penteriani et al., 2020). Our results showed that bears select riparian and riverine environments because the water resource was among the most influential variables in predicting the bear’s potential distribution, as confirmed by Ansari and Ghoddousi (2018).

Distances to human settlements and roads have an important influence on the habitat suitability of bears in the DNP. Bears stay away from human communities until they are forced to look for food. During the winter, they may come upon abandoned grain and agricultural storage facilities. Furthermore, our results confirmed that the human impact, as for example the presence of roads, contributes negatively to the prediction of the bear’s habitat (Mateo-Sánchez et al., 2014; Farashi et al., 2016; Coogan et al., 2018; Almasieh et al., 2019). If on one hand the roads are increasing human-bear interactions, on the other hand, the road favours connections with other bear populations that are scattered or unevenly distributed (Kaczensky et al., 2003; Bischof et al., 2017). Finally, it is important to remember that the roads represent one of the most important causes of death for bears (Herrero et al., 2005; Mohammadi et al., 2021a). The roads used by shepherds can be a facility for the bear’s seasonal migration and foraging for food (pers. obs. by authors). A large amount of a bear’s diet consists of vegetables (Nawaz, 2008; Nawaz et al., 2019). The field-based observation showed that agricultural cultivation and pastoral farming are the primary land use practices also in the zones close to the DNP borders. Large herds of domestic cattle grazing freely inside the bears’ primary habitat create competition for resources. The Himalayan area is strongly susceptible to the effects of human impacts, mainly on land use, negatively influencing the bear presence and distribution (Dar et al., 2021). Throughout their historical range, human land use practices have played a significant role in the bear decline (Zedrosser et al., 2011; Wolf and Ripple, 2017), negatively impacting bear habitat (Dixo et al., 2009), survival (Cushman et al., 2018; Kaszta et al., 2019), dispersal pattern (Cushman et al., 2014; Mateo-Sánchez et al., 2014), genetic diversity (Kaszta et al., 2020), and ecological connections in habitat patches (Scolozzi and Geneletti, 2012; Gao et al., 2013).

The bears migrate following the change in seasons, spending the summers on plateaus at higher elevations and the winters in the valleys, where they hibernate. Cattle grazing and fuelwood collection are two main human-caused disturbances of bear habitat that commonly occur in the valleys. Our findings revealed that DNP covered 73% of the bears’ habitat regions. However, other tiny, well-located, and suitable places were discovered along the DNP’s borders; these areas should be conserved since they exhibit a high degree of environmental similarities to those found in the park. On the one hand, we discovered a robust network between the bear-rich areas of the DNP, but on the other hand, because of anthropogenic activity, the areas outside the DNP were not strongly connected to one another (e.g., Fida et al., 2024). This situation means that bears in DNP may survive in isolated populations with minimal genetic diversity (Nawaz, 2008). Several studies have shown that maintaining core habitats and promoting migration across appropriate habitat patches are necessary to ensure the appropriate survival of vertebrate populations (Margules and Pressey, 2000; Mech and Hallett, 2001). In order to resolve human-bear conflicts, spatial analyses, such as landscape connectivity, are crucial to improving effective management plans to protect bear population and distribution (Lecis and Norris, 2003; Goursi et al., 2021; Mohammadi et al., 2021b; Rehan et al., 2024).

With a predicted loss of 56,501 km² (17.52%) by 2050, mostly due to human activity, Pakistan’s bear habitat is predicted to be seriously threatened (Su et al., 2018). A research-based conservation plan is required to monitor the bear population and the effects of the factors affecting its survival. One such approach would be to preserve the landscape corridors (Bellemain et al., 2007). Although there is a healthy bear population in the park, there may be serious risks related to ongoing construction projects, infrastructure, and human access to bears’ primary habitat. In order to support the unchecked growth of tourism in PAs, it should be banned to construct infrastructure in the DNP and its neighbouring areas that could disrupt bears’ natural habitat. The road that passes through DNP separates the bear’s habitat, and when summer traffic increases, the road becomes busier, causing habitat fragmentation and population isolation. We should establish safe movement corridors, underpasses, and tunnels to allow bears and other wildlife to freely roam across the appropriate ecosystems found along the major road. Our findings imply that the DNP limits should be expanded to include areas where appropriate habitat and pathways were discovered. The Pakistani government recently declared two locations within the bear range to be protected areas: Nanga Parbat National Park and Himalaya National Park. Preservation and restoration of species movement corridors should be the primary focus of conservation efforts to preserve genetic diversity and gene flow among bear populations (Straka et al., 2012).