Ecosystem productivity affected the spatiotemporal disappearance of Neanderthals in Iberia

Fauna, culture and chronology datasets

A geo-referenced dataset of chronometric dates covering the late MIS 3 (55–30 kyr cal bp) was compiled from the literature (dataset 1). The dataset included 363 radiocarbon, thermoluminescence, optically stimulated luminescence and uranium series dates obtained from 62 archaeological sites and seven palaeontological sites. These chronological determinations were obtained from ten palaeontological levels and 138 archaeological levels. The archaeological levels were culturally attributed to the Mousterian (n = 75), Châtelperronian (n = 6) and Aurignacian (n = 57) technocomplexes. A number of issues can potentially hamper the chronological assessment of Palaeolithic technocomplexes from radiocarbon dates, such as pretreatment protocols that do not remove sufficient contaminants or the quality of the bone collagen extracted. Moreover, discrepancies in cultural attributions or stratigraphic inconsistencies are commonly detected in Palaeolithic archaeology. Information regarding the quality of date determinations and cultural attribution or stratigraphic issues is provided in the Supplementary Information.

Our dataset also included the presence of herbivore species recovered from each archaeo-palaeontological site (hereafter referred to as local faunal assemblages (LFAs)), their body masses and their chronology. The mean body mass of both sexes, for each species, was obtained from the PHYLACINE database⁵³ and used in the macroecological modelling approach described below (see ‘Carrying capacity of herbivores’). For visual representation purposes, the herbivore species were grouped into four weight categories: small (<10 kg), medium (10–100 kg), medium–large (100–500 kg) and large (>500 kg). The chronology of the occurrence of each herbivore species was assumed to be the same as the dated archaeo-palaeontological layer where the species remains were recovered. Thus, to estimate the chronological range of each species in each region, all radiocarbon determinations were calibrated with the IntCal20 calibration curve⁵⁴ and OxCAL4.2 software⁵⁵. The BAMs were run to compute the upper and lower chronological boundaries at a CI of 95.4% of each LFA (see ‘Chronological assessment’ for more details). One of the purposes of the current study was to estimate the potential fluctuations in herbivore biomass during the stadial and interstadial periods of the late MIS 3. Accordingly, the time spans of the LFAs were classified into the discrete GS and GI phases provided by Rasmussen et al.⁵¹.

Geographic settings

The Iberian Peninsula locates at the southwestern edge of Europe (Fig. 1). It constitutes a large geographic area that exhibits a remarkable diversity of ecosystems, climates and landscapes. Both now and in the past, altitudinal, latitudinal and oceanic gradients affected the conformation of two biogeographical macroregions with different flora and fauna species pools: the Eurosiberian and Mediterranean regions^13,46. In the north, along the Pyrenees and Cantabrian strip, the Eurosiberian region is characterized by oceanic influence and mild temperatures in the present day, whereas the Mediterranean region features drier summers and milder winters (Fig. 1). Between the Eurosiberian and Mediterranean regions, there is a transitional area termed Submediterranean or Supramediterranean. Lastly, the Mediterranean region is divided into two distinctive bioclimatic belts: (1) the Thermomediterranean region, located at lower latitudes, with high evapotranspiration rates and affected by its proximity to the coast; and (2) the Mesomediterranean region, with lower temperatures and wetter conditions (Fig. 1).

Previous studies have shown that zoocoenosis and phytocenosis differed between these macroregions in the Pleistocene^13,46. However, flora and fauna distributions changed during the stadial–interstadial cycles in the Iberian Peninsula, which suggests potential alterations in the boundaries of these biogeographical regions. The modelling approach used in this study to estimate the biomass of primary consumers is dependent on the reconstructed NPP and the herbivore guild structure in each biogeographical region. To test the suitability of the present-day biogeographical demarcations of the Iberian Peninsula during MIS 3, we assessed whether the temporal trends of NPP and the composition of each herbivore palaeocommunity differed between these biogeographical regions during the MUPT.

Chouakria and Nagabhusan⁵⁶ proposed a dissimilarity index to compare time series data by taking into consideration the proximity of values and the temporal correlation of the time series:

where S₁ and S₂ are the time series of data, u and v represent the values of S₁ and S₂, respectively, and p is the length of values of each time series. CORT(S₁, S₂) belongs to the interval (−1,1). The value CORT(S₁, S₂) = 1 indicates that in any observed period (t_i, t_i+1), the values of the sequence S₁ and those of S₂ increase or decrease at the same rate, whereas CORT = −1 indicates that when S₁ increases, S₂ decreases or vice versa. Lastly, CORT(S₁, S₂) = 0 indicates that the observed trends in S₁ are independent of those observed in S₂. To complement this approach by considering not only the temporal correlation between each pair of time series but also the proximity between the raw values, these authors proposed an adaptive tuning function defined as follows:

where

In this study, k was 2, meaning that the behaviour contribution was 76% and the contribution of the proximity between values was 24%⁵⁷. Hence, f(x) modulates a conventional pairwise raw data distance (d(S₁,S₂)) according to the observed temporal correlation⁵⁶. Consequently, dCORT adjusts the degree of similarity between each pair of observations according to the temporal correlation and the proximity between values. This function was used to compare the reconstructed NPP between biogeographical regions during MIS 3 in the Iberian Peninsula. However, two different biogeographical regions could have experienced similar evolutionary trends in their NPP, even though their biota composition was different. Therefore, this analysis was complemented with a JSI to assess whether the reconstructed herbivore species composition in each palaeocommunity differed among biogeographical regions during the late MIS 3. The JSI was based on presence–absence data and was calculated as follows:

where c is the number of shared species in both regions and a and b are the numbers of species that were only present in one of the biogeographical regions. Therefore, the higher the value the more similar the palaeocommunities of both regions were.

Chronological assessment

Pivotal to any hypothesis of Neanderthal replacement patterns by AMHs is the chronology of that population turnover. To this end, we used three different approaches to provide greater confidence in the results: BAMs, the OLE model and SPD of archaeological assemblages. As detailed below, each of these approaches provides complementary information about the MUPT.

First, we built a set of BAMs for the Mousterian, Châtelperronian and Aurignacian technocomplexes in each region during the MIS 3. As stated above, we compiled the available radiocarbon dates for Iberia between 55 and 30 kyr cal bp. However, not all dates or levels were included in the Bayesian chronology models. Radiocarbon determinations obtained from shell remains were incorporated in the dataset (dataset 1); however, the local variation of the reservoir age was unknown from 55 to 30 kyr bp. Because of uncertainties related to marine reservoir offsets, all BAMs that incorporated dates from marine shells were run twice: including and excluding these dates. All of the archaeological levels with cultural attribution issues or stratigraphic inconsistencies were excluded. The Supplementary Note provides a detailed description of the sites, levels and dates excluded and their justification. All BAMs were built for each technocomplex using the OxCAL4.2 software⁵⁵ and IntCal20 calibration curve⁵⁴.

Bayesian chronology models were built for each archaeological and palaeontological level. Then, the dates associated with each technocomplex were grouped within a single phase to determine each culture’s regional appearance or disappearance. Our interest was not focused on the chronological duration of the Mousterian, Châtelperronian and Aurignacian cultures, but on the probability distribution function of the temporal boundaries of these cultures in each region. Thus, this chronological assessment aims to provide an updated chronological frame for Neanderthal replacement by AMHs in Iberia. For this reason, we did not differentiate between proto- and early Aurignacian cultures, since both are attributed to AMHs.

In each BAM, we inserted into the same sequence the radiocarbon dates associated with a given technocomplex within a start and end boundary to bracket each culture, which allowed us to determine the probability distribution function for the beginning and end moment of each cultural phase⁶. The resolution of all models was set at 20 years. We used a t-type outlier model with an initial 5% probability for each determination, but when more than one radiocarbon date was obtained from the same bone remain, we used an s-type outlier model and the combine function. The thermoluminescence dating likelihoods were included in the models, together with their associated 1σ uncertainty ranges. When dates with low agreement (<60%) were detected, we evaluated two alternatives: first, these outlier dates were automatically down-weighted in the iterative Markov Chain Monte Carlo runs; second, these outliers were removed from the sequence and the models were re-run. To further assess the robustness of the chronology obtained for the MUPT in each biogeographical region, we performed two additional sensitivity tests. First, each model was run four times and the outputs were compared. Second, we repeated each simulation after removing the youngest and oldest dates in each technocomplex and region. It should be noted that the time intervals obtained with the boundary function were larger than those obtained with the date function, given that the date function informs about the duration of the phase, whereas the boundary function informs about when the phase started and stopped being created. All of the codes and outcomes of these sensitivity tests are available in the Supplementary Note.

Different studies have used the same Bayesian modelling approach to assess the MUPT in different regions of Iberia and Europe. However, this approach is subject to certain assumptions and limitations. Since different sites are included in a single phase, the BAMs assume that all dates are related, continuous and dependent on each other, which may result in relatively short time intervals, thus potentially biasing the estimation of the first and last appearance of each culture.

OLE models are used to determine the extinction time of species and, in recent years, this approach has proved to be robust in the estimation of the first and last appearance of Palaeolithic cultures^58,59,60. Unlike Bayesian statistical approaches, the OLE method does not penalize outlying data; therefore, the time intervals for the start and end of each culture are larger. OLE uses the first and last chronological occurrences of a given species or culture to assess how long it continued after the oldest or most recent confirmed occurrence⁵⁸. The mathematical formulation of the model is available elsewhere⁵⁸. Previous studies recommended using the five to ten oldest and youngest dates of each technocomplex to compute the first and last appearance, respectively^58,59. We used the median and range at 95.4% CI of each date obtained with the IntCal20 calibration curve⁵⁴ as input. To address the uncertainty of the chronometric determinations, each date from within each of the ten associated date ranges was randomly drawn from a normal distribution and this was used instead of the calibrated median dates⁵⁹. Such a randomly generated set of ages was assessed using the OLE method and the whole procedure was repeated 10,000 times⁵⁸. In the Thermomediterranean region, the Châtelperronian is represented by only one archaeological assemblage (for details, see Supplementary Note 1). Hence, the use of the BAMs is more appropriate in this case than the OLE approach.

OLE is a well-established methodology for estimating the first and last appearance of species, but this procedure shares some limitations with the BAMs. OLE assumes a continuation of the cultures after or before the latest or earliest currently known occurrences⁵⁸. This assumption is not valid, for example, in cases of rapid extinction, such as those derived from catastrophic events or large-scale migrations. Moreover, as in the BAMs, the OLE method does not reflect discontinuities in the fossil record. To overcome these shortcomings and to provide further support to the chronological aspects of this study, we also carried out an SPD analysis.

SPD is commonly used as a demographic proxy in Palaeolithic research. However, this method is subjected to certain assumptions^52,61. SPD analyses assume that population size has a positive correlation with the number of dates, sites or assemblages, and that the intensity of research and preservation of archaeological sites is nearly uniform across the region of study. These assumptions are not commonly met, so some filtering steps are necessary to ensure that SPD is a robust method with which to infer population dynamics. Following the recommendations of previous research⁶¹, several filtering processes were used. First, dates with large error ranges were removed by eliminating all chronometric determinations with a coefficient of variation ≥0.05 (ref. ⁶²). Second, to ensure that each occupational unit was not overrepresented, SPD estimations were based on the chronological distribution of each archaeological assemblage, and dates from the same archaeological level were merged using the R_Combine function from OxCal before calibration⁶³. Third, radiocarbon determinations obtained from shell remains were excluded because of the uncertainties related to the marine reservoir offsets, and the remaining dates were calibrated with the IntCal20 calibration curve⁵⁴. These filtering steps overcome some of the limitations of this method, but caveats are still necessary⁶¹. Therefore, outcomes obtained from the SPD have not been used to inform population size in this study. Instead, they have been used as an additional method to assess the duration and particularly the frequency of occupation of each culture.

Palaeoclimate reconstructions and validation

The model used to estimate NPP in the current study required the following climatological input data: monthly temperature (C°), precipitation (mm per month), incoming shortwave radiation (Wm⁻²) and rainy days (days per month). These palaeoclimate parameters were obtained from the HadCM3B-M2.1 coupled general circulation model with active atmosphere, ocean and sea ice components⁶⁴. This dataset extends back 60,000 years from 0 bp at 0.5° resolution on a monthly timestep in the Northern Hemisphere⁶⁴. Despite the high temporal resolution of this model, biases are frequently found when comparing global palaeoclimate models with local or regional palaeoclimate reconstructions⁶⁵. Therefore, we assessed the accuracy of the HadC3B-M2.1 simulations in the regions of interest and performed a bias correction before using this dataset in further palaeoecological modelling.

Assessing the accuracy of palaeoclimate models is challenging due to the low availability of observational datasets⁶⁴. Pollen assemblages constitute one of the most common palaeoenvironmental proxies used to evaluate the accuracy of palaeoclimate simulations^64,65. Several empirical palaeoclimate reconstructions are available for the Last Glacial Maximum and the mid-Holocene, but they are less common for MIS 3. Following previous studies, pollen-based palaeoclimate reconstructions were made with weighted averaging regressions⁶⁶. Weighted averaging regressions consider that plant species were more abundant near the climatic conditions most adapted to⁶⁷. The predictive functions used to compute the mean annual temperature (MAT) and mean annual precipitation (MAP) from the palynological record were derived from a training set of modern pollen taxa obtained from the Eurasian Modern Pollen Database version 2 (ref. ⁶⁸). This database contains pollen taxa recovered from more than 2,000 European and Asian localities and associated climatological conditions (Supplementary Fig. 4a). This training set was used to obtain temperature and precipitation transfer functions based on pollen subsets using weighted averaging regression techniques. Prediction errors were simulated by bootstrap crossvalidation (number of boot cycles = 500). Then, the transfer functions were applied to the fossil pollen recovered from 93 palynological assemblages from 51 dated archaeological levels that cover all of the biogeographical regions of Iberia (Supplementary Fig. 4b). The percentage of each pollen taxa was obtained from the species count whenever possible and from the published palynological diagrams in most cases (dataset 2).

Before estimating the MAT and MAP from the palynological record, some filtering criteria were applied to the dataset. First, only the palynological assemblages obtained from levels dated between 55 and 30 kyr cal bp were retained for analyses. Second, palynological assemblages with <100 pollens recovered and species with low representation (<5%) were excluded. Since terrestrial and aquatic taxa are unevenly affected by climatic conditions, it is recommended to also exclude all non-terrestrial pollen taxa, including ferns and non-pollen palynomorphs⁶⁷. After these modifications, the percentages of each taxon were readjusted and we ensured that all of the pollen taxa recovered from the archaeological record were in the training dataset of extant pollen species. After the crossvalidation process, the correlation coefficient (r²) used to estimate the MAT from the palynological record was 0.77 and the root-mean-square error (RMSE) was 4.14. However, the correlation coefficient was lower (r² = 0.50) and the RMSE higher (RMSE = 294.4) for the MAP (Supplementary Table 3). It is important to bear in mind that the palaeoclimatic reconstructions performed from the palynological record were not used in further palaeoecological reconstructions (for example, NPP estimations) but to check the predictions made from the HadCM3B-M2.1 coupled general circulation model and to assess whether bias corrections were necessary.

Different correction techniques have been proposed to rectify the observed biases in climate simulations, but the delta method has been shown to be the most robust approach to correct palaeoclimate simulations⁶⁵. According to the delta method, the bias in a specific area is obtained from the difference between present-day observed and simulated values. Therefore, bias-corrected temperature (T) in a particular area (x) at some time (t) in the past is estimated as follows:

where T_sim is the bias-corrected temperature in a particular area (x) during a specific moment in the past (t) according to the difference between the observed (T_obs) and predicted (T_sim) temperature for the present day.

Likewise, the bias-corrected precipitation (P_sim) for the past is estimated as follows:

To assess the suitability of this bias correction method, the estimated values of MAT and MAP from the palynological record were compared with MAT and MAP obtained from the HadCM3B-M2.1 coupled general circulation model before and after using the delta correction method. Present-day values of MAT and MAP were obtained from the Climate Research Unit version 4 dataset⁶⁹. The difference between the observed and predicted values of MAT during MIS 3 in the archaeological sites with pollen samples was 1.44 °C on average before using the bias correction method and 0.43 °C on average after using the delta correction method (Supplementary Fig. 5). The mean difference between the observed and predicted values of MAP was 1.74 mm per month before using the delta correction method and 0.82 mm per month after performing the bias correction. Accordingly, differences between the observed and predicted values were, on average, closer to zero after completing this bias correction (Supplementary Fig. 5). Moreover, there was a positive correlation between the observed and predicted values of MAT (P < 0.001; r = 0.68) and MAP (P < 0.001; r = 0.67) during MIS 3 (Supplementary Fig. 6). Accordingly, the bias-corrected values of temperatures and precipitations are in good agreement with the empirical reconstructions. These outcomes show the suitability of the delta correction method and the good correspondence between the observed and predicted rainfall and temperature values obtained from the HadCM3B-M2.1 coupled general circulation model⁶⁴ once the bias correction procedure had been performed.

Net primary productivity

NPP was estimated with the Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS) version 4.0. model. LPJ-GUESS (https://web.nateko.lu.se/lpj-guess/) is a process-based dynamic vegetation model that simulates the structure and composition of the land cover in terms of plant functional types (PFTs)⁷⁰. Each PFT is characterized by specific bioclimatic niche, growth form, leaf phenology, photosynthetic pathway and life history traits. The population dynamics of each PFT are determined by the competition for light, space and soil resources in each of a number of replicate patches for each simulated grid cell⁷⁰. There are different versions of LPJ-GUESS⁴⁷ and various PFTs can be incorporated into the model⁷¹. In the current study, we used the standard global PFT set described by Smith et al.⁷⁰. The model simulates the vegetation dynamics in multiple grid cells or patches to incorporate variation due to stochastic processes⁷⁰; however, these patches are simulated independently of each other, whereby the outcomes obtained do not change if the model is run in a specific patch or in many adjacent patches⁷⁰. Following Allen et al.⁷¹, the model was run without the nitrogen cycle nor nitrogen limitation because nitrogen deposition is unknown for the Pleistocene. For each grid cell, 100 replicate patches with an area of 0.1 ha were simulated. For details on the global parameters used, see Supplementary Table 4.

In the present study, LPJ-GUESS was run in cohort mode to estimate NPP in the surrounding area of each archaeological and palaeontological site providing information about herbivore guild composition and human subsistence strategies during the MUPT in Iberia. The NPP accrued was allocated to leaves, roots and woody PFTs according to the specific allometric relationships of each PFT, resulting in height, diameter and biomass growth⁷⁰. NPP was estimated at the end of each simulated year, but the model operates on a daily timestep and requires daily climate conditions. Monthly palaeoclimate drivers can be provided as input, and the model interpolates daily values from monthly values. Thus, the input climate variables for LPJ-GUESS were monthly temperature (C°), precipitation (mm per month), incoming shortwave radiation (Wm⁻²), and rainy days (days per month) for each time slice. These input climate data were obtained from the HadCM3B-M2.1 coupled general circulation model⁶⁴ after performing the delta bias correction procedure outlined above. The model also needs, as input data, the atmospheric carbon dioxide concentration (ppm) and the specific soil classes for each grid cell to incorporate texture-related variables affecting the hydrology and thermal diffusivity of the soils⁷⁰. CO₂ values were obtained from Lüthi et al.⁷² and soil types were obtained from Zobler⁷³.

All simulations were initialized with bare ground conditions (no biomass) and the model was spun up for 500 years until the simulated vegetation was in approximate equilibrium. This spin-up phase used monthly temperature, precipitation, incoming shortwave radiation and rainy days between 55 and 54.5 kyr bp⁶⁴. Thereafter, the model was run at a monthly resolution spanning the period between 54.5 and 30 kyr bp. The accuracy and robustness of the NPP estimations were assessed in two steps. First, we compared the observed and predicted values of NPP for the present day in the regions of interest. Second, we evaluated the sensitivity of the NPP estimations to the uncertainties in the palaeoclimate variables.

Regarding the simulations of NPP for the present day, we used the same modelling protocol described above, but the input data for the spin-up phase consisted of repeated sets of monthly temperature, precipitation, incoming shortwave radiation and rainy days data from between ad 1901 and 1930 obtained from the Climate Research Unit version 4 dataset⁶⁹. After that, the model simulated NPP for the period between ad 1982 and 1998 to allow a more straightforward comparison with observational datasets obtained from Imhoff et al.¹⁹. The results obtained show a significant positive correlation between the observed and predicted values of NPP (P < 0.001; r = 0.78) and the mean average error was 0.071 (Supplementary Fig. 7). These results provide confidence about the performance of this dynamic vegetation model and show its suitability for simulating NPP in the regions of interest.

As mentioned previously, the temperature and precipitation values obtained from the HadCM3B-M2.1 coupled general circulation model were in good agreement with the local empirical reconstructions made with the palynological record. However, it was still necessary to assess the extent the uncertainties in the palaeoclimate simulations could affect the estimated NPP. To this end, we additionally estimated NPP from the temperatures and precipitations obtained from two alternative general circulation models at 0.5° spatial resolution^74,75. Running the LPJ-GUESS to estimate NPP for 25,000 years (from 55,000–30,000 years bp) was a computationally intense task, so this sensitivity test focused on five specific stadial and interstadial phases: GI-12, GS-12, GS-9, GI-8 and GS-5. We selected these specific stadial and interstadial phases because they coincide with the end of the Mousterian in different regions of the Iberian Peninsula and because they represent moments of important fluctuations in NPP. The results obtained showed that the estimated NPP can vary up to 0.06 kg km⁻² yr⁻¹ on average in the Eurosiberian region, 0.07 kg km⁻² yr⁻¹ on average in the Mesomediterranean region, 0.05 kg km⁻² yr⁻¹ in the Supramediterranean region and 0.05 kg km⁻² yr⁻¹ in the Thermomediterranean region when the input palaeoclimate data are obtained from alternative general circulation models (Supplementary Table 5). However, it is important to note that the magnitude of NPP fluctuations between stadial and interstadial periods, as well as the regional differences in the estimated NPP remained constant (Supplementary Fig. 8). Therefore, this sensitivity test reinforces the robustness of the reconstructed temporal trends of NPP and the estimated differences in NPP between biogeographical regions.

Carrying capacity of herbivores

To date, numerous models have related the trophic structure of different ecosystems to their NPP. These studies were built on the empirical evidence for the relationship between NPP and herbivore abundance across diverse terrestrial ecosystems^76,77,78,79 (Supplementary Fig. 9). We obtained a predictive equation to estimate the total herbivore biomass (THB) that could be sustained by a given ecosystem with the data provided by different studies that cover a broad range of terrestrial ecosystems (dataset 3). An MM-type estimator for linear models was used to obtain the predictive function and to compute the 95% CI for the estimation. The obtained predictive function was:

where both THB and NPP are expressed in g m⁻² yr⁻¹ (Supplementary Table 6). The biomass of a herbivore population (B) is commonly estimated by multiplying its population density (D) by the mean adult body mass of the species. Accordingly, the THB of all herbivore species in a given ecosystem could also be expressed as follows:

where n is the number of herbivore species in the community, D_i is the population density of species i expressed in ind km⁻² and W is the mean body mass of both sexes in kg. Damuth⁸⁰ demonstrated that population density changes allometrically with body size across different ecosystems:

where D is expressed in ind km⁻², c is a constant and W is the mean body mass in kg. Among narrow taxonomic groups, the relationship between body mass and population density may differ from the exponent of −0.75 (ref. ⁸¹). Nevertheless, recent studies provide compelling evidence for the relationship of −3/4 on larger scales, demonstrating that, across broad taxonomic levels, the scaling factor of −0.75 fits the empirical observations⁸². Assuming this allometric relationship, the availability of trophic resources and the structure of regional communities, previous studies showed that herbivore abundances can be estimated^83,84,85. Following these earlier works, D can be substituted in equation (8):

Furthermore, c can be estimated as follows:

Accordingly, the biomass (B) of a specific herbivore population species (i) in a given ecosystem can be estimated with equation (12):

To validate this model, we compiled data of 516 extant herbivore population densities in protected areas covering a wide range of different ecosystems from the TetraDENSITY database⁸⁶ and Hatton et al.⁷⁸ (dataset 3). The body mass of each herbivore species was obtained from the PHYLACINE database⁵³. In 95.3% of the national parks and reserves, the estimated herbivore population densities were significantly correlated (P < 0.05) with the observed herbivore abundance, with correlation coefficients (r) ranging between 0.49 and 0.98 (Supplementary Fig. 10). The only two localities where there was no correspondence between the observed and predicted densities were tropical ecosystems with a low number of identified species because of the lack of density surveys (dataset 3). When the observed and predicted values were analysed together across all ecosystems, there was a significant positive correlation between the observed and predicted values (P < 0.001; r = 0.65) (Supplementary Fig. 11). These results confirm the validity of this macroecological modelling approach.

This approach assumes that the biomass of each herbivore population depends on the bottom-up processes of food chain regulation driven by NPP, the specific herbivore guild composition in a given ecosystem and the allometric relationships between body mass and population density. Even in glacial times, one would expect to find a positive relationship between NPP and THB, but megafauna could have different grass exploitation efficiencies than extant large herbivores, as previously suggested to explain the productivity paradox of the mammoth steppe ecosystems²³. Moreover, this approach does not consider other factors that affect herbivore abundances, such as species migratory movements or predatory pressure. For these reasons, the herbivore biomass derived from this modelling approach should not be interpreted as an estimation of the actual herbivore biomass, but rather as the carrying capacity of each herbivore population species in a specific region. It is worth noting that the current study did not aim to estimate precisely the biomass of each herbivore population, but rather to investigate to what extent MIS 3 climate changes affected the potential biomass availability for secondary consumers.

This modelling approach was applied to each MIS 3 herbivore palaeocommunity in the biogeographical regions of Iberia. A palaeocommunity was determined by the LFAs found in a specific biogeographical region during a particular time span without significant species turnover⁸⁷. However, different gaps may affect the number and type of herbivore species that conform to a given palaeocommunity (for example, uneven sampling effort or taphonomical biases). A method for partially correcting these drawbacks in palaeoecology is the minimum census technique. According to this method, the species composition of each palaeocommunity is not inferred from the species occurrences in each specific time interval (that is, stadial and interstadial phase), but from the chronological range of each species in each biogeographical region. However, this technique does not solve the issues related to the completeness of the species lists that conform to each palaeocommunity. To further assess the totality of the number of species in each palaeocommunity, a rarefaction analysis was performed. Rarefaction (interpolation) and prediction (extrapolation) curves are commonly used to assess species richness according to the sampling effort. This method assumes that when the number of samples increases, the species richness approaches an asymptote. Accordingly, the observed species richness was compared with the expected richness in a sample size of 100 LFAs in each palaeocommunity and a bootstrap method was applied (n = 500) to obtain the 95% CI for each diversity estimate. The outcomes obtained show that the actual species richness in each palaeocommunity was close to the asymptote (Supplementary Fig. 12) and, consequently, increasing the sample size of LFAs would not increase the number of species in each palaeocommunity (Supplementary Table 7) significantly.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com