The populations of many species of large mammals occur in small isolated and fragmented habitat patches in the human-dominated landscape. Maintenance of habitat connectivity in fragmented landscapes is important for maintaining a healthy population of large mammals. This study evaluated the landscape patches and their linkages on two carnivores (leopard and Himalayan black bear) and seven prey species (northern red muntjac, chital, sambar, wild pig, Himalayan goral, rhesus macaque, and langur) between Chitwan National Park (CNP) and Annapurna Conservation Area (ACA) by using the least-cost path approach and the Linkage Mapper tool in ArcGIS. A total of 15 habitat patches (average area 26.67 ± 12.70 km²) were identified that had more than 50% of the total studied mammals. A weak relation among the habitat patches was found for chital and sambar (Cost-weighted distance CWD: Euclidean distance EucD >100), showed poor connectivity between the habitat patches, while the ratio of CWD and EucD was low (i.e., low least-cost path) between the majority of the patches for muntjac, wild pig and leopard hence had potential functional connectivity along the landscape. Similarly, a low least cost path between the habitat patches located in the mid-hills was observed for Himalayan goral and Himalayan black bears. Furthermore, the multi-species connectivity analysis identified the potential structural connectivity between the isolated populations and habitat patches. Therefore, these sites need to be considered connectivity hotspots and be prioritized for the conservation of large mammals in the landscape.

Occurrence points collection

The study blocks were divided into four distinct blocks, namely A, B, C, and D considering landscape features, the main courses of rivers, and topographical features (Figure 2C). Block A covers the BCF, part of CNP, and surrounding areas of BCF (Kabilas, Jugedi, Kerabari, Chaukidanda, Simaldhap) up to the Mahabharat range (it runs closely parallel to the Chure range and separates the Terai with the Hill region, i.e., mid-hill) of Chitwan district. Block B covers human-dominated mid-hill landscapes such as Devghat, Bandipur, Abu Khairani Rural Municipalities, and Vyas Municipality of Tanahun district. It follows the Seti and Trishuli River basins along with mid-hills. Block C covers the Bhimad Municipality, parts of Rishing Rural Municipality, Ghiring Rural Municipality, Magde Rural Municipality and Shuklagandaki Municipality of Tanahun District, and Rupa Rural Municipality of Kaski District along the Seti River basin. Block D covers Bharatpokhari, Nirmalpokhari, Pumdibhumdi, Panchase, Lumle, Ghandruk, Landruk, Deurali, and the Australian Camp area. This block harbors four types of forests: national forest, community forest, protected forest (Panchase), and conservation area (Annapurna).

Transects were laid for the collection of presence points of selected species (nine species of mammals) in the landscape (Figure 1). The presence points of ungulates and primates were collected on the basis of direct sighting whereas, the presence points of the carnivores were collected on the basis of signs left by them (e.g., scats, scrapes, pugmarks, scent spray, etc.). Transect size and length were determined based on forest patch size. After identifying forest patches using a topo-base map (Esri, 2017), transects were overlaid, with patches selected based on diameter; patches less than 2 km in diameter were excluded. Transects (150 out of 164) were systematically laid out according to patch size and accessibility in four blocks (31 in A, 35 in B, 38 in C, and 46 in D). Inaccessible areas (14 transects) due to deep river gorges, steep mountains, and swampy lands were excluded. Transect lengths ranged from 1.18 to 7.84 km, with a minimum 500 m separation in regular forest patches, varying in scattered habitats like Mid Hills (Figure 2C, Table S1). We also collected the presence of those mammals opportunistically from other possible sites of the study area (e.g., croplands, river banks). These presence coordinates were recorded by using the Global Positioning System (GPS- Garmin eTrex 10). The collected occurrence data were spatially filtered in 30 m by using the Spatially Rarify Occurrence Data tools of SDMtoolbox 2.0.0 in ArcGIS (Brown, 2020; Kaboodvandpour et al., 2021). The filtered data were converted into .CSV format for Maxent modelling (Table 1). The large mammals whose presence locations were less than 25, were removed from further analysis.

Environmental variables

To minimize the risk of over-fitting the model and develop the most parsimonious model, the environmental variables were selected based on field knowledge, experts’ suggestions, and an extensive literature review of studied large mammals (Dickman & Marker, 2005; Mishra, 1982; Rather et al., 2020; Watts et al., 2019). The slope and terrain ruggedness index (TRI) were extracted by using the digital elevation model (DEM) in ArcGIS 10.8 (ESRI, 2019). The classified image from Landsat (acquisition date 2020-03-17) (Landsat 8, Operational Land Imager (OLI)) was used for calculating the Euclidian distances to the nearest forest, grassland, water sources, developed area or human settlements, and cropland. We classified the images into eight different classes (Water sources, barren areas, grassland, riverine forest, Sal-dominated forest, mixed forest, cropland, and developed area) by using supervised classification based on the ground-truthing points (Adhikari et al., 2022). Among these classified eight classes, we merged riverine forest, Sal-dominated forest, and mixed forest as a single forest layer. We extracted water sources, grassland, forest, cropland, and developed areas from the available data and calculated Euclidian distances in ArcGIS 10.8 to be used as environmental variables for modeling. The Normalized Difference Vegetation Index (NDVI) is the most popular and used to quantify the greenness of the vegetation, and vegetation density and detect the changes in plant health using red and near infra-red bands of a remotely sensed image (Pettorelli et al., 2011; USGS, 2022; Yengoh et al., 2015), hence we selected NDVI as one environmental layer for mammals. Additionally, the modified Normalized Difference Water Index (MNDWI) is calculated by using the green and Short-wave Infrared (SWIR) bands and it enhances the features of open water. MNDWI also minimizes the features of developed areas that are correlated with open water in other indices (Xu & Guo, 2014; Xu, 2006). Furthermore, the Normalized Difference Built-up Index (NDBI) is a ratio that minimizes the effects of terrain brightness differences and atmospheric effects (Zha et al., 2003). Two spectral bands NIR and SWIR are used to enhance the build-up or developed area, thus differentiating built-up over the natural area. The values of each environmental variable were extracted at presence locations (Table 1). For the layer of prey richness of leopard, the suitability map of preys was calibrated as 0 for absent and 1 for the present of the species based on mean equal test sensitivity and specificity logistic threshold. Then, these layers were   combined as a single layer.

A total of 13 environmental variables were used for the modelling (Table 2). The variables were differed on the basis of nature of the mammals (Table 2). The selected variable layers were converted into ASCII format with the same resolution, extent and projection system. The spatial resolution of 30 m and UTM 45 N projected coordinate system was used for the modelling.

Habitat suitability models

Maxent develops a model based on a series of features (environmental variables) (Phillips et al., 2006). Two types of data (occurrence data and environmental layers) were used for processing in the Maxent program (Phillips et al., 2006). The CSV file of the occurrence points in the samples menu and all selected variables layers in ASCII format in the environmental layers’ menu bar were loaded for analysis. The replicates and replicated run type were fixed 25 and subsample respectively. The Maxent model ran with 25 iterations and 1000 background points with 70 % of the points used as training data and 30 % points used as validation of the model. The output of the model was logistic. The performance of the model was evaluated based on AUC values of the receiver operator characteristic (ROC) plot analysis (Phillips, 2008; Phillips et al., 2006; Phillips & Dudík, 2008). The value of the predicted suitability ranges from 0 to 1. The logistic probability of suitability was further regrouped as 0 – 0.2 = unsuitable, 0.2 – 0.4 = moderately suitable, 0.4 – 0.6 = suitable, and 0.6 – 1 = highly suitable (Ansari & Ghoddousi, 2018; Kogo et al., 2019). All the spatial analysis and classification were performed in ArcGIS 10.8 (ESRI, 2019). We used these results of habitat suitability to identify the habitat patches of the species and preparation of resistance layer. 

Landscape resistance 

The resistance or cost map was prepared using a raster habitat suitability map (Figure S1). Every cell on the map has a numeric value that indicates the cost that should be paid to pass through each cell (Bagli et al., 2011; Morovati et al., 2020). The cost map was developed by inverting the value of habitat suitability using the following formula (Almasieh et al., 2019; Morovati et al., 2020).

Cost = 100 × (1 - habitat suitability)                                                               

The lower cost is assigned to highly suitable areas whereas the highest cost is for the habitats with low suitability (Almasieh et al., 2019; Morovati et al., 2020).

Identification of habitat patch

The continuous probability of occurrence was converted to binary predictions of presence and absence based on average equal sensitivity and specificity threshold. The predicted maps of all species were combined to identify the species richness of an area. The habitat patches were defined based on the number of species predicted in that area. About 50% of species’ present areas with 5000-pixel size were defined as the patch (Sahraoui et al., 2017).

Modelling connectivity

The least-cost path (LCP) algorithms were used to identify a path or corridor (or linkage) between two geographical locations (Adriaensen et al., 2003; Unnithan Kumar & Cushman, 2022). The program Linkage Mapper 2.0.0. (McRae & Kavanagh, 2011) was used to identify the LCP for the movement of mammals from one patch to another. The Linkage Mapper identifies the closer patch, develops the networks between the patches, and calculates the least-cost distance and paths (McRae & Kavanagh, 2011). The lower cost-weighted distance is regarded as the strong corridor between two patches. The lower value of the least cost path (LCP) is regarded as the lower resistance for the movement of the animals i.e., it is highly suitable (Unnithan Kumar & Cushman, 2022). Two metrics were calculated to show the quality of each linkage. One is the ratio of cost-weighted distance (CWD) and Euclidean distance (EucD) that separate each pair of habitat patches. If the ratio of CWD and EucD is equivalent to 1, it is regarded as the highest possible quality linkage (Dutta et al., 2016). The second metric is the ratio of cost-weighted distance (CWD) and the length of the least-cost path (LCP). This provides the average resistance encountered along the optimal path between the habitat patches. The least-cost path of each species was identified and then, combined to find the single multi-species corridor between the patches using the raster calculator tools of ArcGIS. The Kernel density estimation method was used to identify the hotspots (Thakali et al., 2015) for the connection of isolated populations of mammals in the patches.