Understanding historic patterns of land use and land cover change across large temporal and spatial scales is critical for developing effective biodiversity conservation management and policy. We quantify the extent and fragmentation of suitable habitat across the continental range of Asian elephants (Elephas maximus) based on present-day occurrence data and land-use variables between 850 and 2015 A.D. We found that following centuries of relative stability, over 64% (3.36 million km²) of suitable elephant habitat across Asia was lost since the year 1700, coincident with colonial-era land-use practices in South Asia and subsequent agricultural intensification in Southeast Asia. Average patch size dropped 83% from approximately 99,000–16,000 km² and the area occupied by the largest patch decreased 83% from ~ 4 million km² (45% of area) to 54,000 km² (~ 7.5% of area). Whereas 100% of the area within 100 km of the current elephant range could have been considered suitable habitat in the year 1700, over half was unsuitable by 2015, driving potential conflict with people. These losses reflect long-term decline of non-forested ecosystems, exceeding estimates of deforestation within this century. Societies must consider ecological histories in addition to proximate threats to develop more just and sustainable land-use and conservation strategies.

Elephant occurrence data

Elephant occurrence locations were initially compiled from the Global Biodiversity Information Facility (https://www.gbif.org/), Movebank (https://www.movebank.org/) and published literature as well as data contributed by the authors based on direct sightings, data logged via tracking devices, and camera traps (n>5000 locations). Records were first checked visually for irrelevant points (e.g., occurrences outside natural continental range, from GBIF) and then refined to include locations representing ecosystems where the species could conceivably flourish, including but not exclusively limited to protected areas. To minimize sampling bias that could result in model overfitting, we further subsampled data to cover the full distribution as widely as possible while eliminating redundant points located within any particular landscape. For instance, thousands of potential redundancies from collar-based tracking datasets were removed by using only one randomly selected data point per individual, per population or landscape. Outliers from the remaining points were removed using Cooks’ distance to eliminate locations that could represent potential errors. The final dataset consisted of 91 occurrence points spanning the years 1996-2015 which served as training data, where all data other than from GBIF and cited literature were contributed by the authors or individuals listed in acknowledgments. QGIS and Google Earth Pro were used to initially visualize and process the data.

Predictor variables

We used the Land-Use Harmonization 2 (LUH2) data products 25 as our environmental variables. The LUH2 datasets provide historical reconstructions of land use and land management from 850 to 2015 CE, at annual increments. The LUH2 data products were downloaded from the University of Maryland at http://luh.umd.edu/data.shtml (LUHv2h “baseline” scenario released October 14th 2016). They contain three types of variables gridded at 0.25° x 0.25° (approximately 30 km2 at the equator): state variables describing the land-use of a grid cell for a given year, transition variables describing the changes in a grid cell from one year to the next, and management variables that describe agricultural applications such as irrigation and fertilizer use, totaling 46 variables. Of these, we selected 20 variables corresponding to all 3 types which were expected to be relevant to elephant habitat use based on knowledge of the species’ ecology 21,22,32,81. Using ArcGIS 10 (ESRI 2017) we extracted each variable between 850–1700 CE at 25-year increments, and annually between 1700–2015.  We separately obtained elevation from the SRTM Digital Elevation Model.

Data analysis

We limited the geographic extent of all analyses to the 13 range countries in which elephants are currently found. We used MAXENT, a maximum entropy algorithm 82, to model habitat suitability using the ‘dismo’ package in R (R Core Team 2017). Resulting raster files were binarized in ArcGIS into suitable and unsuitable habitat with a pixel size of approximately 20 km² as a cutoff threshold. As there is no commonly accepted threshold type 84, to ensure that the specific choice of threshold did not affect the observed trends, we initially used three possible thresholds: 0.237, representing ‘maximum test sensitivity plus specificity,’ 0.284 corresponding to ‘maximum training sensitivity plus specificity,’ and 0.331 representing ‘10th percentile training presence’. Unless otherwise stated, for subsequent analyses we show only results using the threshold of 0.284, where everything below this threshold was classified as ‘unsuitable’ and everything above it was classified as ‘suitable’. The resulting binary maps were re-projected using the WGS84 datum and an Albers Equal Area Conic projection.

Polygons representing the known elephant range were digitized from Hedges et al. 2008 from the category labelled as “active confirmed”. We refer to the areas within these polygons as “current range,” and refer to areas outside them as “potential range”. We compared the total extent of suitable habitat within and outside the current elephant range, quantifying changes over time. Country-level analyses were conducted for all countries except Indonesia and Malaysia where the Bornean and Sumatran ranges were treated separately in recognition of the distinct subspecies in these two regions. We ranked each region based on the percentage of the current range within that region as well as the proportion of the estimated elephant population found within it, and calculated the ratio of these ranks.

We calculated the total change in extent of suitable habitat by subtracting the area of suitable habitat available in 2015 from the area available in 1700, as major changes were observed within this period. We also specifically quantified the percentage of suitable habitat found within a 100 km buffer of the current range polygons in both years. We then calculated fragmentation statistics using the program FRAGSTATS v.4.2 88. These metrics characterize changes to the spatial configuration of habitat in addition to their absolute extent. We used a ‘no sampling’ strategy with the search radius and threshold distance set to 61 km (approximately three-pixel lengths) based on the movement and dispersal capacity of elephants.

See associated paper for references, tables and figures.

Spreadsheet data: Microsoft Excel / Google sheets / Open Office

GIS data: ArcGIS / R / QGIS

Rdata or scripts: R / R Studio / Notepad ++