Humans have reshaped the distribution of biodiversity across the globe, extirpating species from regions otherwise suitable and restricting populations to a subset of their original ranges. Here, we ask if anthropogenic range contractions since the Late Pleistocene led to an under‐representation of the realized niches for megafauna, an emblematic group of taxa often targeted for restoration actions. Using reconstructions of past geographic distributions (i.e., natural ranges) for 146 extant terrestrial large‐bodied (>44 kg) mammals, we estimate their climatic niches as if they had retained their original distributions and evaluate their observed niche dynamics. We found that range contractions led to a sizeable under‐representation of their realized niches (i.e., niche unfilling). For 29 species, more than 10% of the environmental space once seen in their natural ranges has been lost due to anthropogenic activity, with at least 12 species undergoing reductions of more than 50% of their realized niches. Eighteen species mainly now be confined to low‐suitability locations, where fitness and abundance are likely diminished; we consider these taxa climatic refugees. For those species, conservation strategies supported by current ranges risk being misguided if current, suboptimal habitats are considered baseline for future restoration actions. Because most climate‐based biodiversity forecasts rely exclusively on current occurrence records, we went on to test the effect of neglecting historical information on estimates of species’ potential distribution – as a proxy of sensitivity to climate change. We found that niche unfilling driven by past range contraction leads to an overestimation of sensitivity to future climatic change, resulting in 50% higher rates of global extinction, and underestimating the potential for megafauna conservation and restoration under future climate change. In conclusion, range contractions since the Late Pleistocene have also left imprints on megafauna realized climatic niches. Therefore, niche truncation driven by defaunation can directly affect climate and habitat‐based conservation strategies.

To evaluate the effect of shifted baselines on estimates of sensitivity to future climate change—measured in terms of potential distribution variation—we used ecological niche models based on different assumptions. The first scenario assumes that species’ distributions are in equilibrium with the environmental conditions across their ranges today and uses climatic information solely from species’ current ranges to calibrate ecological niche models, which we call current‐based models. The second scenario relies on climatic information from species’ natural ranges, that is, simulating a scenario where species have never experienced heavy anthropogenic stressors (Faurby & Svenning, 2015). The latter considers that fundamental niches tend to be conserved over time (Peterson, 2011) so that occurrence records from different periods should provide additional information on species’ climatic tolerances from environments without contemporary counterparts (Faurby & Araújo, 2018; Lima‐Ribeiro et al., 2017; Martínez‐Freiría et al., 2016). We name this approach a natural‐based model, to adhere to the terminology originally proposed by the authors of the PHYLACINE 1.2.1 dataset (Faurby et al., 2018).

For each species, we sampled random points within the species’ current and natural ranges, proportionally to its range size (Table S2), and used ecological niche models to generate potential distribution maps. To do so, we used MaxEnt, a presence‐background method in which the species’ distribution is an unknown probability along with the full background points, that is, non‐negative values that add up to one (Elith et al., 2011). MaxEnt is robust to the presence of a moderate level of locational error and still provides useful predictions of species’ environmental preferences (Graham et al., 2008). The values of predictor variables at localities within natural and current ranges restrict the unknown distributions so that the average and variance values of environmental predictor should be, therefore, close to empirical values (Graham et al., 2008; Merow & Silander, 2014). However, the complexity of the fit to the observed values can be adjusted by transformations on the original predictor values (“feature classes”) (Muscarella et al., 2014). In this work, we compared two combinations of feature classes: (1) MXS: Maxent Simple (only linear and quadratic features); (2) MXD: Maxent Default (linear, quadratic hinge, product, and threshold features, based on MaxNet package) (S. Phillips, 2017; S. J. Phillips et al., 2017). The combination of feature classes exhibiting the highest accuracy metric (F_pb, calculation explained in Supporting Information section “Model accuracy assessment”) was then selected for final projections.

Because the performance of ecological niche models is affected by the spatial distribution of background points (Barbet‐Massin et al., 2012), we used a stepwise approach to select and partition our background data. First, we defined species‐specific background extents and built a preliminary BIOCLIM (Busby, 1991) habitat suitability model to constrain our background data to regions considered less suitable, that is, suitability <0.3 (Engler et al., 2004). Then, we partitioned our presence‐background data by the latitude and longitude lines that divide occurrence localities into blocks in a checkerboard‐like fashion. Optimal band sizes were considered those that (i) exhibited smaller spatial autocorrelation, based on Moran's I, and (ii) minimized the difference in the number of records between blocks (Roberts et al., 2017; Velazco et al., 2019). Blocks were then alternatively used for fitting and evaluating the model, and evaluation metrics were summarized across iterations (Andrade et al., 2020; Muscarella et al., 2014). Continuous suitability surfaces (values ranging from 0 to 1) were then created to reflect the relationship between species’ occurrences and their environment. Values closer to zero indicate lower predicted suitability, whereas values closer to 1 suggest high environmental suitability.

Because we used range maps as a source of environmental information, the outcome of our ecological niche models represents the environmental conditions most frequently observed across species’ known distributional limits. Therefore, our results should not be interpreted in terms of probability of occurrence per se, nor can be directly translated into any abundance‐related metric. Such broadly defined climate envelopes, however, are meant to provide an initial assessment of species climatic suitability at the continental scale and are useful to investigate macroecological relationships between biotas and environments (Sales et al., 2019, 2020a, 2020b). We stress that our results should not be explored at face value in conservation assessments at local spatial scales.

Climate information from current and natural ranges was, therefore, used as input in our ecological niche models and was compared in terms of future potential distribution areas, for each species (model details and parameterization are presented in the Supporting Information section “Modeling framework”). To do so, the ecological niches inferred from current‐based and natural‐based models were projected onto climate forecasts, dated to the year 2090 and based on the 6th Assessment Report of the International Panel of Climate Change (IPCC, 2021). Two representative concentration pathways were considered: an optimistic level of emission of greenhouse gases (ssp126), and a more extreme scenario (ssp585). While the first scenario represents the “best case” future from a sustainability perspective, with temperature increases of less than 2°C, the latter assumes “no climate policy” and predicts potential increases of almost 5°C in global mean temperature. Both forecasts, therefore, anticipate increases in temperature that are above the 1.5°C thresholds suggested to avoid the negative effects of climate change on global ecosystems and human well‐being (IPCC, 2018).

Within each climate change scenario, we chose three distinct climate models to encompass uncertainty on future projections—namely BCC‐CSM2‐MR, CanESM5, and MIROC6—using a stepwise procedure aimed at minimizing similarity (Sanderson et al., 2015). The resulting projections of those three climate forecasts were combined to create a single consensus map of the potential future distribution for each species at each climate change scenario. Only cells that were predicted as suitable by all climate models were included in the final consensus map. Models were evaluated using methods that do not rely on true absence, namely the Jaccard’s (1908) and the Sørensen's (1948) similarity indexes, in addition to the F_pb, a proxy of the F‐measure based on presence‐background data (Li & Guo, 2013), according to the equations in Table S3. To restrict our analysis to regions likely accessible to species, we further limited access to future suitable areas according to species‐specific dispersal constraints.

File List: 

1) Current_based_future_binary.rar maps of species' future potential distribution, calibrated with environmental information from species' current ranges (available as Current_PHYLACINE.rar)

2) Current_PHYLACINE.rar maps of species' current realized distribution, obtained from PHYLACINE dataset

3) Present_natural_future_binary.rar maps of species' future potential distribution, calibrated with environmental information from species' present natural ranges (available as Present_natural_PHYLACINE.rar)

4)Present_natural_PHYLACINE.rar maps of species' present natural distribution, obtained from PHYLACINE dataset

Relationship between files:  Present natural ranges depict the locations species might have occupied in the present, had they not been affected by human presence. Current ranges depict the locations species are known to occupy in the present.