Potential risk zone for anthropogenic mortality of carnivores in Gandaki Province, Nepal

Abstract Anthropogenic pressures in human‐dominated landscapes often contribute to wildlife mortality. Carnivores are especially vulnerable to human‐induced mortality due to the perceived threat to livestock and humans. Despite having widespread conservation implications, carnivore mortality data have been largely underutilized within Nepal. This study utilized Maxent to identify high‐risk areas and explore the contribution of habitat attributes associated with carnivore mortality using the casualty database within the Gandaki province of central Nepal. We categorized the risk to carnivore species in three taxonomic groups, Felid, Viverridae, and Herpestidae, and identified a 3704‐km2 area within the province at high risk for carnivore casualty. The middle mountains were the riskiest physiographic zone, and the Annapurna Conservation Area represented the largest risk zone among the four protected areas. Agricultural land was the most problematic area in terms of carnivore casualty. The human population was positively associated with high‐risk areas and the number of casualties, whereas protected area cover had a negative association. This study identified that the common leopard was at the highest risk of mortality and therefore would benefit from the implementation of an action plan and species‐specific conservation strategies, especially within identified high‐risk zones. An expansion of protected areas in the middle mountain region would serve to greatly reduce carnivore casualty. Species distribution modeling can be further used with national‐level spatial and temporal mortality data to identify the most prominent casualty times and pinpoint potential casualty locations throughout the country.

Understanding the attributes of mortality and the severity of risk in specific areas can be vital to conservation endeavors, especially in areas of high HWC. For example, studies regarding the circumstances of mortality are an important consideration in understanding population dynamics and the risks associated with the loss of breeding pairs (Bekoff, 2001;Mateo-Tomas et al., 2012;Seiler & Helldin, 2006;Taylor-Brown et al., 2019). The long-term survival of adults, especially females, is crucial for maintaining a viable population of carnivores (Knight & Eberhardt, 1985;Purvis et al., 2000;Weaver et al., 1996).
Mortality data can be utilized to determine the sustainability of a given population and can provide insights on the methods needed to minimize overall mortality rates (Goodrich et al., 2008).
To better understand HWC and the spatiotemporal interactions between humans and wildlife, the species distribution model (SDM) has been employed to map areas of importance to carnivore survival (Kalle et al., 2013;Kramer-Schadt et al., 2013). The usage of SDM in the field of predicting risk for wildlife can help conservation endeavors by narrowing the focus to include only the economic and management efforts which are required to develop conservation plans for specific sites within a given timeframe (Mateo-Tomás et al., 2012). Maxent is a widely used approach for SDM (Fitzpatrick et al., 2013;Phillips et al., 2006) and utilizes sample background location data contrasted against presence location data to estimate potential occurrence zones (Phillips et al., 2006). It has demonstrated the best predictive power across all sample sizes and has shown superior performance compared with other modeling software even when sample sizes are low (Elith et al., 2011;Wisz et al., 2008). This approach has previously been used in various studies to map suitable habitats for wildlife under the influences of anthropogenic, environmental, and topographic variables (Bai et al., 2018;Panthi et al., 2019;Sharma et al., 2020;Thapa et al., 2018). Though Maxent has been used to assess and map risk areas as well as predict mortality risk zones for various wildlife around the world (Garrote et al., 2018;Mateo-Tomás et al., 2012), its application to carnivores in Nepal is novel.
In Nepal, most of the studies on the HWC have focused on risks to human lives, livestock, crops, and infrastructure (Acharya et al., 2016;Aryal et al., 2014;Koirala et al., 2012;Sapkota et al., 2014;Sharma et al., 2021). However, dimensions of risk and casualty to wildlife have been largely understudied. This study, therefore, aimed to map potential risk zones for carnivore mortality using SDM. We aimed to identify a potential risk zone for four carnivore groups and analyze the attributes associated with their mortality. We hypothesized that high-risk areas for these carnivores are associated with agricultural lands and anthropogenic structures (e.g., road, settlements) in non-protected regions. We also hypothesized that the risk area for the common leopard (Panthera pardus) is greater than the risk to smaller carnivores.
The Gandaki province includes four protected areas (Annapurna Conservation Area (CA), Manaslu CA, Dhorpatan Hunting Reserve (HR), and Chitwan National Park (NP)). The study area has forest coverage equaling 7138.3 km 2 (32.6%) which is home to more than 55 mammalian species . Some of the major carnivore species found in the province are felid species (common leopard (Panthera

| Data collection
The study area was divided into eleven blocks in accordance with administrative districts. The study blocks represented five physiographic zones and included five protected areas of the province, hereby capturing the heterogeneity of physiography and protected/ non-protected regions. The division forest offices of each district, wildlife rescue center, and the protected area offices were visited for secondary data collection regarding wildlife casualties between January 2019 and February 2021. Information on incident attributes such as the name of affected species, location of mortality (GPS coordinates), cause of mortality, age of individual, and season of mortality was recorded. The casualty location specified in the database of the divisional forest, protected areas (PA), and wildlife rescue center was verified through field visit and key informant interviews.

| Selection of variables for modeling
We collected information on topographic, vegetation, and anthropogenic variables (Table 2) from several sources to use in maximum entropy modeling. The environmental variables considered in our models have been previously utilized for identification of suitable habitat (Sharma et al., 2020) and risk zone mapping (Karki & Panthi, 2021).
Digital elevation model (DEM) having 30-m resolution was downloaded from the Web site of United States Geological Survey (https://earth explo rer.usgs.gov/). Slope and aspect were derived from the DEM using ArcGis software (ESRI, 2017). Shape files of water sources were downloaded from the Geofabrik Web site (https://www.geofa brik.de/data/shape files.html) and converted to distance raster file using ArcGIS (ESRI, 2017). Since the climatic variables were not available in 30-m resolution, elevation was taken as a proxy of temperature. Data associated with forest cover were downloaded from the Global Forest Change (GFC) Web site (Hansen et al., 2013). The Enhanced Vegetation Index (EVI) time-series data of 2018 and 2019 which were derived from images of Landsat 8 were analyzed with the help of the Google Earth Engine.
The anthropogenic variables were also used in our models. The shape files available on the Geofabrik Web site (https://www.geofa brik.de/data/shape files.html) were used for information on the location of paths and roads. The information on locations of settlements was derived from the Department of Survey, Nepal. ArcGIS was utilized to create distance raster files of paths, roads, and settlements (ESRI, 2017). Land-use and land cover (LULC) data were downloaded from the International Centre for Integrated Mountain Development Web site (ICIMOD; http://www.icimod.org) (Uddin et al., 2015).
Information on population density was obtained from humdata Web site (https://data.humda ta.org/) and processed in ArcGIS. All variables were pre-processed in ArcGIS (ESRI, 2017) to convert in ASCII format and final spatial resolution of 30 m.

| Data analysis
The descriptive mortality data regarding the family, cause, age, and season were analyzed for the significant differences and association between the variables using chi-squared test of homogeneity and goodness of fit. The study focused on maintaining accuracy on prediction of high-risk zones for various groups of carnivores. Efficiency of model prediction is improved by accounting uneven sampling bias (Kramer-Schadt et al., 2013), and hence, spatial filtering in oversampled region was deemed useful (Phillips TA B L E 2 Initial variables used in VIF test stepwise elimination process to be further used in Maxent modeling.  (Pearce & Ferrier, 2000). Even though AUC is widely used model evaluation parameter, it is severely criticized by researchers (Lobo et al., 2008). To overcome this limitation, it is recommended to use AUC in combination with methods such as AIC (Akaike information criterion) or TSS (True skill statistics). Using R software (R Core Team, 2020), we used threshold-dependent method, TSS (sensitivity + specificity −1) which places equal weight on model sensitivity and specificity, with its value ranging from −1 to 1 (Allouche et al., 2006).
The threshold of maximum sum of sensitivity and specificity is recommended to generate information on distribution using presence-only data (Liu et al., 2013). This threshold was used to convert the continuous probability map into a binary map consisting of high-and low-risk probabilities. Areas with probability values above the threshold were categorized as high-risk zone whereas the areas with lower probability value than the threshold were categorized as low-risk zone. The high-risk zones obtained for each group of carnivores were not mutually exclusive and could overlap with each other; therefore, the overall carnivore risk zone was obtained by creating union of all four high-risk zone layers of carnivore groups in ArcGIS. The high-risk shapefile was then intersected over land use/ land cover, district, protected areas, and physiographic zones within the province to obtain the intersected data.
We performed Canonical Correspondence Analysis (CCA) (TerBraak & Smilauer, 2002), which is a multivariate constrained ordination technique. This method extracts major gradients among combinations of explanatory variables in a dataset. We selected this method because we had categorical district wise data and the independent variables were consistent within the sample site.
CCA was utilized to determine the association of high-risk area (in km2) within each district and the number of casualty in each district, with several environmental and anthropogenic variables (forest cover, livestock holding, human population, protected area cover) within each district. These variables are frequently mentioned in literatures (Distefano, 2005;Manral et al., 2016;Peterson et al., 2010;Sangay & Vernes, 2008) influencing the intensity of humancarnivore conflict. The data regarding high-risk area (km 2 ) for each district were extracted from Maxent (union of high-risk areas of four groups of carnivores) using ArcGIS, whereas district wise data on other variables such as forest cover, livestock holding, and human population were extracted from opendatanepal Web site (https:// opend atane pal.com/).

| Statistics of mortality
Among the total studied cases (n = 232), number of felid mortality was highest (51% out of total mortality cases), followed by Viverridae

| Total risk zone for carnivores
We found that a total of 3704 km 2 (17% of total area of province) was potential high-risk zone for carnivore mortality in the Gandaki province of Nepal ( Figure 3). The Tanahun district demonstrated the highest risk followed by the Syangja, Gorkha, and Kaski districts, each representing over 10% of total risk area (Annex 2). Similarly, we found that the overall carnivore mortality was highest in the middle mountain range (85% of total risk area), whereas 229 km 2 (6% of total risk area) was found inside protected areas. Additionally, agricultural lands were riskiest for carnivores, followed by forested areas, encompassing 65% and 28% of total risk area, respectively (Annex 2).

| Risk zone for leopard and small felid
A total of 2322 km 2 (11% of total area of province) was identified as potentially high risk for leopard casualties (Figure 3) (Annex 1). Out of the total risk-prone areas for leopard inside the province, agricultural lands represented 54%, followed by forested land at 40%.
The riskiest district for leopard casualties was Tanahun district (22% of total high-risk area) followed by Syangja (17%), Kaski (15%), and Gorkha (15%) districts. The result of a jackknife test indicated that variables contributing most common leopard risk were elevation, distance to road, and distance to settlement (Figure 4). Similarly, a total of 660 km 2 (3% of total area of province) inside the Gandaki province was concluded to be under high risk for small felid casualties, with the highest risk areas falling within agricultural lands (85%) ( Figure 3) (Annex 1). The Tanahun and Kaski districts represented 26% and 20% out of total high-risk area, respectively. Forest cover, distance to road, and distance to settlement were found to highly influence risk models (Figure 4).

| Risk zones for Viverridae and Herpestidae
A total of 1589-km 2 area (7% of total area of province) was considered at high risk for Viverridae and 717 km 2 (3% of total area of province) area at high risk for Herpestidae ( Figure 3) (Annex 1). Again, agricultural areas represented the highest risk, 89% and 80% of total risk area for the two groups of carnivores, respectively. Tanahun (20%) and Nawalparasi (17%) established themselves the riskiest districts for Viverridae, whereas Tanahun (23%) and Kaski (23%) were found to be riskiest for Herpestidae. Land use/land cover and distance to building highly influenced Herpestidae risk zone, whereas Viverridae risk zone was influenced by distance to road and land use/land cover (Figure 4).

| Accuracy of models
The average accuracy of the four models was good with AUC ≥ 0.89 and TSS ≥ 0.51. The AUC, TSS, and threshold used to convert the continuous probability map to binary high-risk/low-risk map of each models are represented in Table 3.

| Factors influencing number of casualty and risk area (km 2 ) within districts
CCA indicated that number of casualty and high-risk area (km 2 ) within the districts were closely related ( Figure 5). Both were positively associated with the human population whereas they revealed negative association with protected area cover within the districts. Monte Carlo 999 permutation test indicated that the analysis was significant (p < .05). CA1 primary axis variables explained 86.26% of variance (Table 4).

| DISCUSS ION
Our America , and cause-specific mortality of world's terrestrial vertebrates (Hill et al., 2019).
A total of 3704 km 2 (17% out of the total area of province) was identified as high risk for carnivore casualties. Districts within middle mountain regions of Gandaki province encompassed the highest risk area (Figure 6). In contrast, the lower risk zones were mostly located in the High Mountain and high Himalayan regions. Higher  have been reports of conflict in the region Pahari et al., 2021), which could easily have turned into wildlife casualty incidents. Low risk for carnivores in the Himalayan region could be explained by the presence of protected area management systems, which promote conservation of wildlife among local people (Schutgens et al., 2019), incentive-based programs (Spiteri & Nepal, 2008), and compensation schemes in the event of livestock casualty, hereby preventing retaliatory killings. In contrast, the middle mountains encompassed the largest risk area, likely due to its large area and presence of densely populated cities with lower conservation incentives. The low-and moderate-risk districts within the study area have lower human densities compared with high-risk districts (UNFPA, 2017). High population densities are often associated with high anthropogenic pressure and hence pose a greater threat to local wildlife (Pietersen et al., 2014). The CCA performed in this study also represented positive association between human population and high-risk zone (km 2 ) for carnivore casualty (p < .05), supporting earlier findings. Anthropogenic mortality in mammals has been found to escalate with increasing human footprints and human-associated impacts on the landscape (Gubbi et al., 2014(Gubbi et al., , 2021Hill et al., 2020).
Results of our study also reinforce positive association between human population and number of carnivore casualty.
The importance of protected areas in preserving biodiversity and sustaining wildlife is well known (Acharya et al., 2016;Paudel & Heinen, 2015). Even though the protected areas cover 45% of the total area of the Gandaki province, only 6% of total risk area was included within protected areas. The CCA performed in this study also revealed negative association between protected area coverage and risk area (r = −0.35, p = < .05) for carnivores. This demonstrates the efficiency of protected area in reducing carnivore mortality.
Protected areas inside Nepal are disproportionately located in high mountains and high Himalayan region (Paudel & Heinen, 2015).
Subsequently, carnivores and other fauna of the middle mountains and lowlands are not adequately protected, with most forced to survive in or near human-dominated forested landscapes (Paudel et al., 2012). Therefore, the main risk area (>94%) fell within unprotected areas. This result supports earlier findings of Acharya et al. (2016), who concluded that human-dominated landscapes and the regions outside protected areas are the human-wildlife conflict hot spots in Nepal.
Larger mammals are likely to be at more risk due to increasing anthropogenic pressure Thapa, 2014). Our study explored common leopard to be at more risk than other small carnivores. Since the common leopard is a large carnivore, it is the only species in our study which is capable of causing significant loss of livestock and has the potential to cause human casualties (Acharya et al., 2016(Acharya et al., , 2017. The mid-hill districts represented the greatest risk areas for common leopard. Several incidents of HWC were recorded from these areas (Acharya et al., 2017;. The incidences of human and livestock casualties have acted as a catalyst toward increasing leopard fatalities due to retaliatory killings. Agricultural lands inside the province represented the most at-risk land-use type for common leopard followed by forest areas and grasslands. The increasing urbanization and fragmentation of landscapes has caused significant reduction in the prey base of big cat species such as leopard (Puri et al., 2020;Schneider, 2001 agricultural lands, which harbor livestock such as cattle, buffalo, goats, and poultry, are increasingly being used by common leopards as hunting grounds (Kabir et al., 2014;Kshettry et al., 2017). Naha et al. (2020) reported high occurrence of conflict when livestock were allowed to graze freely within multi-use areas. Forested areas were also established as quite risky for common leopard, due in part to the species' need for large home ranges. Leopard often regularly patrol their range in search of a mate and to maintain their territory (Stein et al., 2011).
Despite being elusive, this behavioral necessity of covering large home range makes the leopard more vulnerable to anthropogenic threats, hereby increasing casualty risk even within the forest and fringe areas (Acharya et al., 2017;Naha et al., 2020).
Distances to roads and settlements were highly important variables to develop risk zone of common leopard. This result concurs within other documented literature (Edgaonkar & Chellam, 2002;Jacobson et al., 2016;Kumbhojkar et al., 2020) which mentioned high numbers of HWC events in areas of high anthropogenic pressure. Many studies have reported roadkill as one of the most significant causes of wildlife mortality (Baskaran & Boominathan, 2010;Sayyed & Mahabal, 2015).
In contrast to leopards, small felids were the group at lowest risk within the province. This may be due to their elusive nature and diet of mice and shrews, which are of low importance to humans.
Small felids in the study (leopard cat and jungle cat) are sympatric species occurring in various ranges within the Indian subcontinent (Mukherjee et al., 2010). These two cats are considered being at similar risk, given their comparable size and overlapping prey preference (Majumder et al., 2011;Rajaratnam et al., 2007). Despite being at the lowest risk overall, this study predicted agricultural areas to represent the most risk for small felid mortality. Since agricultural areas harbor mice, rats, and shrew, they are presumed to attract small felids. In agricultural countries such as Nepal, most of the households in rural areas are built within close proximity of agricultural lands in order to make the fields and agricultural work more accessible and efficient. When these small felids wander through agricultural area in search for the prey species, they are bound to stumble upon humans, which may result in HWC (Inskip & Zimmermann, 2009).
Hence, the vulnerability of small felids in the agricultural area is understandable.
Our study found that members of the Viverridae family were also at higher risk in agricultural lands. Civets are mostly omnivorous, but feed primarily on fruits and berries (Khan et al., 2019). These animals range agricultural lands in search of food and are consequently F I G U R E 6 District wise carnivore risk zone categories where risk area greater than 500 km 2 represents high-risk district, risk area between 500 and 100 km 2 represents medium risk district, and risk area less than 100km 2 represents low-risk district at risk from humans. Land use and distances to roads were recognized as the most important variables for prediction of Viverridae risk zones. Viverridae are nocturnal species, which forage at night (Bu et al., 2016), and as such, they are found upon road ways after dark, resulting in incidents of roadkill. Species within Herpestidae are more likely to occur in grassland and forestlands inside complex burrow systems (Mahmood & Nadeem, 2011). Distances to building and land-use land cover were prominent risk variables for their mortality, indicating areas with high anthropogenic pressure. However, casualties of mongoose species may also be associated with their reputation as invasive and pest species (Ćirović et al., 2011;Hays & Conant, 2007

ACK N OWLED G EM ENTS
Ministry of forest, environment, and soil conservation, Gandaki province, Nepal, funded this research. We would like to thank Pokhara zoological park and wildlife rescue center, Institute of forestry, Gandaki province forest directorate, Division forest offices, Annapurna conservation area office, Manaslu conservation area office, Dhorpatan hunting reserve office, and Chitwan national park/ buffer zone office for their support.

CO N FLI C T S O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
The dataset that is associated with this study is available at dryad

A N N E X 2
Total high-risk area (%) for each taxonomic group in various sub-layers. The high-risk area for each group of carnivore is not mutually exclusive and can overlap with each other. Each category represents % of total high-risk area within the province whereas the protected area category represents % of total high-risk area within the protected area in province