Potential Distribution of the Critically Endangered Chinese Pangolin (Manis pentadactyla) in Different Land Covers of Nepal: Implications for Conservation

Abstract

Anthropogenic activities have driven many wildlife species towards extinction. Among these species, the geographic distributions of many are poorly documented, which can limit the effectiveness of conservation. The critically endangered Chinese pangolin (Manis pentadactyla) is experiencing population decline throughout its range due to land-use changes and illegal trade for food and traditional medicine. Using distribution records and maximum entropy modeling, we predicted the potential distribution of the Chinese pangolin across Nepal. Most suitable potential habitats of the Chinese pangolin occurred in forest areas of the mid-hill region in central and eastern Nepal, followed by cultivated land. Almost all potential suitable habitats of the Chinese pangolin occurred outside of protected areas, and most of them were encroached upon by cultivated land, human settlements, and infrastructure developments. The results from this study provide baseline information on the potential suitable habitats of the Chinese pangolin in Nepal, which helps to develop site- and species-specific management plans and to identify priority areas to minimize the current threats to the pangolin and enhance the stewardship of species conservation.

1. Introduction

Many wildlife species are at risk of extinction in response to anthropogenic activities, including climate change, habitat alteration or loss, biological invasion, infrastructure, or combinations of these and other factors [1,2,3]. Consequently, vertebrate extinction risk is increasing, and for mammalian species, extinction rates have accelerated during the last several decades [4,5]. However, the extinction rates of species and the relative effects of factors causing extinction risk vary across regions, particularly in relation to changes in land cover [6]. Conservation of species depends on understanding of their ecological needs, which first requires knowledge of species’ distributions and habitat associations [7,8]. Identifying species’ potential distribution is a prerequisite for species conservation [9,10] and can play a vital role in ecological restoration by, for example, connecting meta-populations across landscapes and population viability analyses [11,12,13,14].

Of the eight pangolin species worldwide, Nepal supports the Chinese pangolin (Manis pentadactyla) and Indian pangolin (M. crassicaudata). In addition to Nepal, Chinese pangolins occur in the Asian countries of Bangladesh, Bhutan, China, India, Lao PDR, Myanmar, Thailand, and Vietnam [15,16]. The Chinese pangolin is listed as critically endangered on the IUCN Red List of Threatened Species and is included in Appendix I of the Convention on International Trade in Endangered Species of Wild Fauna and Flora and in the protected mammal species of Nepal’s Wildlife Protection Act 1973 [17,18]. The primary cause for their threatened status is poaching for illegal trade of scales and meat [16,17,18,19,20,21,22,23,24,25,26].

Chinese pangolins reportedly prefer forests over other land types [24,27,28,29]; however, the extent of some of the forest area in Nepal is declining due to conversion to cultivated land and settlements, dams and roads [30], timber harvest, human-caused forest fires, and livestock grazing [31]. The resulting habitat fragmentation can also isolate Chinese pangolins to small patches [32] and potentially alter their behavior [33,34].

There are limited data on the population and habitat of the Chinese pangolin in Nepal [31]. Generally, concerned entities have prioritized conservation in areas where species occurrence has been confirmed. Studies are very localized in a few locations, and the effects of current land cover are under-studied. The Pangolin Conservation Action Plan for Nepal (2018–2022) also focuses more study on the habitat selection of Chinese pangolin [31]. Here, we used species distribution modeling (SDM) to study the potential distribution of the Chinese pangolin in existing land cover types and to provide the baseline information for developing an area- and species-specific conservation strategy in Nepal. Maximum entropy modeling (Maxent; [35]) is widely used to estimate species distribution from presence-only data and is particularly well suited for modeling threatened species with limited available data [36,37,38]. In addition, maximum entropy models provide flexible choices of thresholds to use with binary outputs [35] in order to quantify and define reserve areas for threatened species such as the Chinese pangolin. This data will be beneficial for establishing corridors and/or reserve designs for the critically endangered Chinese pangolin, like other species whose populations are declining [39]. The aim of this study was to identify the potential habitats of the Chinese pangolin within existing land covers throughout Nepal in order to guide the prioritization of future research, monitoring, and conservation.

2. Material and Methods

2.1. Study Area

Nepal comprises 147,181 km² (80°04′ to 88°12′E and 26°22′ to 30°27′N) and is bordered by China and India. Nepal is divided into five physiographic regions of increasing elevation: Tarai, Siwalik, mid-Hill, mid-Mountain, and High Mountain (Himal), with an overall elevational range of 60 to 8848 m [40]. These regions have correspondingly different climates, with that of the Tarai considered tropical, that of the Siwalik considered subtropical, those of the mid-Hills and mid-Mountain considered temperate, and that of the High Mountain consisting of subalpine and alpine climates. At present, Nepal has 77 districts within seven provinces. The official names of all of these provinces are not assigned, except for Bagmati, Gandaki, Karnali, and Sudur Pashchim. We collected Chinese pangolin presence data from 14 districts of eastern to western Nepal, including the Taplejung, Ilam, Dhankuta, Shankhuwasaba, Solukhumbu, Sinhupalchok, Sindhuli, Kavrepalanchok, Ramechhap, Dolakha, Bhaktapur, Kathmandu, Makwanpur, and Gorkha districts (Figure 1).

2.2. Data Collection and Management

We collected presence-only data on Chinese pangolin using transects, plots, and opportunistic observations. We established 137 transects (Figure 2) whose lengths varied from 500 to 1000 m (Mean $\pm \mathrm{SD}=554.7 \pm 141.93 \mathrm{m}$). The transects were in forests, shrublands, and in cultivated lands, and the length of each transect was determined by the physiography of the locality.

Signs of the Chinese pangolin’s occurrence were recorded up to 50 m on either side of the transect. In addition, we found 956 presence-only data from opportunistic surveys where the establishment of transects was not possible. All data were collected during 2008–2019 from eastern to western Nepal in different physiographic regions, from the Siwalik to Mid-mountain regions (Figure 1). Overall, we obtained 3066 occurrence records at elevations ranging from 378 to 2406 m. The Chinese pangolin’s home range size is <1 km²; therefore, we removed more occurrence points of this species through spatial filtering. After spatial filtering for 1 km resolution, only 110 presence locations were available for further analysis (Figure 1), and finally, the distance between two occurrences points was >1 km. We confirmed that these data were of Chinese pangolins because Indian pangolins are known to occur only in a small area of southern Nepal at elevations below 500 m, whereas Chinese pangolins are more widely distributed and reportedly occur below 2000 m [18,41]. To avoid confusion between Indian and Chinese pangolins, we used only direct sightings of Chinese pangolins in areas below 500 m elevation and, whenever possible, confirmed their occurrences through direct communication with key informants. In areas with elevations greater than 500 m, we used presence locations of Chinese pangolins based on direct sightings and indirect observations such as burrows, footprints, and scratches see [29].

2.3. Spatial Modeling

We prepared potential spatial layers of Chinese pangolins by initially selecting 19 bioclimatic variables at an arc resolution of 30 (ca. 1 km²) from Worldclim (www.worldclim.org; [42]) as the environmental layers. From these, we selected seven environmental variables, including Mean Diurnal Range (Bio_2), Iso-thermality (Bio_3), Temperature Seasonality (Bio_4), Mean Temperature of the Warmest Quarter (Bio_10), Precipitation of the Driest Month (Bio_14), Precipitation of the Warmest Quarter (Bio_18), and Precipitation of the Coldest Quarter (Bio_19) (Table S1), based on their statistical importance determined by jackknife analysis (Figure S1) for Chinese pangolin presence data and their correlation values (|r| < 0.70; Figure S2) to avoid model overfitting [43]. We used these environmental layers to estimate the potential distribution of the Chinese pangolin with maximum entropy modeling software (Maxent version 3.3.3 K; [37]). We constructed our Chinese pangolin distribution model with presence-only data and bio-climatic variables.

Maxent uses presence-only data to determine the potential area by using a set of variables using the presence of model species that are satisfied by environmental variables [44]. The software produces a map of habitat suitability from lowest to highest with values ranging from 0 to 1, respectively. Maximum entropy models are commonly used for SDM due to their continuous probabilistic output, threshold flexibility, and applicability for small sample sizes [35,38].

We ran Maxent using 75% of the presence data for calibration and the remaining 25% for model validation [45] with the bootstrapping procedure with 20 replications and a maximum of 500 iterations. We evaluated model validation and accuracy by calculating the Area Under the Curve (AUC) using Receiver Operating Characteristics (ROC) analysis [46]. However, because of criticism of the use of AUC for model evaluation [47], we also evaluated model performance using the True Skill Statistic (TSS) package: \tab TSS\cr; [48]. The TSS values range from -1 to +1, with the greatest value indicating good model performance with perfect discrimination while zero or less indicated no better performance than random. We calculated all of these statistical values using the program R Version 3.3.3 [49]. We converted Maxent probability scores from the Maxent output into the suitable and unsuitable habitats using the Maximum Training Sensitivity and Specificity Logistic Threshold using the raster calculator in QGIS. Maximum Training Sensitivity and Specificity Logistic Threshold is a prominent method for presence-only data because it produces the highest sensitivity values [50]. We used the suitable habitat raster layer to identify the extent of the potential distribution of Chinese pangolins occurring within existing land covers.

2.4. Land Cover Classification

To extract the land cover of this area, we used terrain-corrected (Level 1) 30 m resolution Landsat 8 Operational Land Imager (OLI) images obtained in 2016 to evaluate major land cover (U.S. Geological Survey; https://earthexplorer.usgs.gov; Table S2). We further verified all images, and a Flash Line-of-Sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) has been used for image processing. For the extraction of land cover, we used Support Vector Machine (SVM) [51,52,53] algorithms. Basically, SVM provides a flexible supervised classifier option with higher accuracy than that of the maximum likelihood (ML) classifier [54,55]. We assigned a Radial Basic Function (RBF) with maximum penalty parameters with a sum of 100. We described the classification equations of each kernel based on the section of ENVI version 5.3 (www.exelisvis.com). We created seven land cover classes, including forest, shrubs, cultivated, settlement, barren, grassland, and water bodies such as dams and irrigation canals (Table S3), using the classification scheme developed by Anderson et al. [56].

We evaluated land cover classification accuracy with a confusion matrix [57]. Overall accuracy (OA), user’s accuracy (UA), and producer’s accuracy (PA) were evaluated using field survey information collected during 2015, 2016, 2017, and 2018, as well as historical records including topographical maps developed by the Survey Department of Nepal [58] with scales of 1:25000 and 1:50000, high resolution Google Earth images, and government records and community recall. To collect ground truth data, we used the Global Positioning System (GPS). Using GPS, we validated the land cover against 2450 randomly selected reference points (at least 350) for each class.

3. Results

Model performance was high, based on the area under the ROC curve (AUC), true skill statistics (TSS), sensitivity, and specificity values of 0.95, 0.81, 0.96, and 0.84, respectively. Almost all (97.25%) of the presence locations fell under the predicted climatically suitable habitats shown by the Maxent results (Figure S3). The quantitative relationship of the ecological niche of the Chinese pangolin between predictive variables and the logistic probability of presence (habitat suitability) ranged from 360–534 for Temperature Seasonality (Bio_4), 1 to 29 $℃$ for Mean Temperature of the Warmest Quarter (Bio_10), and 377–1114 mm for Precipitation of the Warmest Quarter (Bio_18) (Figure S4). Climatic variables Bio_4, Bio_10, and Bio-18 had the greatest contributions to model development (Figure S5).

The model predicted about 28,768 km² of suitable habitat for the potential occurrences of Chinese pangolins in Nepal (Figure 3). However, of the potential distribution of Chinese pangolins, 12,713 km² (44%) was in modified land cover, mainly in cultivated land, human settlements/infrastructure developments, roads, landslides, and irrigation canals and dam constructions (Figure 4,

Table 1. Chinese pangolin’s (Manis pentadactyla) predicted suitable bioclimatic range in different land covers according to Nepal’s physiographic division; the suitable habitat is encroached upon by settlements, cultivated land, barren land, and water bodies in Nepal.

Land Cover Type	Chinese Pangolin’s Potential Suitable Habitat (Km2)					Total Area (Km2)
	Tarai < 300 m	Siwalik 300–1000 m	Mid-hill 1000–3000 m	Mid-mountain 3000–5000 m	High mountain > 5000 m
Settlements	17.54	10.68	245.60	14.01	0.00	287.83
Cultivated	260.64	578.83	9868.28	1571.90	0.04	12279.69
Forest	189.69	2683.96	8841.47	2223.88	6.57	13945.56
Shrub	1.88	93.97	877.04	255.10	0.37	1228.36
Barren	24.27	143.43	179.67	15.22	0.13	362.72
Grassland	0.13	25.24	323.27	168.10	2.22	518.96
Others	7.71	21.91	86.65	28.64	0.31	145.22
Total	501.87	3558.01	20421.98	4276.84	9.64	28768.35

). Of the seven land cover classes, forest cover (48%) and cultivated land (43%) comprised the greatest amounts of potential distribution, followed by shrubs (4%), grasslands (2%), barren land (1%), settlements (1.0%), and others, including water bodies (< 1%) (Figure 4, Table 1).

Among areas of potential distribution, the greatest area (71%) was in the mid-Hill region between 1000 to 3000 m of elevation from eastern to central Nepal (Table 1, Figure 5), followed by the mid-Mountain region (15%) with 3000 to 5000 m of elevation, and 14% in the Tarai and Siwalik regions at below 1000 m of elevation (Figure 5, Table 1). Most of the potential distribution (94%) was predicted outside of existing protected areas (Figure 6).

The Bagmati Province supported the largest area of potential distribution (41%), followed by Province 1 (40%), Gandaki Province (18%), and Province 2 (2%). Comparatively, more potential distribution occurred in forested areas (19%) of total suitable areas in both Province 1 and Bagmati (Figure 4,

Table 2. Chinese pangolin’s (Manis pentadactyla) predicted suitable bioclimatic range in different land use in Nepal according to province. Suitable habitats are encroached upon by settlements, cultivated land, barren land, and water bodies in Nepal.

Land Cover Type	Chinese Pangolin’s Suitable Habitat (Km2)				Total Area (km2)
	Province 1	Province 2	Bagmati Province	Gandaki Province
Settlement	35.45	0.26	187.68	64.44	287.83
Cultivated	5452.41	29.61	4576.03	2223.10	12281.16
Forest	5358.27	527.41	5506.59	2551.55	13943.83
Bush	190.04	55.46	815.11	167.95	1228.57
Barren	120.12	27.33	192.32	22.91	362.68
Grassland	69.24	2.48	341.76	105.54	519.02
Others	53.70	3.55	72.23	15.78	145.25
Total	11279.23	646.11	11691.72	5151.29	28768.35

), followed by 9% in the Gandaki Province. Among the total suitable areas that are in cultivated land, Province 1 supported the most (19%), followed by the Bagmati Province (16%) and Gandaki Province (8%) (Figure 4 and Figure 5, Table 2). The remaining land covers comprised < 1% of the remaining potential distribution.

4. Discussion

Our Chinese pangolin distribution model performed well and indicated the potential distribution throughout Nepal. We acknowledge that our distribution model was based on the available locations concentrated in the mid-hills of Nepal. It is difficult to obtain occurrence records for low-populated endangered species such as the Chinese pangolin. The knowledge of the ecology and distribution of the Chinese pangolin in Nepal is limited. However, research on captive breeding of the Chinese pangolin has revealed that temperatures between 18–27 $℃$ are suitable for their survival [59]. The low variation in the Temperature Seasonality of the predicted area of Chinese pangolin habitation in the mid-hill regions of Nepal may support a potential habitat for this species. However, the Chinese pangolin does have a limited ability to regulate body temperature, and they also have a low metabolic rate [60]. The extreme temperature in some months may affect the activities of Chinese pangolins, although they can adjust by staying inside the burrow, where they can dissipate heat through the soil and the burrow entrance [61,62,63]. The physiography in the mountain region is steeper; Chinese pangolins prefer to construct burrows between 30$°$–60$°$ steep slopes [64], which may help to maintain a stable diurnal temperature inside the burrows. Under such circumstances, even when air temperatures outside of the burrows are extreme, Chinese pangolins are potentially able to perform normal activities because of the more stable diurnal and seasonal temperatures inside burrows [65].

The forest areas are reported to be a resource for termites and ants [66,67,68,69], common prey species of Chinese pangolins. We found a greater potential distribution of Chinese pangolins in forest lands. This has been corroborated by the findings of other studies in Nepal [27,28,29]. Nepal is rich in wildlife diversity as a consequence of community-based forest management systems [70,71]. The extent of forests and shrubs in the hills and mountainous regions of Nepal is increasing [72], therefore resulting in more extensive refuge areas for both Chinese pangolins and their prey species [73].

The largest area of potential habitats for Chinese pangolins was found in cultivated areas. These areas are rich in ant and termite populations [74]. This is because, in practice, the cultivated lands in Nepal are bordered by trees and shrubs which can also support occurrences of the Chinese pangolin [75,76,77]. However, these areas have been heavily encroached, resulting in habitat fragmentation and the reduced occurrence of the Chinese pangolin [29]. The frequent disturbances by farmers in cultivated land limit burrows for shelter and movement of individuals (HBK, personal observation). Thus, the forests near cultivation have often been observed to be found as highly suitable habitats for Chinese pangolins [29]. Grass- and shrublands, which are important grazing lands for livestock, have a lower suitability than forested lands as Chinese pangolin habitats [78]. This could be due to the low prey availability, high disturbance pressure, and the negative effect of livestock grazing on Chinese pangolins [29].

Though the Chinese pangolin’s greatest potential distribution occurs in hilly and mountainous regions, these areas are currently fragmented by human infrastructure development [79]. The most common infrastructures, including settlements and unplanned road construction, are threatening the mammalian species in Nepal, including Chinese pangolins [29,32,80]. These actions may restrict the distribution of Chinese pangolins. For example, a study in Eastern China documented that 50% of the potential habitat for this species was lost between 1970 and 2000 [81]. Moreover, the cover of potential distribution of Chinese pangolins lies outside the protected areas of Nepal. This may be one of the major threats for the poaching of Chinese pangolins in Nepal. The conservation of the species may also be further curtailed by Nepal’s species conservation strategies/programs, which are mainly focused on protected area management activities. People outside the protected areas use Chinese pangolins meat and scales for their medicinal value (S. Tamang, Sindhupalchok, personal communication; [23,24,82]). Therefore, governmental consideration of additional protected areas, as well as conservation easements to protect Chinese pangolin habitats, is warranted.

Some action towards the conservation of Chinese pangolins is being undertaken. The government has launched the Pangolin Conservation Action Plan for Nepal (2018–2022) [31]. This Action Plan is designed to initiate immediate actions to identify potential habitats and the regular monitoring of populations and suitable habitats, in part to mitigate poaching and illegal trade. Our modeling provides the site-specific baseline information that can be used to prioritize areas for implementation of this Action Plan. In addition to the actions identified, we recommend that the government prioritizes the community-based conservation strategies for the Chinese pangolin. Community-based action, such as active participation of local people in government policy making, resource use, and involvement in awareness programs, would increase the effectiveness of species conservation [83,84,85]. In addition, the forests adjoining cultivated areas in the mid-hills are important for the Chinese pangolin; we recommend limiting anthropogenic disturbances in these areas for the betterment of Chinese pangolin population.