Modeling population extirpation rates of white-bellied and giant pangolins in Benin using validated local ecological knowledge

Pangolins are globally threatened by unsustainable hunting for local use and illegal international trade, plus habitat loss. In Benin (West Africa), white-bellied and giant pangolins have experienced a contraction in their distribution areas and population decline during the last two decades. To better understand the factors underlying declines in these species, we investigated extirpation rates of populations over the last 20 years. Because pangolins are elusive species difficult to monitor by standard methods, the status of populations has been assessed through a local ecological knowledge (LEK) approach. We collected information on persistence or extirpation status of pangolins from 156 localities. A binomial model was built to predict population persistence probability as a function of past and ongoing landscape changes, initial abundance (1998), human


| INTRODUCTION
Pangolins are mammals that have been facing high extinction risk due to unsustainable levels of poaching for local/regional wildlife markets, exacerbated by international demand from Chinese Traditional Medicine markets (Challender et al., 2020;Heinrich et al., 2016). Although the international trade has been mainly focused on Asian species, it has recently shifted to African pangolins due to the severe decline of Asian pangolin populations (Hu et al., 2020). During the last decade (2010-2019), more than 500,000 individuals of African pangolins were illegally traded internationally (Challender et al., 2020), with the white-bellied (Phataginus tricuspis) and giant (Smutsia gigantea) pangolins being the most traded species (Challender et al., 2020;Emogor et al., 2021;Ingram et al., 2019). In addition to the deleterious effects of the international trade (Challender et al., 2020;Ingram et al., 2019;Zanvo et al., 2021), African pangolins have undergone unprecedented degradation of their habitats, especially in western Africa where >80% of the forest cover vanished over the last century (Aleman et al., 2018). In light of this concerning situation for pangolins, it is likely that their populations are declining in many parts of West Africa and there is a paucity of quantitative data for assessing temporal trends in population size or occupancy and, in turn, identifying areas for prioritizing conservation actions. An exception is Benin, where a recent study (Zanvo et al., 2020) reported declining trends in distribution area for two pangolin species.
Endangered on the IUCN Red List of Threatened Species (Nixon et al., 2019;Pietersen et al., 2019), the whitebellied and giant pangolins are respectively classified as Endangered and Critically Endangered on the Red List for Benin (Akpona & Daouda, 2011). They are threatened by local trade in bushmeat and traditional medicine markets (Zanvo et al., 2021), and by international trafficking Zanvo et al., 2021). Local Ecological Knowledge (LEK)-based surveys suggested declining trends in distribution area for both species in Benin. The white-bellied pangolin underwent a contraction of its occurrence range by one third and the giant pangolin experienced a quasi-total extirpation with 93% of its occurrence range lost over the last two decades (Zanvo et al., 2020). Habitat degradation and overexploitation were the major drivers of population decline according to local hunters (Zanvo et al., 2020). The two species are currently distributed in human-dominated landscapes within patchily distributed habitats (Zanvo et al., 2020), and this habitat fragmentation is likely to impact their persistence, as reported for other mammals (see Horv ath et al., 2019;Schneider, 2001).
Understanding drivers of species distribution across space and over time is an essential step for establishing spatial conservation priorities (Addison et al., 2013;Mouquet et al., 2015), particularly in a context of rapid anthropogenic habitat changes, as observed in Benin (Djaouga et al., 2021;Mama et al., 2013;Padonou et al., 2017). In this context, the long-term persistence of pangolins may be negatively affected by infrastructure development (e.g., road networks), deforestation, and land cover changes, but also positively by the natural permeability of the landscape to dispersing individuals (e.g., watercourses). Given that overexploitation is a likely driver of pangolin decline in Benin (Zanvo et al., 2020) where wildlife markets provide economic incentives to local hunters (Zanvo et al., 2021), the persistence of a population may depend on its accessibility to hunters (Brodie et al., 2015) and/or the level of the demand (e.g., proximity to a wildlife market). Within such a matrix of degraded habitats, protected areas should act as refuges for wildlife (Macedo et al., 2019;Pacifici et al., 2020). However, the efficiency of protected areas in Benin has not really received any attention with regard to pangolin conservation.
Despite high conservation concern for pangolins in Benin (Akpona & Daouda, 2011), little is known about key features of their ecology (Akpona et al., 2008;Zanvo et al., 2020). This hampers informed decision-making, essential to reverse the ongoing declining trends observed in pangolins. Pangolins are cryptic and elusive animals and using conventional monitoring methods to assess population status across space and over time is therefore challenging (Willcox et al., 2019). This is even more difficult in developing countries, where scientific research has been facing financial constraints (Field et al., 2005;Gaillard, 2010) and long-term scientific data are scarce (Early-Capistr an et al., 2020). In such context, it becomes important to find innovative and affordable solutions to survey endangered species such as pangolins. LEK is one such alternative approach recently used to study Asian and African pangolins (Nash et al., 2016;Zanvo et al., 2020). Some authors have raised concerns about the uncritical use of LEK in environmental and conservation studies (Davis & Ruddle, 2010;McKelvey et al., 2008). When population status is assessed via LEK, it is therefore important to validate the approach by comparing its output with independent evidence of occurrence or population size, particularly when modeling is involved (Bélisle et al., 2018;Gilchrist et al., 2005). Our study combines LEK with the presence of pangolin scales as independent evidence of species occurrence.
In Benin, the decline of the two pangolin species has been well documented from a previous study (Zanvo et al., 2020), but this general pattern can mask spatial heterogeneity in population trends and extirpation/ colonization dynamics. The aim of this study is to understand why some pangolin populations have been extirpated within a period of two decades , while the others are still extant. This question is of paramount importance to design future local conservation actions and to assess the effectiveness of the current reserve network. Particular attention has been given to temporal landscape changes, to test the hypothesis that habitat degradation is one of the major divers of pangolin decline in Benin, an assertion frequently expressed in local interviews (Zanvo et al., 2020).
The extirpation status of the focal populations come from the LEK survey conducted by Zanvo et al. (2020) and so it is important to know whether reliable statistical models can be built from such data. Fortunately, an important result of our study is to show that modeling extirpation from LEK data accurately predicted population persistence in an independent set of populations for which pangolin presence is known with certainty. We explored the likelihood of population extirpation in relation to landscape changes and indirect indices of anthropogenic pressure to account for the two main factors identified in pangolin population decline according to interview surveys (Zanvo et al., 2020). The effectiveness of protected areas in preventing or mitigating population decline has been tested as well. It is well-established that the persistence of a population increases with its size; therefore, the contribution of LEK-derived initial population size (abundance index at 1998) has been accounted for in the models for both giant and white-bellied pangolins.

| DATA COLLECTION
We conducted LEK surveys in 312 villages across Benin from April 2018 to March 2019 and the collection of independent evidence of presence (scales) was extended to June 2020. The villages were selected within 25 Â 25 km grid cells (see Zanvo et al., 2020). We received verbal consent of local authorities, but also of each local hunter involved in the focus groups after we had explained the objectives of the study. We gave full assurance of anonymity to each participant. We conducted surveys using: (i) a poster illustrating 12 different species, including 9 species native to Benin (including the two target pangolins) and three species not native to Benin to not influence the respondents' answer about the target species and (ii) an identification guide including different body parts of the four African pangolins to ease the species identification. We carried out the data collection using focus groups (8-15 people) with local hunters following three steps: (i) poster-based presence/absence of pangolins within the local environment, (ii) identification of African species occurring in the local environment using the identification guide of pangolin species, and (iii) local hunter's perceptions of pangolin abundance prior and up to 1998 and at the time of the interviews in 2018 using a pre-established questionnaire (Zanvo et al., 2020). We selected these periods of time in order to reduce the bias related to the LEK degradation over time (Aswani et al., 2018) but also these periods of time were included in the decades of significant deforestation rate in Benin (Padonou et al., 2017). During interviews, hunters provided their perceptions of pangolin abundance at both time periods using the following: "missing"-probability of observing an individual during the target period is zero in the local environment; "scarce"-probability of observing an individual during the target period is near to zero in the local environment; "medium"-high probability of observing an individual every month in the hunting season (dry season in Benin) for the target period; "abundant"-high probability of observing many individuals every month during the target period; and "highly abundant"-high probability of observing many individuals every day during the target period, in the local environment (e.g., municipal or legislative election). For the initial abundance (1998), we referred to a local event at the target period of time with the help of our local guides who at completed at least the primary school and were able to located over time local event and make it understandable to local hunters.
To ensure the reliability of local hunter responses, direct evidence of recent pangolin presence (scales collected after 2013) was collected (if possible) with focus group participants, and also with non-participants identified as holders of direct evidence within each sampled village using the snowball technique. We used the presidential election in 2016 as reference time period with the help of participants to accurately determine the time period when each direct evidence was harvested. Participants and direct evidence donors were men and adult hunters (>18 years old, with at least one experienced hunter >50 years old), and volunteers with at least 10 years residency in the village (Kotschwar Logan et al., 2015). They were selected through informal interviews with villagers.
To assess the indirect influence of spatial configuration of the wildlife markets on the persistence of pangolins in Benin, we systematically and simultaneously georeferenced all the wildlife markets occurring in the districts where we were able to record presence, either historic (pre-1998) or contemporary of at least one pangolin species. We applied a snowball technique to identify all the wildlife markets in each district.

| Model, variables, and tested hypotheses
From the original survey data, we discarded the villages where pangolins have never occurred, and the villages in which both species had already disappeared in 1998 according to local hunters. Among the 312 villages surveyed by Zanvo et al. (2020), we selected those (N = 156) with at least one species present by 1998 ( Figure 1). We recorded presence of white-bellied pangolin in all 156 villages and giant pangolin in only 52 villages. Of the 156 villages where white-bellied pangolin was recorded through LEK, 52 of them also provided direct evidencebased observations. We used LEK data (156 observations combining both species) to calibrate the model, whereas direct evidence-based observations (52 observations) acted as a validating set. The modeled variable was binary: 0 if a pangolin population had been extirpated (absence in 2018), 1 if it has persisted until 2018 (presence in 2018). Note that scales are evidence of recent presence (individuals hunted between 2016 and 2020) and cannot be used to validate earlier, historical reported presence. We considered abundance levels in 1998 and 2018 as the initial and current abundances of pangolins, respectively, for running our models. To model the probability of persistence (Pr2018), the following explanatory variables were considered in order to test hypotheses related to human exploitation, landscape changes, initial population abundance, and interspecific differences.
Species: this binary variable referred to our focal species: the white-bellied pangolin (coded 0) versus giant pangolin (coded 1). Large-bodied species are more vulnerable to increased mortality rates (as induced by human exploitation in our case) than small-bodied species and as a result, they are more likely to decline, and eventually to be extirpated. Indeed, large body size is positively correlated with extinction risk in mammals (Liow et al., 2008;Lyons et al., 2004). The giant pangolin is by far the largest of our focal species: it is 12-14 times heavier than the white-bellied pangolin (see Kingdon & Hoffmann, 2013). It would have been preferable to model each species independently, but there is not enough recent occurrence data for the giant pangolin to run a reliable model for this species alone ( Figure 1b). Analyzing the two species jointly allows comparing their average extirpation probabilities but makes the debatable assumption that, otherwise, the two species are sufficiently close in their habitat requirements to have identical responses (regression coefficients) to the tested variables. The two species are sufficiently close in their habitat requirements so that it is very unlikely that they display an opposite response to an important variable (e.g., a positive correlation with deforestation for one species and a negative one for the second). However, observing a significant species effect in the model could be the result of different degrees of sensitivity to habitat changes.
ab1998 represents the initial (circa 1998) population abundance index as given by LEK surveys. This is an ordinal variable scaled as follows: 1 (scarce), 2 (medium), 3 (abundant), and 4 (highly abundant). All else being equal, the persistence of a population is positively related to its size (Lande, 1993), irrespective of the nature of extinction, whether it is deterministic (extinction resulting from a declining population size) or stochastic (extinction resulting from the temporal fluctuation in population size). Therefore, we hypothesized that the current persistence of a pangolin species in a local environment is strongly dependent on its abundance level two decades ago (1998).
distPA is the Euclidian distance (km) from each village to the nearest protected area under the management of government officials. If the protected areas are efficient in maintaining viable populations of pangolins, we hypothesized that the probabilities of persistence of the white-bellied and giant pangolins are higher near protected areas, acting as sources for emigration.
distRoad is the Euclidian distance (km) from each village to the nearest asphalted road. Given that previous studies reported correlations between road density and the high extraction of wildlife (Constantino, 2016;Espinosa et al., 2014), we hypothesized a negative relationship between distRoad and the probability of a pangolin population to still be extant by 2018.
distMarkets is the Euclidian distance (km) from each village to the nearest wildlife market (bushmeat and traditional medicine markets). The northern markets are drawn from Djagoun et al. (2013), the others were identified during our LEK survey. Wildlife markets provide economic incentives to local hunters in tropical Africa (see, Ingram et al., 2021) that could exacerbate, in their occurrence areas, the decline of taxa with the highest commercial values, such as pangolins. We therefore hypothesized that the probability of persistence increases for both pangolin species with the distance to the nearest wildlife market.
distWater is the Euclidian distance (km) from each village to the nearest main watercourse (permanent). The white-bellied and giant pangolins may occur in dense woodlands with nearby water courses (Kingdon & Hoffmann, 2013;Zanvo et al., 2020). We hypothesized that the probability of spatial persistence of both species is higher near watercourses.
dft i is the total deforestation area (km) around each sampled village within a buffer size of radius i (see below) between 2001 and 2019. The deforested areas were extracted from the Global Forest Change database (https://google.earthengine.app/view/forest-change; Hansen et al., 2013). For each buffer size, we clipped for each village the extent of cumulated deforested areas from 2001 to 2019 using ArcGis 10.8.1 (Esri France); 2001 is the oldest year available, while being sufficiently close to 1998 to provide a meaningful trend over our period of interest. The ecological requirements of white-bellied and giant pangolins are strongly related to forest habitats, particularly primary and secondary forests with large trees (see details in the metadata computation section; Supplementary information). Accordingly, we hypothesized that the strength of population decline, and in turn the probability of extirpation, is positively related to the rate of deforestation.
LSI i is the change in a habitat suitability index between 1998 and 2015 within a buffer of radius i. The suitability index is based on the land-use categories included in the Global Land cover database (http://maps. elie.ucl.ac.be/CCI/viewer/download.php). According to the literature on pangolin habitat requirements (Akpona et al., 2008;Angelici et al., 1999;Fischer et al., 2002;Kingdon & Hoffmann, 2013;Pages, 1968;Pagès, 1975;Zanvo et al., 2020), each land-use category received a suitability score (Table S1) that we used to compute a suitability index, SI, over the area of interest (see Supplementary information for details). We applied the same scoring to the white-bellied and giant pangolins in our study area, given that LEK-based past and current distribution revealed an overlapping of occurrence in habitat for both species (Zanvo et al., 2020). For each village, SI was computed within a buffer of radius i for 1998 (SI i,1998 ) and 2015 (SI i,2015 ) to obtain the land-use suitability change index over time (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) LSI i = SI i,1998 À SI i,2015 . A positive value indicated that the suitability of the landscape increased over time, and a negative one the opposite. Land-use data were not available after 2015, but the latter is sufficiently close to 2018 to provide meaningful trends over our period of interest. We hypothesized a positive correlation between pangolin population trends and habitat suitability trends and therefore that population extirpation would be more likely if LSI is negative.
Two of the explanatory variables ( dft i and LSI i ) depended on a buffer size for their computation. The buffer was the area defined by a circle centred on the village of interest. Five buffers of radius i have been considered: 5, 10, 15, 20, and 25 km. The choice of a minimum 5 km buffer size referred to 7.7 km mean hunting effort previously determined in West Africa (Alexander et al., 2015), and the maximum 25 km buffer size was based on the extent of our sampling unit (25 Â 25 km).

| Statistical analysis
Correlations between variables (Table S2) were computed and, as only small or moderate values were observed (<0.4), coefficient estimates are unlikely to have been biased by collinearity. The following generalized linear model with binomial error and logit link was fitted to the calibration data set for each buffer size: Pr2018 $ species + ab1998 + distPA + distRoad + distMarkets + distWater + dft i + LSI i : We ran for each buffer size i (5, 10, 15, 20, and 25 km radius length) a model counting for both two species through the "species" variable. Based on the Akaike's Information Criterion (AICs) of models, we selected the buffer size providing the best fit.
For 52 villages, direct evidence (scales) of the recent presence of white-bellied pangolin was available allowing us to test the accuracy of the model predictions based only on LEK observations (n = 156). To do this, we estimated the probability of persistence of white-bellied pangolin for each village where scales were available using the model retained in the above step and setting the variable "species" at zero, where values >0.5 predicted persistence, while values <0.5 predicted local extinction. Even if we were not able to test the reliability of LEK-based model prediction for the giant pangolin due to the lack of direct evidence, the LEK was collected simultaneously from the same hunters and following the same protocol for both white-bellied and giant pangolins during surveys. We then consider that the LEK was accurate for both species.
By combining the calibration (n = 156) and the validation (n = 52) sets, we then fitted the model using all N = 208 observations, to assess whether the explanatory variables contributed significantly to the model and acted in the anticipated way. We also fitted a second regression model using only the observations of white-bellied pangolin (N = 156 observations).
Based on the modeled trends of species status over the last 20 years, projections were made for the next 20 years following two possible scenarios of landscape changes. In both cases, population abundance is updated by the recorded abundance in 2018 (e.g., ab2018). The first scenario assumes that the landscape will not change in the future, meaning no more deforestation (dft = 0) and no land-use changes (LSI = 0). The second scenario assumes that the observed trends over the past 20 years will be maintained for the next 20 years. This implies that the deforestation and land-use changes will continue and the same values of dft and LSI for the last 20 years (1998-2018) will be obtained for the 2018-2038 period. All the other predictors are kept un-changed. We used the values of persistence probability to calculate the number of occurrence villages globally and per quartile (lower quartile = first 25% predicted probability of persistence and upper quartile = last 25% predicted probability of persistence) in 2018 and 2038 in order to visualize the trends over time and the variability among villages.
Statistical analyses and mapping of distributions (occurrence in 2018 and the persistence probability of pangolins) were carried out in R version 4.0.5 (R core development team 2021) and ArcGis 10.8.1 (Esri France), respectively.

| RESULTS
Model selection considering buffer size revealed that the binomial GLM model developed using a 15 km buffer had the lowest AIC value (Figure 1; Supplementary  information).
The LEK-based model was highly accurate in predicting presence of white-bellied pangolin across the villages where actual presence was also confirmed with scales. Only one village was misclassified, therefore resulting in 97% correct classification (Figure 2; Supplementary information).
In the model using all the data (N = 208), all the variables that contributed significantly (p < .05) to the fit acted as we anticipated (regression coefficients of the anticipated sign), and only two variables were not significant (distMarkets and distWater; Table 1). The variables contributing the most to the model were species identity and initial abundance (both with p < .001). Among 156 occurrence villages of white-bellied pangolin in 1998, our model predicted a disappearance of the focal species in 33% of them in 2018 against 69% and 63% following the scenarios, future change mirroring change in the 1998-2018 time period and no future change respectively by 2038 (Figure 2). In the upper quartiles, almost all the villages (85%-90%) will continue to be occurrence areas of white-bellied pangolin by 2038 unlike to the villages in the lower quartiles where an almost total disappearance of the focal species will be observed throughout and whatever the scenarios. The probability of persistence for the white-bellied pangolin displayed strong spatial heterogeneity with higher values of probability of persistence (.586-1.0) in villages close to protected areas (upper quartiles) and the lower values of probability (.000-.011) in villages relatively far from protected areas (lower quartiles) in 2018 and 2038 (Figure 3). In contrast, for the giant pangolin, persistence probability was low throughout the area studied (Figure 2b). The same interspecific difference in spatial variability was observed when persistence was predicted for 2038. For both species, assuming  either a status quo in habitat suitability or a habitat change continuing at the ongoing rate had a negligible effect on the predicted decline ( Figure 2). The occurrence areas of the two species was predicted to decrease in the coming years, and the projection was particularly alarming for the giant pangolin, which was predicted to be virtually extinct from Benin by 2038. These results were qualitatively unchanged when the model was fitted to the white-bellied pangolin alone (Table 1). Two hotspots of high presence probability were observed: near the forest patches between 6 5 0 0 00 and 7 5 0 0 00 (Lama forest reserve and the sacred grove Gnanhouizoumè) and a central block of protected areas located between 8 5 0 0 00 and 10 0 0 0 00 (Mont kouffé-Wari Maro-Ouémé supérieur forest reserves) (Figure 3). The spatial configuration of predicted persistence is almost unaffected by choosing one scenario of landscape change over the second one. Predicted spatial patterns of persistence based on the model fitted with only the white-bellied pangolin were virtually the same as those obtained with the two species model ( Figure S3). As expected from the high rate of decline observed in Figure 2, the predicted probabilities of persistence until 2038 of the giant pangolin in the two localities where it was reportedly present in 2018 is almost nul (p < .000001), whatever the scenario of landscape change.

| Spatio-temporal patterns in extirpation rates
Effective conservation implies informed decision-making, which is not an achievable objective without relevant data. African pangolins and, in particular, the giant pangolin in Benin (Zanvo et al., 2020) is a data poor species at great extinction risk. Although LEK is subject to some controversies, since its international recognition by the Convention on Biological Diversity, its combination with scientific knowledge is increasingly encouraged in all disciplines of ecology, including modeling (Bélisle et al., 2018) to guide conservation decision-making. Our study, using LEK combined with direct evidence, partly fills an ecological knowledge gap and provides data about the spatio-temporal dynamics of pangolin populations that should prove useful to conservation management ensuring the long-term persistence of pangolins in Benin.
The GLM model revealed a high probability of pangolin population persistence in localities where abundant populations were reported two decades ago. This result is coherent with the expectation that small populations are more likely to go extinct than large ones, whatever the extinction processes and drivers.
Considering the model structure (no interaction between species and other variables), the low number of giant pangolin occurrences, and the validation test performed for the white-bellied pangolin only, detailed inter-specific comparison of model outputs should be done with caution. However, our study empirically supports the prevailing notion that giant pangolin populations have been extirpated at a higher rate than whitebellied pangolin populations. This is an expected pattern considering the positive correlation between body size and extinction risk in mammals (Liow et al., 2008;Lyons et al., 2004). Several factors can explain the sensitivity of large species to human impacts, such as they generally occur at low density and have low reproductive rates (Ripple et al., 2017). In our analysis, population abundance was accounted for; therefore, our results point out a higher sensitivity of the giant pangolin to hunting. Besides body size, differences in life-history features between species may explain interspecific variability in rate of decline. For instance, it has been suggested (Newton et al., 2008) that in Vietnam, terrestrial pangolin species faced a greater risk of hunting than more arboreal, elusive species. As the giant pangolin is terrestrial and the white-bellied pangolin is arboreal, our results are coherent with this hypothesis. Besides a higher sensitivity to hunting, giant pangolin might decline at a higher rate than white-bellied pangolin because of a higher sensitivity to habitat changes. This could suggest that giant pangolin is more specialized in its ecological t requirements than the white-bellied pangolin. According to our model, the probability of persistence of a pangolin population increases with distance to roads. In Benin, pangolins are subject to heavy poaching pressure, in particular to supply local trade occurring in wildlife markets , and the negative effect of an asphalted road network on the persistence of pangolins is in line with previous studies, which demonstrated that road development is a driver of high extraction of wildlife in natural habitats (Constantino, 2016;Espinosa et al., 2014). Asphalted roads might also facilitate access to pangolin populations and, hence, increase hunting pressure (Brashares et al., 2001). Another possible scenario is that road density is an indirect measure of other human perturbations (e.g., habitat fragmentation), which could negatively affect the persistence of pangolin populations. The combined effects of the above scenarios would likely have contributed to the extirpation of pangolin populations in areas with heavy road density. However, the impact of roads on biodiversity is a combination of varied direct and indirect effects, and noise disturbance (Ibisch et al., 2016) might be involved in observed declines.
Contrary to our expectation, proximity to wildlife markets had no effect on the probability of population persistence. For pangolins, in contrast to ungulates , the proximity to wildlife markets did not appear to increase hunting pressure. In support of this observation, studies have reported that in Benin, it is not uncommon for pangolins to be transported over long distances (200-300 km) before being sold, which may blur more localized impacts . It is known that in Benin, pangolins are subject to heavy poaching pressure, in particular to supply local trade occurring in wildlife markets (Zanvo et al., 2021). The proximity of asphalted roads may facilitate the access of hunters to markets where they can sell pangolins. This could reflect the complex relationship existing between a wildlife trade system and its surrounding landscape, as observed for instance in Ghana (McNamara et al., 2015).
Our analyses pinpointed the negative impacts on pangolins of deforestation and reduction of suitable habitats over the last two decades, confirming that habitat degradation has contributed to the decline of pangolin populations. Other studies have also pointed out that land-use changes led to local extinction in mammals (Magioli et al., 2021;Visconti et al., 2011). Given the ecological requirements of white-bellied and giant pangolins, strongly linked to tree-dominated habitats (Akpona et al., 2008;Hoffmann et al., 2020;Jansen et al., 2020;Pages, 1968;Pagès, 1975;Zanvo et al., 2020), the fast change of land-use in Benin, which converted 53% of forests, woodlands, and tree savannas into farmlands and degraded savannas in just two decades (1990-2010Padonou et al., 2017), could explain pangolin population decline. An evidence-based ecological relationship between habitat and species loss has been demonstrated elsewhere (Horv ath et al., 2019). In our model, we were not able to assess whether or not this threshold has already been reached in some areas, and this point deserves further attention. However, the fact that the models predict for the next two decades almost the same decreasing trends with and without deforestation and land-use changes ( Figure 3) suggests that habitat degradation is probably not the primary factor of recent and anticipated population declines. This result does not entirely support local knowledge/perception, which ranked deforestation as the primary driver of pangolin decline in Benin (Zanvo et al., 2020). However, considering the significant deforestation rate and land use change observed in Benin during the 20th century (Djaouga et al., 2021;Padonou et al., 2017), it is likely that habitat degradation and defaunation of large-bodied animals (rare) changed the hunting pattern that has recently been focusing on medium and small animals as pangolins (see Djagoun et al., 2022). The latter probably exacerbated by recent international demand in pangolins in West Africa (Emogor et al., 2021). Molecular-based inference of demographic history of white-bellied pangolins in the Dahomey Gap suggested that the species experienced a drastic reduction of its effective population size at least two centuries ago . The possibility that the temporal changes in pangolin habitat suitability have been imperfectly captured in our approach cannot be entirely discarded given our models were based on LEK subject to some biases and the sampling unit extent (25 Â 25 km) was approximated during the focus groups. However, in light of the available information about pangolin habitat requirements, it seems unlikely that major landscape features have been overlooked in our habitat suitability index and future efforts should be devoted to better modeling of micro-habitat requirements.
Apart from the above-mentioned deleterious spatial factors for pangolins, our model showed significant probability of persistence near protected areas for the whitebellied and giant pangolins in Benin. Not surprisingly, this result supports previous studies that have highlighted the cornerstone importance of protected areas in biodiversity conservation (Macedo et al., 2019;Pacifici et al., 2020). Forest reserve patches have been the last strongholds for the conservation of pangolins in Benin: mitigating the combined effects of habitat changes (road development, deforestation, land-use changes) and overhunting. The fact that no significant effect of proximity to watercourses (a proxy for forest gallery) was observed, suggests that proximity to gallery forests is not a sufficient condition to prevent populations from declining even if many local hunters evoked recent observations near watercourses within protected areas. Some level of protection seems to be required as well. The role of protected areas is more perceptible in the Lama forest reserve (6 5 0 0 00 and 7 5 0 0 00 ), which is the best managed forest reserve in Benin and in the central block of protected areas, including the Mont Kouffé-Wari Maro-Ouémé supérieur forest reserves (8 5 0 and 9 30 North) established on an uneven relief making difficult access to certain marginal lands (mountains and rocks), unsuitable for agriculture. It is also worth noting that, at the northern limit of the occurrence range of the white-bellied pangolin in Benin (circa 11 N, Zanvo et al., 2020), the species persists mainly in or at the vicinity of Alibori forest reserve, emphasizing the efficiency of protected areas in sustaining populations when unprecedented degradation of natural landscape occurs (see Arouna et al., 2016). Protected areas are not a panacea in preventing population extirpations (Macedo et al., 2019;Pacifici et al., 2020), particularly if they are too small or not well managed (Brashares et al., 2001). In case of management gap, the resistance of protected areas to anthropogenic landscape changes (conversion of forests, woodlands, and tree savannas into farmlands; Padonou et al., 2017) and poaching (see Zanvo et al., 2021) could be compromised. In Benin, forest reserves have a large extent of management gaps so poaching of pangolins occurs regularly in protected areas, emphasizing that such areas only provide partial protection. In addition, it is possible that some forest reserves are too small to ensure long-term viability of pangolin populations, even in the absence of any detrimental human impacts.

| Caveats
We have used direct evidence of presence (scales) to validate LEK data to understand the drivers of spatial persistence of the white-bellied and giant pangolins in Benin. These results demonstrate that LEK could be a valuable decision-informing tool in the management of wildlife (Bélisle et al., 2018;Gilchrist et al., 2005). The land-use suitability change index over time is among the most important variables that explain the persistence of pangolins in Benin, and therefore it would be helpful to improve its reliability using a scoring approach that integrates both the type of ecosystems and the microhabitat requirements of each species. In contrast with occurrence data, limited information exists about pangolin abundance in Benin (Akpona et al., 2008); and, more generally, assessing population status of pangolins is challenging (Willcox et al., 2019). Therefore, designing a protocol based on direct observations to validate LEK-based abundance indices is expected to be faced with many constraints. In the same vein, and because it is likely that population density and overexploitation are negatively correlated (Brodie et al., 2015), it would be important to have a direct estimate of hunting pressure on pangolin populations, and, if possible, the use of less indirect proxies (e.g., encounter rate in occurrence habitats although acquiring reliable data on these topics is challenging).

| Management implications
Our results strongly suggest that the combined effects of overexploitation and human-driven landscape changes have resulted in numerous pangolin population extirpations. In contrast to the likely unavoidable extinction of giant pangolins, the white-bellied pangolin is declining but at a spatially heterogeneous rate so that the great majority of the population still benefits from supportive conditions within protected areas, and it is expected to persist over the next decades. However, even halting deforestation and habitat degradation will not be enough to prevent the decline of white-bellied pangolin populations. The projected probability of persistence is low for several populations. In spite of this, the populations with the highest predicted persistence over time should be considered for reliable spatial prioritization, an essential step for long-term conservation planning (Grantham et al., 2020). In particular, management efficiency of protected areas with management gaps should be improved by reinforcing daily combat against poaching, alongside rigorous law enforcement against illegal trade to help reduce local/international trade of pangolins. Habitat protection will also help improve microhabitat requirements and substantially increase the possibility of longterm conservation of pangolins in Benin.

ACKNOWLEDGMENTS
We are grateful to the Public Forest Service of Benin for providing the research permits to conduct the fieldwork. We thank the Beninese local authorities who gave us a prior consent and without whom this study would not have been possible. We are grateful to all the local people who participated in the surveys. We thank Helen Nash for her useful revision of an early draft of the manuscript. Data collection was funded by RADAR-BE Project, Jeune Equipe Associé à l'IRD. Stanislas Zanvo was financially supported by the ARTS-IRD program and Rufford Foundation.

DATA AVAILABILITY STATEMENT
Data collected and generated during this study are available as supplementary material to this paper.