Using natural history collections to investigate changes in pangolin (Pholidota: Manidae) geographic ranges through time

Georeferencing

For the NHM locality data, all georeferencing was carried out by one of us (EB) using the GBIF Best Practice guidelines (Chapman & Wieczorek, 2006). Each location was found in Google Maps, and decimal latitude and longitude coordinates for the midpoint of the locality were recorded. We also recorded an extent (in m) from the midpoint to the furthest place that could still be considered to be within the locality as a measure of error (i.e., the point-radius method; Wieczorek, Guo & Hijmans, 2004), along with a certainty score (0, 25, 50, 75 or 100%) describing how confident we were that the locality truly fell within the recorded extent. For the GBIF data we followed the same protocol, except where localities had already been georeferenced we used the existing coordinates and just estimated the extent (in m).

Before analyses we checked all records for data quality. We excluded specimens with extents >1,000 km and specimens from zoos or where we suspected the locality data were the shipping, rather than collection, location. Specimens with localities more than 100 km outside their current geographic range were also checked. Where this was a clear-cut taxonomic error we corrected the taxonomy, for example any Manis species found in Sri Lanka must be Manis crassicaudata as no other members of the genus are found there. If there was no obvious taxonomic error we assumed this was a recording error and omitted the specimen from further analyses. Details of excluded and modified specimens are shown in Tables S1 and S2.

To incorporate georeferencing error, for each specimen locality we created a point-radius polygon from the error extent for that locality, using the R packages sf and sfe (Pebesma, 2018; Curry, 2019; Fig. 1). The point-radius polygon is a circle around the locality with a radius equal to the error extent of the locality.

Area of habitat (AOH) maps

We used the known habitat preferences, altitudinal limits and 2019 geographic range polygons all sourced from the IUCN Red List website (IUCN, 2020) to create AOH maps for each pangolin species. Trends in species AOH over time can be useful in obtaining estimates of their population decline under the IUCN Red List criteria and in assessing levels of fragmentations of pangolin habitats (IUCN, 2020; Brooks et al., 2019). No elevation data could be found for Phataginus tricuspis and P. tetradactyla so just the habitat preferences were used to create the AOH maps for both of these species. Minimum altitudinal limits were not present for either Smutsia species, so the value was assumed to be 0 (as is the case with all other species where this data was available; Brooks et al., 2019). We downloaded a 2018 land cover classification raster (i.e., gridded data) at 100m horizontal resolution from the Copernicus Climate Change Service version 2.1.1 (Buchhorn et al., 2019) and an elevation raster from WorldClim (Fick & Hijmans, 2017) at 1 km resolution. Following the methods outlined by Brooks et al. (2019), we used the raster package in R (Hijmans, 2018) to crop both the elevation and land-use rasters for each species using the IUCN range maps to define the extents. We then used the raster calculator in QGIS version 3.14 (QGIS.org, 2020) to create binary rasters using the known habitat preferences and altitude limits of each species, with the following formula (where XX is a species): (1)${\displaystyle \left(\text{“}\mathrm{XX}\text{\_}\mathrm{elevation}\text{”}<0\right)\ast 0+\left(\left(\text{“}\mathrm{XX}\text{\_}\mathrm{elevation}\text{”}>=0\right)\mathrm{AND}\left(\text{“}\mathrm{XX}\text{\_}\mathrm{elevation}\text{”}<=2015\right)\right)\ast 1+\left(\text{“}\mathrm{XX}\text{\_}\mathrm{elevation}\text{”}>2015\right)\ast 0.}$

A value of 1 in the binary raster is therefore suitable habitat and 0 is unsuitable. These binary elevation and land-use rasters were then resampled and multiplied together for each species in R to create the binary suitable/unsuitable AOH maps.

Anthropogenic data

Geographical analyses are often based on summarising data within grid cells of a defined size, then using values for each grid cell as input data for downstream analysis (e.g., Derrick, Cheok & Dulvy, 2020). We downloaded historical and current human population size data as population count and population density from the History Database of the Global Environment (HYDE v3.2.1; Goldewijk et al., 2017) at 5′ or 0.083 degree resolution. For land-use change, we downloaded historical land-use states data at 0.25 × 0.25 degree spatial resolution from the Land-Use Harmonization (LUH2) project (LUH2 v2 h; http://luh.umd.edu/data.shtml; Hurtt et al., 2020). These data represent the fraction of each 0.25 × 0.25 degree grid cell covered by each land-use type. The land-use types we focussed on were forested primary land (primf; all forested areas including scrublands and savanna, previously undisturbed by any human activities post 850), non-forested primary land (primn; all non-forested areas including grasslands and wetlands, previously undisturbed by any human activities post 850) and urban land (urban). Both HYDE and LUH2 data were based on the high data-driven land-use reconstructions from HYDE (Goldewijk et al., 2017).

For the human population size and the land-use data, we calculated the change in population size or land-use cover (but not both) for each 0.25 × 0.25 degree grid cell between 2015 (present-day) and 1850, 1900 or 1950 in turn. We used 2015 as our present-day baseline as these were the most recent data available (Goldewijk et al., 2017; Hurtt et al., 2020), and used 1850 as our earliest historical date as most specimens were collected between 1850 and 1950 (ideally we would extract values for the year each specimen was collected, but we do not have year data for all specimens). For each specimen and time period, we then extracted the mean change in every population size or land-use variable across the 0.25 × 0.25 degree grid cells occupied by the specimen’s point radius polygon. We used these values in our analyses below. Note that specimen MVZ_MVZ:Mamm:125554 was excluded from these analyses because it comes from an extremely small island on the P’eng-hu peninsula and there were not sufficient land-use or human population data within this locality for overlaps to be calculated.

All data required to repeat the analyses are available on the NHM Data Portal (data.nhm.ac.uk; Buckingham et al., 2019).

Correlates of overlaps with AOH maps

For each species of pangolin we calculated overlaps between specimen localities and present-day AOH polygons in two ways (Fig. 1). (1) Percentage locality overlap: the percentage of specimen localities (including their error extents) that fall within their species AOH polygons; (2) Percentage area overlap: the mean percentage area of specimen point-radius polygons that overlap with AOH range polygons. We calculated these percentages for specimens with certainty scores ≥ 50% and extents <50 km with and without duplicates, and for all specimens with certainty scores >0%, with and without duplicates. We also grouped species by continent, ecology and IUCN Red List status and calculated the percentage overlaps for these groupings. Note that Smutsia temminckii had an extremely complex AOH map because of its broad range, meaning that we were unable to extract overlaps from it due to computational issues (lack of available memory). We therefore used the IUCN map for this species to estimate overlaps rather than its AOH.

To test for correlates of overlaps among specimen localities and present-day range maps, we first fitted binomial generalised linear models (GLMs), with the number of specimens whose localities (including their error extents) were within their species AOH polygons (success) and the number of specimens whose localities were not within their species AOH polygons (failure) for each species as the response variable (i.e., a binomial response where the number of successes and failures were jointly modelled). We used the collection year, continent, whether the species was primarily terrestrial or arboreal/semi-arboreal, and IUCN Red List status, as predictors, and used standard model checks for GLMs (Q-Q plot, histogram of residuals, residuals vs. linear predictors, response vs. fitted values) to assess model fit. These analyses have very low power; most have only eight data points, one for each species, except the collection years analysis which has 106 data points, one for each year × species combination.

Next, we repeated these analyses using percentage area overlap for each specimen as the response variable. We used the species, collection year, continent, whether the species was primarily terrestrial or arboreal/semi-arboreal and IUCN Red List status, as predictors (note that we could not fit multiple regressions because these resulted in rank deficient matrices because species can be written as a combination of the other predictors). Because percentage area overlap is a proportion, we fitted beta regression models after first transforming percentage area overlap to remove 0s and 1s using the following equation (Douma & Weedon, 2019): (2)${\displaystyle x_i^t=x_i\left(n-1\right)+0.5∕n}$

where $x_i^t$ is the transformed value of x_i, and n is the total number of observations in the dataset. We then fitted our models using the betareg function in the betareg R package (Cribari-Neto & Zeileis, 2010) using a logit link function and with both fixed and variable precision (ϕ). Estimators from beta regression can be biased, particularly at low sample sizes (Douma & Weedon, 2019; Grün, Kosmidis & Zeileis, 2012), so we repeated our models using the bias correction and bias reduction methods built into betareg (Grün, Kosmidis & Zeileis, 2012). We compared these results to our original models to ensure that bias was not influencing our conclusions.

Anthropogenic drivers of range change or contraction

If species have changed their ranges in response to anthropogenic changes in their habitats, we expect to see negative correlations between measures of overlap in past and present-day ranges and changes in anthropogenic variables. To test this, we used beta regression as described above to determine whether percentage area overlap was correlated with changes in human population size or land-use, between 2015 and 1850, 1900, or 1950 in turn. For human population size we used changes in log population count and log population density as predictor variables. For land-use change we investigated changes in the percentage cover of forested primary land, non-forested primary land and urban land as predictor variables.