Studying past community dynamics can provide valuable insights for anticipating future changes in the world's biota. However, the existing fossil record is too sparse to enable continuous temporal reconstructions of wholesale community dynamics. In this study, we utilise machine learning to reconstruct Late Quaternary community structure, leveraging the climate-trophic structure relationship. We followed a four-stage approach: 1) identify and map trophic structure units (TSUs) at the global scale based on the guild richness and composition of extinct and extant terrestrial mammal species weighing over 3 kg; 2) train a random forest model to predict TSUs based on contemporary climatic conditions; 3) hindcast the global distribution of TSUs using climatic conditions as reconstructed over the past 120,000 years; and 4) compare TSU hindcasts against elements of community trophic structure as estimated with the fossil record. Models project significant shifts in the geographical distribution of community trophic structures, with more pronounced changes occurring during the Pleistocene-Holocene transition. These shifts exhibit regional variations, particularly in Eurasia and North America, where the models project notable reductions in the distribution of less complex trophic structures over the last 24,000 years. Hindcasts partially identified the alterations in community structure seen in the fossil record, demonstrating a match between the observed and predicted times of change in mammal community structure (between 24 and 8ka BP). However, climate-driven projections of trophic guilds diverged from fossil records during the Holocene. While the fossil record indicated a decrease in the number of grazers and carnivores, our models projected an increase in these numbers. Characterising community-wide responses to climatic changes is essential to address key questions about past and future impacts of such changes. Although further research is needed to refine the models, our approach offers a perspective for addressing the complex interactions among climate and trophic community organisation across ecosystems.

Our approach consists of four stages: 1) identify and map in a 1 × 1° grid cell surface the community TSUs that emerge from the raw data on species distributions, with species reclassified into guilds so that patterns of co-occurrence among guilds are examined instead of species co-occurrences; 2) train a random forests model to predict the community TSUs in every grid cell given contemporary climatic conditions; 3) use the trained model to hindcast the global distribution of community TSUs based on past climate data; and 4) assess the usefulness of hindcasts by examining the match between hindcasted community TSU transitions and observed community transitions in the fossil record. 

Identification of community TSUs

Following Mendoza and Araújo (2022), we characterized the trophic structure of mammal communities in every 1 × 1° grid cell based on the number of species from each trophic guild that occurred therein. This resulted in a data matrix (with n trophic guilds as columns and m grid cells as rows), where each cell represents a point in an n-dimensional ‘trophic space’ defined by the number of species from each trophic guild. We then performed a fuzzy-clustering algorithm on this trophic space to identify those well-defined clusters bringing together cells with a similar pattern in their trophic structure. This search of repeated patterns conforms to the expectation that communities, even if continuously distributed in a geographical space, will not be continuously distributed in the trophic space, rather forming well-defined clusters (here referred to as trophic structure units – TSUs) (Mendoza and Araújo, 2019, 2022).  

The optimal number of clusters was determined by tracking changes in the AMDi through a series of clustering analysis runs with increasing numbers of user-defined clusters (from 2 up to 15 clusters). In fuzzy clustering, clusters with elements showing a high degree of membership are considered well-defined clusters, whereas clusters with many elements with unclear membership are considered diffuse or fuzzy. The AMDi measures the average membership degree of the samples making up the clusters. It ranges from 0 to 1, being always zero when the artificial samples are randomly distributed and close to 1 when clusters all the samples have a high degree of membership for a cluster. Because of the many samples and dimensions in our data, there are many possible clustering solutions, so we selected the solution with the highest AMDi after 100 replicates. Please refer to Mendoza and Araújo (2022) for further details on the conceptual and methodological aspects of the AMDi approach.

We identified and mapped the community TSUs using the global distributional ranges of non-marine mammals obtained from the IUCN (IUCN, 2021) and Phylacine v.1.2.1 (Faurby et al., 2018, 2020) databases. For both, the AMDi approach showed a peak for 6 user-defined clusters, indicating the occurrence of six community TSUs (see the supplementary material). We used the two datasets to assess the model’s sensitivity to baselines with and without human influence, and to data availability. While one model makes a hindcast of community TSUs based on the contemporary distribution of mammals (including reduced distributional ranges and excluding extinct species, hereafter ‘Mendoza and Araújo dataset’), the other makes it based on a distribution recreating a world without recent human-induced extinctions (hereafter ‘Phylacine natural dataset’). Besides, the latter includes a lower number of species and trophic guilds (Mendoza and Araújo: 9 trophic guilds of 5272 terrestrial mammal species vs. Phylacine natural: 7 trophic guilds of 1315 terrestrial mammal species with body mass above 1kg) to recreate the possible influence of the incomplete fossil preservation of small mammals and the inability to infer some trophic guilds from fossil evidence (i.e., pollinators). 

Figure 1 shows the distribution of the six TSUs and their main characteristics using the Phylacine Natural database. From this point onwards, we focus on the Phylacine-based model since it provides greater comparability with the available fossil record. This database contains the pre-human ranges for extant and extinct mammal species derived by combining the IUCN-based distributions of species unaffected by human pressure with historical distributions of species, fossil co-occurrence, and known range modifications caused by humans. Furthermore, it includes extinct mammals and trophic guild categories that are inferable from taphonomic evidence (i.e., Browser, Grazer, Mixer-feeder [Herbivore], Carnivore, Mesocarnivore, Omnivore, or Invertivore). Results using the second baseline (Mendoza and Araújo dataset) are included as Supplementary material.

Model training and evaluation

The relationship between the occurrence of community TSUs (in both baselines) and bioclimatic variables plus net primary productivity (NPP) was modelled using a ‘Random Forests classifier.’ This is a machine learning algorithm widely used in classification problems given its extraordinary capacity to deal with highly correlated variables and fully exploit the information contained in them (Boulesteix et al., 2012). Random Forests (RF) grows many classification trees using different subsets of samples and predictors. To classify a new sample, each tree gives a classification, and the RF model chooses the classification having the most votes over all the trees in the forest (Breiman, 2001).

We used the default options of the RandomForest package (Liaw & Wiener, 2002) to train RF classifiers (500 trees, 4 variables at each split). However, we did not evaluate the model’s performance using the out-of-bag score provided by the package. When data are spatially autocorrelated, there is a risk that measures of model evaluation are inflated (Segurado et al., 2006). Therefore, to limit the effect of spatial autocorrelation on the evaluation of predictions, we performed a 36-fold cross-validation with geographically independent test data. The procedure consisted of splitting the globe into 36 bins of 10º longitude (from 180º West longitude to 180º East longitude) and fitting models that included the cells from all bins except the one used for testing (see also Benito, 2021). At the end of this sequential process, the Cohen’s Kappa coefficient (Carletta, 1996) was estimated across the 36 partial models to estimate the overall model performance. 

The (pre-industrial) modern-era data for 17 bioclimatic variables and NPP were extracted from Beyer et al.'s dataset (Beyer et al., 2020), which is a 1° resolution dataset covering climate and NPP for the last 120ka at a temporal resolution of 1-2ka (1ka time steps between the modern era and 22,000 BP, and 2ka time steps between 20,000 BP and 120,000 BP). This climatic reconstruction was performed by combining HadCM3 climate simulations across the last 120,000 years with high-resolution HadAM3H simulations across the last 21,000 years and modern-era instrumental data. Complementary, the Biome4 global vegetation model (Kaplan, 2003) was used to reconstruct the net primary productivity across the same time intervals. Bioclimate variables from this dataset are equivalent to the 19 WorldClim (Fick & Hijmans, 2017) used in macroecological analyses, but excluding BIO2 (Mean Diurnal Range) and BIO3 (Isothermality) which are not available. 

Hindcasting

Trained models were used to hindcast the global distribution of TSUs in a 1º × 1º grid cell surface at time intervals of 1ka during the Holocene (0-14ka BP) and 2ka during the Late Pleistocene (14-120ka BP). Bioclimatic variables and NPP for these 68 time-intervals were extracted from Beyer et al.'s dataset. 

The hindcasting produced 68 maps describing the global distributions of the six community TSUs over the last 120,000 years. To assess the global redistribution of trophic structures, we evaluate the temporal changes in the number and position of grid cells occupied by each TSU. However, as the RF produces classifications based on the number of trees (= votes), assessing changes based on category labels alone may be inaccurate. This is particularly true when the most probable community TSU in a grid cell is determined by a difference of few votes (i.e., TSU1 = 4500 vs. TS2 = 4400). Thus, to complement the categorical-based analysis, we employed a symmetric Procrustes analysis to evaluate how similar are TSU distributions at two intervals using the number of votes. This analysis searches for common structures between two datasets describing the same objects and measures how similar they are (Legendre & Legendre, 2013) (ranging from 0 to 1, where a value closer to 1 indicates that both structures are alike). In our case, the Procrustes analysis made a cell-by-cell comparison between present-day (baselines) and the past distributions in terms of the number of votes obtained by each TSU. As such, we created temporal series of 67 values describing the temporal decay in distribution similarity from the present to 120ka in the past.

To identify the time intervals with abrupt changes in temporal trajectories, we employed the function 'breakpoints' of the strucchange package (Zeileis et al., 2002). This function uses dynamic programming to find breakpoints that minimize a linear model's residual sum of squares (RSS) with m+1 segments; and the Bayesian Information Criterion (BIC) to find the optimal model as a trade-off between RSS and the number of segments (Zeileis et al., 2003). In turn, it can detect changes in the direction or magnitude in time series regression models by estimating the minimum breakpoints – or coefficient shifts – occurring along a time series based on piecewise linear models. Here, breakpoints in temporal trajectories allow identifying periods associated with sharp transitions in the relative frequency of each community TSU (when using the number of grid cells) or the similarity to the present day (when using the Procrustean correlation coefficient).

Model usefulness

The usefulness of hindcasts should ideally be assessed based on a comparison between projected trends with observed ones (Araújo et al., 2019). However, such independent model evaluation can be difficult owing to the sparse and incomplete fossil record (Fraser et al., 2021). Besides, there is a mismatch between the resolution of global hindcasts (1º) and the spatial coverage of fossil sites (frequently less than 10km²). As a first step towards independent temporal testing of the usefulness of hindcasts, we used the available dated fossil record to compare: (1) trends in similarity decay, (2) periods of significant community reshuffling, and (3) changes in the number of species of each trophic guild. Notice we use the term “usefulness” rather than “accuracy” because the quality of the fossil data is so limited that, whatever the result of such modelling-data comparison, it should be interpreted as an exploratory rather than a confirmatory analysis.

The first two comparisons were performed using temporal series of the hindcasted and observed decay in similarity. We obtained the first time series using the Procrustes analysis (as described below) on the grid cells covering the four Eurasian regions and the Australian continent. For the latter, we computed the similarity between contemporary (around 0.5ka BP) and past regional communities using the Bray-Curtis index, which estimates similarity based on the number of species per trophic guild that the two communities share. Both measures can be considered ‘equivalent’ in the sense that both quantify how similar were regional paleo-communities in comparison to current ones. To determine whether the hindcast periods of community reshuffling resemble the ones observed in the fossil record, we applied the ‘breakpoints’ function and compare the years (and confident intervals) in which sharp transitions are detected in the time series.

The third comparison to evaluate model usefulness was done between observed and projected changes in trophic guild richness. Projected changes in the number of species per trophic guild were computed as weighted averages of the number of species per guild that characterize each TSU, weighting by the number of grid cells that each TSU was projected to occupy at a given regional and temporal bin. For example, to compute the regional richness of each trophic guild in the Mediterranean around 14ka BP, we assess the type and number of grid cells occupied by each TSU (= the weighting factor). Then we estimate the average number of species per trophic guild that characterizes every type of community TSU in the present day (= trophic profiles; Figures 1 and SM2-3). Finally, we computed the weighted average of these average numbers using the number of grid cells as the weighting factor.