Global biogeographical regions reveal a signal of past human impacts

www.ecography.org ECOGRAPHY Ecography Page 1 of 13 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Manuel J. Steinbauer Editor-in-Chief: Miguel Araújo Accepted 15 October 2023 doi: 10.1111/ecog.06762 2023 1–14 2023: e06762 © 2023 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos Ecologists have long documented that the world’s biota is spatially organised in regions with boundaries shaped by processes acting on geological and evolutionary timescales. Although growing evidence suggests that historical human impact has been key in how biodiversity is currently assembled, its role as a driver of the geographical organisation of biodiversity remains unclear. Using non-volant terrestrial mammals we set up a bioregionalization procedure focused on two datasets, one describing the current ranges of terrestrial mammals, and another describing their potential natural ranges in absence of historic anthropogenic land use. We then quantified the relative importance of anthropogenic land use (5000 and 2000 years ago, and present time) to predict the current and natural biogeographical regions across the Earth. In general, past and present human land use were important predictors of current bioregions but did not largely contribute to predict natural bioregions. Past anthropogenic land use seems to have left an imprint on the taxonomic differentiation of some of the largest biogeographical realms, whereas land use at present stands out as a driver of the taxonomic differences between medium-sized subregions, i.e. within and among continents. Our findings suggest that anthropogenic actions during the last millennia have had a farreaching effect on the spatial organisation of the Earth’s non-volant mammals. Keywords: bioregionalization, historic anthropogenic land use, terrestrial mammals Introduction How is the world’s biodiversity organized, and why do large-scale patterns of taxonomic diversity change through natural geographic regions? These questions have attracted the attention of naturalists since the early 19th century (von Humboldt 1806, Wallace 1876, Ricklefs 2004, Daru et al. 2017, Dinerstein et al. 2017). A key step in understanding the organisation of biodiversity is to identify the assemblage of regions based on their shared organisms (Wiens 2011). Wallace was among the first to propose that the world’s fauna is organized hierarchically in broad regions Global biogeographical regions reveal a signal of past human impacts Marta Rueda ✉1,2, Manuela González-Suárez 3 and Eloy Revilla2 1Departamento de Biología Vegetal y Ecología, Universidad de Sevilla, Sevilla, Spain 2Departamento de Biología de la Conservación, Estación Biológica de Doñana (EBD-CSIC), Sevilla, Spain 3Ecology and Evolutionary Biology, School of Biological Sciences, University of Reading, Reading, UK Correspondence: Marta Rueda ([email protected]) Research article 13 Page 2 of 13 shaped by geographic and climatic factors (Wallace 1876). About 150 years later, the development of multivariate analytical techniques together with an increase in the availability of global species distribution databases have led to revaluations of Wallace’s proposal (Kreft and Jetz 2010, Procheş and Ramdhani 2012, Holt et al. 2013, Rueda et al. 2013) and the improvement of our understanding about the extrinsic determinants driving the main dissimilarities among bioregions (Ficetola et al. 2017). Broadly, what we know is that the interplay of multiple factors has jointly contributed to the formation of bioregions, but their prominence varies across the globe (Riddle and Hafner 2010, Ficetola et al. 2017, Mazel et al. 2017). Processes acting deeper in the past, like plate tectonics, are the most important in explaining the separation between strongly divergent biogeographical realms. Factors representing present-day ecological barriers, like climatic heterogeneity, determine the separation between less dissimilar regions, i.e. within and among continents; while mountains – the third factor in importance – operate at all levels of biogeographical differentiation (Ficetola et al. 2017). The potential role of anthropogenic impact in shaping bioregions has been largely ignored because extant biogeographical regions have traditionally been assumed to reflect the natural organisation of biodiversity that resulted from ecological, historical, and evolutionary processes acting over millions of years (Lomolino et al. 2010). Yet, rising evidence showing that human-mediated species introductions are already affecting extant biogeographic regions (Capinha et al. 2015, Bernardo-Madrid et al. 2019) challenges this view. Furthermore, accumulating evidence shows that Quaternary human activities induced shifts in the plant and animal communities we see today (Lyons et al. 2016), and there is little doubt that modern humans have been a major driver in the extinction of large mammals during the late Pleistocene and early Holocene (Sandom et al. 2014, Smith et al. 2018). Moreover, archaeological, palaeoecological and genetic data suggest that cumulative human transformation of ecosystems over millennia has resulted in dramatic changes in composition, community structure, richness and genetic diversity of a diverse array of organisms across taxonomic groups (Ellis et al. 2013, 2021, Boivin 2016, Mottl et al. 2021). This body of evidence raises the question of whether historical anthropogenic pressures may have been large enough to leave an imprint detectable today on biogeographical assemblages globally. To answer this question, we set up a bioregionalization and analytical procedure focused on two datasets; one describing the current distribution ranges of global terrestrial mammals, and another set inferring their present natural ranges, which represents estimates of where species would potentially live without anthropogenic pressures (Faurby and Svenning 2015, Faurby et al. 2018). Using both types of distribution ranges we build hierarchical bioregionalizations and compute models to identify the variables (plate tectonics, climatic heterogeneity, mountains, and anthropogenic land use over the Holocene, specifically 5000 and 2000 years ago as well as present time) that best predict the assemblage of species in bioregions. We hypothesize that if anthropogenic land use has had a notable effect on the biogeographical assemblages and boundaries as we know them, then 1) the biogeographic configurations built using current distribution ranges will differ from those based on natural ranges; and 2) human land use should predict the biogeographic patterns obtained using current ranges, but not patterns obtained using natural ranges. Material and methods Generating presence–absence matrices using range distributions Current and natural distribution ranges of terrestrial mammals were obtained from PHYLACINE 1.2 (Faurby et al. 2018), which contains range maps for all 5831 known mammal species that lived since the last interglacial (~ 130 000 years ago until present). Ranges for current extant species contained in PHYLACINE are IUCN (2016) distribution maps, whereas present natural ranges represent estimates of where species would potentially live if they had never experienced anthropogenic pressures, and – in the case of extinct species – had not gone extinct (Faurby and Svenning 2015 for details). Note that Faurby and Svenning (2015) did not use present or past human land use variables to estimate PHYLACINE natural ranges, thus, natural ranges were not defined considering the predictors we test here avoiding circularity. Both current and natural ranges in PHYLACINE are projected to Behrmann cylindrical equal area rasters with a cell size of 96.5 by 96.5 km. Using these data, we generated two presence–absences matrices, one for the current and one for the natural ranges, in which every row represents a grid cell and every column a species. From these matrices we excluded 1) Homo species; 2) bats, because their high dispersal capacity and their almost global distribution make it difficult to establish clear biogeographical boundaries (Wallace 1876); 3) marine mammals, including polar bears, and other fully aquatic species (https://marinemammalscience.org/ science-and-publications/list-marine-mammal-species-sub- species for a full list of marine species not included); 4) cells containing less than 50% of land area to approximate equalsize samples; and 5) cells containing fewer than five species to reduce potential distortions caused by having few taxa (Kreft and Jetz 2010). These exclusion criteria rendered a total of 14 087 cells representing 3960 terrestrial mammal species for the current matrix, and 14 151 cells representing 4306 for the natural matrix. Given that differences between the biogeographical configurations of current and natural bioregions could be driven by the huge distribution ranges typical of extinct megafauna, we also built presence–absence matrices using only the extant current mammals (i.e. not including fossil ranges) that have had range contractions (i.e. have and estimated natural distribution range in PHYLACINE different from the IUCN current range). Overall, we obtained 536 extant mammals whose distribution ranges have changed. Note that a few extant species became extinct on some continents as part 16000587, 0, Downloaded from https://nsojournals.onlinelibrary.wiley.com/doi/10.1111/ecog.06762 by Universidad De Sevilla, Wiley Online Library on [14/02/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 3 of 13 of the megafauna extinction event but survived on others. Consequently, we adjusted distribution ranges for those continentally extirpated species. To do that, first we did a literature review (Barnosky et al. 2004, Koch and Barnosky 2006) to find those species that had become extinct prehistorically on some continents but survived on others (Eurasia was considered a single continent) (Supporting information for species selected). Then, we excluded any portion of the species’ range found on that continent. Building biogeographical regionalizations To each of the presence–absence matrices we applied a machine-learning algorithm referred to as affinity propagation (AP hereafter) (Frey and Dueck 2007) to build biogeographical regionalizations at different biogeographical resolution. We chose AP because it can compress massive datasets very efficiently (i.e. with lower error in comparison with other clustering methods, Frey and Dueck (2007) for details), and its good performance in hierarchical bioregionalization procedures (from few grid cells to large realms) has already been demonstrated in a previous work (Rueda et al. 2013). The AP algorithm works by detecting special data points called exemplars, and by a message-passing procedure it iteratively connects every data point to the exemplar that best represents it until an optimal set of exemplars and clusters emerges. In contrast to algorithms in which exemplars are found by randomly choosing an initial subset of data points, AP takes as input measures of ‘similarities’ between pairs of data points (grid cells here) and simultaneously considers all the points as potential exemplars. The optimal set of exemplars is the one for which the sum of similarities of each point to its exemplar is maximized. We first used the presence–absence matrices to calculate pairwise similarities between pairs of cells, and selected Hellinger distance as a similarity index. The Hellinger distance is a modification of the Euclidean distance (Legendre and Gallagher 2001) and is used to avoid the ‘double-zero’ problem, i.e. when two sites that have no species in common are assigned the same distance as two sites that share the same species; and the ‘species–abundance paradox’, which frequently occurs when two sites share only a small fraction of all the species in the same regional pool. This is expected to be a particular problem at the margins of biogeographical regions where sites may be quite different from one another rather than in the centre of a region where sites are likely to be very similar in their species assemblages (Legendre and Legendre 1998, Gagné and Proulx 2009). We then used a protocol based on a successive application of AP (Rueda et al. 2013) to obtain a biogeographical upscaling from the smallest possible bioregions (i.e. the highest biogeographical resolution) to the largest ones. We performed an initial AP analysis involving all grid cells of the similarity matrix. This first AP run generates the optimal solution of the highest resolution bioregions, while also identifying its exemplars. We obtained 1128 clusters and exemplars for the current distributions and 1053 for the natural distributions. Then, using the exemplars as the new units of analysis, we again conducted an AP, i.e. we calculated a new similarity matrix and re-run a new AP. This process was repeated until a small and coherent number of large clusters emerged; specifically, four runs were needed to get nine (with the current distributions) and eight (with the natural distributions) large bioregions (Supporting information). Finally, to obtain maps of each clustering result, we classified each grid cell (row) of every presence–absence matrix according to the cluster to which they were assigned in its corresponding AP analysis. AP analyses were performed using the ‘APCluster’ R package (www.r-project.org, Bodenhofer et al. 2011). Assessing differences in biogeographical configurations The degree of similarity between clusters (i.e. bioregions) of current and natural bioregions was estimated using the Jaccard index. For that, we used the cluster_similarity function of the ‘clusteval’ R package (www.r-project.org, Ramey 2012). This function computes the similarity between two clusterings of the same data. Jaccard index ranges between 0 (no similarity) to 1 (perfect match). Predictors To each of the cells of the presence–absence matrices we assigned a mean value of several predictors. We considered variables previously tested as determinants of biogeographical boundaries related to climatic heterogeneity, orographic barriers, tectonic movements, and instability of past climate (Ficetola et al. 2017) plus variables related to anthropogenic land use over the Holocene. In particular, we tested four climatic variables: annual total precipitation, mean annual temperature, seasonality in temperature, and seasonality in precipitation. All climatic variables were extracted from the WorldClim dataset (Fick and Hijmans 2017) up-scaled at a 96.5 km resolution. These variables represent both average conditions and variability within years and have been shown to be determinants of vertebrate distributions (Sexton et al. 2009). Climatic conditions have strongly shifted over the Quaternary period, and have been shown to play a role in the present-day species distributions, endemism, and assemblages (Araújo et al. 2011). To test for the potential effect of past climate change or stability, we calculated the average velocity of climate change since the Last Glacial Maximum (LGM; ~ 22 000 years ago) and since the Mid-Holocene (MH; ~ 6000 years ago) (Sandel et al. 2011). For that, we used the mean annual temperature and annual total precipitation for the MH and the LGM calculated by means of the model for interdisciplinary research on climate (MIROCESM) (Watanabe et al. 2011). Plate tectonics are responsible of the long-term isolation of the biotas, whereas mountain ranges represent major barriers to dispersal for most mammals (Lomolino et al. 2010, Ficetola et al. 2017). Plate tectonics were obtained from Bird (2003). Each grid cell was assigned the tectonic plate to which 16000587, 0, Downloaded from https://nsojournals.onlinelibrary.wiley.com/doi/10.1111/ecog.06762 by Universidad De Sevilla, Wiley Online Library on [14/02/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 4 of 13 it belongs. We removed those cells that were represented by more than one tectonic plate to avoid potential wrong ‘cooccurrence’ patterns, i.e. mixing faunistic elements separated by different plates even though the species do not actually co-occur. We used the GTOPO30 to calculate the difference between minimum and maximum elevation per grid cell to obtain a metric of elevational variance aimed to capture the complexity of topographic changes. Finally, past anthropogenic land use was obtained from Ellis et al. (2013). This dataset consists of two global spatially explicit reconstructions two main models of land use across the Holocene, HYDE and KK10 models. We chose the more realistic KK10 model, which assumes that humans use land more intensively when population density is high and land scarce (Kaplan et al. 2011). In counterpoint, the HYDE model omits land-use intensification and predicts that except for the developed regions of Europe, human use of land was insignificant in every biome and region before 1750 (start of the industrial revolution). The KK10 model dataset contains raster information for human land use for 10 different time periods; 8000, 5000, 3000, and 2000 years ago, and 1000, 1500, 1750, 1900, 1950, and 2000 of the present era. We first selected four different time spans: 8000, 5000, and 2000 years ago, and current (year 2000 of the present time), but the descriptor of land use 8000 years ago was not included in the random forest models due to its low variability at global level (Supporting information). Evaluating potential drivers of the biogeographical configurations We used random forest classification models (Breiman 2001) to assess the factors that may predict the classification of cells bioregions and to estimate the relative importance of the predictors. Note that random forest models were only carried out for large and medium-sized bioregions or subregions, as some of the predictors (e.g. plate tectonics) do not make ecological sense to explain smaller bioregions. Random forest is a machine learning method based on a combination of a large set of decision trees. Each tree is trained by selecting a random set of variables and a random sample from the training dataset (i.e. the calibration dataset). The accuracy of the models is given by the out-of-bag (OOB), an estimate of the misclassification rate that represents an unbiased error rate of the model that is calculated by counting how many cases in the training set are misclassified and dividing the number by the total number of observations. First, we split every dataset into training set (70% data, used to fit the model) and test set (30% data, used to evaluate the performance of the model). Second, with the training set we performed cross-validation to evaluate different hyper parameters of the model using grid random search (Supporting information). Random forests main hyper parameters include the number of decision trees in the forest (ntree) and the number of features considered by each tree when splitting a node (mtry). Values for ntree ranged between 7505 to 4604, while values for mtry ranged between 3 and 7 (Supporting information). Third, we run the Random Forest classification models using the training dataset and, finally, we estimated the accuracy of our models by checking the predictions against the actual values in the test dataset. Random forest models were computed using the R package ‘RandomForest’ (Liaw and Wiener 2002). Additionally, we also defined classification models using XGBoost (R package ‘xgboost’, www.r-project.org, Chen et al. 2023) to compare algorithm performance in terms of accuracy. Here, we report results from random forest as it had better overall accuracy. Random forests are able to disentangle interacting effects and identify non-linear and scale-dependent relationships that often occur at the scale of the analysis performed here among multiple correlated predictors (Cutler et al. 2007). Evidence from genomic studies suggests that variable importance measures may show a bias towards correlated predictor variables (Nicodemus et al. 2010), hence we eliminated those predictors showing a moderate to high Pearson’s correlation (r > 0.50, Supporting information). As a result, our final models included seasonality in temperature and precipitation as representatives of climate heterogeneity; velocity of climate change since the LGM and the MH to the present as representatives of past climate change; plate tectonics; variation in elevation representing mountains; and human land use at the present, 5000, and 2000 years ago. Human land use 5000 and 2000 years ago showed a large correlation (r = 0.87), so we decided to perform separate models for each time period. Globally the correlation between human land use at present and 2000 years ago is relatively low (r = 0.29); however, we also decided to perform independent models for both predictors because Pearson’s coefficient using modified t-test (below) showed that the correlation between both variables is influenced by spatial autocorrelation. Since our models include spatial data, the selection of predictors with random forest could be partly influenced by spatial autocorrelation. To assess this point, we performed correlations between predictors using the method of t-test modification (Clifford et al. 1989, Dutilleul 1993) to incorporate spatial structure in the Pearson’s correlations. This is a hypothesis testing procedure based on a modification version of the correlation coefficient that takes into account the spatial structure of the predictors. This analysis was computed using the function modified.ttest of the R ‘SpatialPack’ package (www.r-project.org, Vallejos et al. 2020). The results showed evidence of spatial autocorrelation in many of the correlations, even for low Pearson’s correlation values (Supporting information). Therefore, to account for spatial autocorrelation we also run random forest classification models incorporating spatial structure as a possible predictor. To do that we first generated Morans’s Eigenvector Maps (MEMs) (Dray et al. 2006), and ranked and reordered them based on their Moran’s I. These eigenvectors can be included in model as spatial predictors. Here we include the first three eigenvectors. MEMs were calculated using the mem function from the ‘spatialRF’ package in R (www.r-project.org, Benito 2021). The algorithm used in the process of regionalization only uses species presence–absences, which means that 16000587, 0, Downloaded from https://nsojournals.onlinelibrary.wiley.com/doi/10.1111/ecog.06762 by Universidad De Sevilla, Wiley Online Library on [14/02/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 5 of 13 humans can only have a role if present and if causing local extinctions, thereby changing the amount of shared species between grid cells. At the largest biogeographical scale two bioregions, Madagascar and Australia, appear characterised by an extremely low historical human land use (Supporting information). Note that in both cases human pressures may be poorly captured by the input layers of the KK10 model used to define the variables of past anthropogenic land use, as the role of pre-historic human impact in the extinction of the megafauna of both Australia and Madagascar is well known (Miller et al. 2005, Hansford et al. 2018). In any case, the ‘false’ absence of human impact could incorrectly give support for historical land use as a cause. Besides, both regions have a 100% endemicity for non-flying mammals, meaning that humans cannot change their clustering without introducing species. Therefore, to assess to what extent these factors might be affecting the results, we also ran the models excluding Australia and Madagascar. Finally, we note that our models are unbalanced regarding the number of climate versus human land use variables, which overemphasize the importance of human impact. To evaluate this possibility, we also fitted models without climate variables. Calculating global and local predictor importance We measured variable importance using the mean decrease in accuracy, which is obtained by permuting randomly each variable and assessing the decrease in classification accuracy of the model (Liaw and Wiener 2002). Thus, the greater the loss of accuracy due to the permutation (or exclusion) of a single variable, the more important that variable is considered, i.e. variables with a large mean decrease in accuracy are more important for classification of the data. Random forest also calculates the local variable importance, which defines the importance of a variable in the classification of a single sample (grid cell here) and therefore shows a direct link between variables and samples (Touw et al. 2012). The local importance score is derived from all trees for which the sample was not used to train the tree. For each sample, the percentage of votes for the correct class in the variable permuted OOB data is subtracted from the percentage of votes for the correct sample in the original OOB data (i.e. the untouched OBB data) to assign a local importance score. The score reflects the impact of a variable on correct classification of a given sample: negative (i.e. the variable does not positively contribute to the classification), 0 (the contribution is neutral) and positive (the variable positively contributes to the classification). Given that local importance values are noisier than global importance ones we run the same classification five times (five per biogeographical scale) and averaged the local importance scores to obtain a robust estimation of local importance values (Touw et al. 2012). We present maps of mean local importance values for present and past human land use. In these, it is important to note that high local importance can reflect both areas where human land use was high and areas where it was negligible, i.e. both are equally valuable for the classification analysis (Supporting information). To facilitate interpretation, in the main text we present maps of local importance only for areas where land used was > 10%. All analyses were performed in the R ver. 4.3.0 (www.rproject.org). Results Bioregionalizations Both for the current and natural bioregionalizations we identified a hierarchical system of biogeographic regions with four levels (Fig. 1). For current mammalian ranges, the broadest delineation of nine large bioregions was mostly consistent with the original maps of Wallace’s realms and subregions (Fig. 1a). For natural mammalian ranges, the eight bioregions obtained showed a different biogeographic arrangement with respect to the nine large current ones, especially in the Northern Hemisphere (Supporting information). Thus, while the natural broadest bioregions show a single Holarctic region, in the current ones this is separated into a Nearctic and a Palearctic region. Differences between current and natural bioregions extended to the subregions or medium-sized bioregions and smaller bioregions (Supporting information) and were confirmed for the bioregionalizations carried out with the 536 species with different natural and current ranges (Supporting information). Figure 1. Current and natural broader regions or realms. Only results for the broader regions are shown here: nine current bioregions (a), and eight natural bioregions (b). The greatest differences between both biogeographical configurations are more relevant in the Northern Hemisphere, where, for the natural bioregions, the Palearctic and the Nearctic are united forming a great extension, the typical biographical boundary between the Nearctic and the Neotropical realms disappears and extends towards the north, and India is no longer part of the Asian region and is included instead within the ‘Saharo-Arabian’ region. 16000587, 0, Downloaded from https://nsojournals.onlinelibrary.wiley.com/doi/10.1111/ecog.06762 by Universidad De Sevilla, Wiley Online Library on [14/02/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 6 of 13 Global importance of predictors of the biogeographical patterns for regions and subregions For the broader bioregions, human land use at present had a moderate relevance in the classification models (Fig. 2a). In contrast, past human land use (5000 and 2000 years ago) was the third most important predictor of the biogeographic structure of current realms behind plate tectonics and seasonality in temperature (Fig. 2c, e), with a moderate, non-negligible predictive value (~ 22 and 21% of cells classified correctly when the model was run just with human land use 2000 and 5000 years ago, respectively, Fig. 3). In broad accordance with our prediction, past and present human land use were poorer predictors of natural bioregions (Fig. 2b, d, f). Results were consistent also if we excluded Australia and Madagascar (Supporting information). The misclassification rate (out-of-bag, OOB hereafter) for these set of models was quite low, never exceeding 3.6%, and the model accuracy ranged around 0.97. For medium-sized bioregions or subregions, past and present human land use were the second predictors in importance for the current bioregions (Fig. 4a, c, d; OOB values in general < 8.6%, model accuracy 0.92–0.94). As hypothesized, the importance of land use at present was smaller for the natural bioregions; although, unexpectedly, past land use remained a strong driver (Fig. 4b, d, e). Results for models accounting for the spatial autocorrelation (Supporting information) showed a similar tendency, especially for the broader bioregions and past human land use (OOB values in general < 2.5%, model accuracy 0.98– 0.99). In fact, for these models the differences for past human impact are sharpened when current broader bioregions are compared with the natural ones. On the contrary, the importance of human land use at present fades in favor of other predictors when the spatial structure is taken into account for both current and natural broader bioregions. For mediumsized bioregions, the importance of plate tectonics decreases and it is no longer the main predictor (OOB values in general < 3.8%, model accuracy 0.96–0.97). When accounting for spatial autocorrelation the importance of human land use at present becomes greater for current bioregions and decreases for natural bioregions. However, the past human land use has a completely opposite effect to that expected since it becomes the main predictor of medium-sized natural bioregions, implying that spurious factors might be acting on this scale. Models without climatic variables had lower accuracy as expected but did not alter our main findings with similar variable importance rankings that support the role of human land use (Supporting information). Local importance of predictors of the biogeographical patterns Maps showing local importance scores revealed differences in the location of grid cells where past and present human land use have been analytically important for the classification of the broader current bioregions (Fig. 5). The cells where human land use at present was most relevant are located in southern Africa, the easternmost boundary between the Palearctic and the Asian regions, and part of the northern boundary between the Palearctic and ‘Saharo-Arabian’ regions (Fig. 5a). The cells where human land use 2000 years ago was more important are located mainly throughout Central America and the Andes mountain range, Europe except for the northernmost part, most of the Asian region, and the boundary between the Palearctic and the Asian regions (Fig. 5b). The cells where human land use 5000 years ago was more important mostly coincided with those of 2000 years ago, but with minor scores of importance (Fig. 5c). Some of these areas coincided for the broader natural bioregions, yet the importance of Central America and the Andes for the classification of natural bioregions weakened (Fig. 5e–f) whereas the importance of the area of southern India increased (Fig. 5d–f). For the medium current bioregions, the grid cells where the human land use at present was most relevant for the classification of bioregions included those highlighted for the broader bioregions plus most part of Europe, southern India, eastern Australia, and a large area of the Nearctic (Fig. 6a). When focusing on the human land use 2000 years ago, the areas affected by human land use were similar to those of the broader bioregions plus a Sub-Saharan area crossing Africa from west to east (Fig. 6b). Again, the cells where human land use 5000 years ago was more important mostly coincided with those of 2000 years ago, but it is noteworthy that the importance of the Sub-Saharan area disappears (Fig. 6c). The location of cells affected by human land use for the medium natural bioregions quite coincided with that for current bioregions but with less importance scores generally in Central America and the Sub-Saharan area (Fig. 6e–f). Discussion Our study shows that human activity during the last few millennia has left its mark on the spatial organisation of the Earth’s biodiversity. First, we show that biogeographical configurations based on current and natural mammal distributions differ, particularly in the Northern Hemisphere. These differences were not solely the result of including extinct megafauna – with huge distribution ranges – in the natural biogeographic regions, as biogeographical patterns for the 536 extant terrestrial mammals with different current and natural ranges confirmed our findings. Second, we find that human land use is an important predictor of current bioregions but does not largely contribute to predict natural bioregions. For example, for the broader regions or realms, past human land use was the third most important predictor behind plate tectonics and temperature seasonality, but was of minor importance on natural regions. The difference between current and natural bioregions for past human land use was even more exacerbated when spatial autocorrelation was considered, such that past human land use became the second 16000587, 0, Downloaded from https://nsojournals.onlinelibrary.wiley.com/doi/10.1111/ecog.06762 by Universidad De Sevilla, Wiley Online Library on [14/02/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 7 of 13 most important factor after plate tectonics in predicting the configuration of the current broader bioregions. This implies that past anthropogenic actions have seemingly left a perceptible biogeographical imprint on the assembly of the global realms we recognize today. Moreover, for the medium-sized subregions, human land use at present was among the first most important predictors, suggesting that recent human activity is already causing changes in the biogeographical assemblages at the subregional scale, i.e. within and among continents. Again, the importance of human land use at Figure 2. Ranking of importance values for the drivers of taxonomic differentiation for the current and natural broader bioregions. Above panels show importance values from models with human land use at present for (a) current and (b) natural bioregions, while the panels below show importance values from models with human land use 2000 years ago for (c) current and (d) natural bioregions, and human land use 5000 years ago for (e) current and (f) natural bioregions. The numerical values associated with each horizontal bar are mean decreased in accuracy values (i.e. how much the accuracy decreases when the variable is excluded from the model) and represents variable importance. Importance was measured by the drop-in classification accuracy after predictor randomization in random forests. Consequently, higher values of mean decreased in accuracy indicate variables that are more important to the classification. Values from 0 to 1 in the x-axis represents relative importance values (i.e. division of each value of mean decrease in accuracy by the largest value) and show a ranking of variable importance. OOB (out-of-bag) represents the percentage of cells misclassified in each model. VCC = velocity of climate change; MH = Mid-Holocene; LGM = last glacial maximum; annual P = annual precipitation; mean T = mean annual temperature; ya = years ago; Var. in elevation = variation in elevation. 16000587, 0, Downloaded from https://nsojournals.onlinelibrary.wiley.com/doi/10.1111/ecog.06762 by Universidad De Sevilla, Wiley Online Library on [14/02/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 8 of 13 present increased as a predictor of current medium-sized subregions when the spatial autocorrelation was considered. The mammalian distributions that we know today are not only a reflection of the most recent human actions but also those exerted during the last few millennia (Faurby and Svenning 2015, Polaina et al. 2020). Recent studies had already shown that human actions in the present, in terms of human-mediated species introductions, can be evident at biogeographical scales (Capinha et al. 2015, Bernardo- Madrid et al. 2019), but here we also show evidence of an effect of past human land use. What has happened from 5000 years ago that has left a noticeable footprint in the Earth's biogeographical regions? This period coincides with the emergence of an ‘intense’ agricultural land expansion and the settlement of densely populated areas in the Mediterranean, southwest Asia, south Asia, and eastern China (Kirch 2005, Ellis et al. 2013, Stephens et al. 2019). Long-term impacts to enhance agricultural productivity, including forest clearing, increased fire frequencies, species invasions, and soil erosion, were already apparent in some regions at this time (Ellis et al. 2013, Boivin et al. 2016), and may have transformed vegetation structure and species composition across many regions, from the Mediterranean to the Tropics (Stephens et al. 2019, Zheng et al. 2021). In fact, reconstructions suggest that > 20% of Earth’s temperate woodlands had been already impacted by humans by the late Holocene (Ellis et al. 2013). Indeed, our results on local importance show the match between some of the biogeographical boundaries and regions of the New World and the distribution – and presumed impact – of ancient Mesoamerican and South American pre- European civilizations (Richard 1997, Dunning et al. 2012). In the Old World, we note the prominence of past human land use covering most of the Asian realm territory, associated to the rise of rice cultivation in eastern, southern, and southeast Asia over the Holocene (Silva et al. 2015, Zheng et al. 2021). Actually, cumulative archaeological data show that the scale of agriculture and land use in these regions was significant; expansion of land area used for livestock and rice Oryza sativa paddy agriculture in southern Asia was sufficient to increase atmospheric methane emissions between 4000 and 1000 years ago (Fuller et al. 2011). Moreover, the atmospheric CO2 decline registered between 1570 and 1620 has been attributed to the uptake by vegetation following the agriculture abandonment caused by the population crash after the arrival of Europeans in America (Lewis and Maslin 2015). Past human land use can also partly explain differences between the southern Palearctic boundary and the Afrotropic below the Sahara. Both can be related to the transformation of European and Near Eastern landscapes during the Roman period (Butzer 2005), and the expansion of sorghum Sorghum bicolor cultivation 3000 to 2000 years ago in Sub-Saharan Africa (Boivin et al. 2016). While we found differences in importance of human land use 5000 and 2000 years ago between current and natural medium-sized bioregions, we were surprised that past human land use was an important determinant of natural bioregions, which were based on distributions of mammals in the absence of human impact. These results could reflect an overlap between human land use and biogeography. If the distribution of humans during the middle and late Holocene reacted to the same drivers as the distributions of terrestrial mammals, human land use could appear as a predictor even if it had no causal effect. In fact, the importance of this variable appears to be a consequence of how human land use was distributed globally 5000 and 2000 years ago, or more accurately where human activity was largely absent. Note that while current human land use is distributed over wide areas of the planet, that of 5000 and 2000 years ago was centralized in specific areas coinciding with the distribution of ancient civilizations. This means that, on a subregional scale, many bioregions can be well predicted by the absence of human land use. Indeed, an additional analysis that excluded subregions with no reported human land use 2000 years ago showed very little importance of past human land Figure 3. Diagrams describing the accuracy in classification of the three most important drivers of the current broader bioregions for models including human land use 2000 and 5000 years ago. Values indicate the percentage of cells correctly classified (1-OOB error) in their bioregion, and are obtained from running random forest classification models for each individual variable, by pairs and for the three most important factors together. Note that the number of variables modifies the parameters of the models and causes there to be a variation in the out-of-bag (OOB) when compared with models containing all variables. 16000587, 0, Downloaded from https://nsojournals.onlinelibrary.wiley.com/doi/10.1111/ecog.06762 by Universidad De Sevilla, Wiley Online Library on [14/02/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Page 9 of 13 use for both current and natural medium-sized bioregions (Supporting information). While our results show perceptible impacts of past human land use on current biogeographical patterns for terrestrial mammals, there are some limitations in the available data that could potentially bias our results. First, IUCN distribution ranges do not represent true areas of habitat, meaning that there may be a fair amount of unsuitable habitat in the ranges resulting in errors of commission (Brooks et al. 2019). Yet, our approach (assigning presences to grid cells) does not assume a given species is found in all areas of a cell, but rather that it could be found in some areas, which for cells of the size considered here (96.5 × 96.5 km) should greatly reduce commission errors. Second, IUCN distribution maps are the most comprehensive global source of range data for terrestrial mammals; however, we are aware that sampling effort, that influences range definition, varies across the world. Unfortunately, that spatial bias is also present in other available sources, like GBIF, in which nation-wide differences in funding and data sharing have led to huge differences in data contributions (Yesson et al. 2007, Beck et al. 2014). Third, while in PHYLACINE (Faurby et al. 2018) it is claimed Figure 4. Ranking of importance values for the drivers of taxonomic differentiation for the current and natural medium-sized bioregions or subregions. Above panels show importance values from models with human land use at present for (a) current and (b) natural bioregions, while the panels below show importance values from models with human land use 2000 years ago for (c) current and (d) natural bioregions, and human land use 5000 years ago for (e) current and (f) natural bioregions. Specifications as in Fig. 2. 16000587, 0, Downloaded from https://nsojournals.onlinelibrary.wiley.com/doi/10.1111/ecog.06762 by Universidad De Sevilla, Wiley Online Library on [14/02/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License