Climate effects on archaic human habitats and species successions – Nature

2Ma simulation

We conducted the 2Ma simulation with the Community Earth System Model (CESM), version 1.2, at an ocean and atmosphere resolution of approximately 3.75° × 3.75°. The model uses bathymetry of the Last Glacial Maximum and time-varying forcings of greenhouse gases¹⁵, ice sheets¹⁵ and astronomical insolation conditions¹⁶. CESM1.2 has a relatively low standard equilibrium climate sensitivity (ECS) of 2.4 °C per CO₂ doubling, which lies outside the likely range of estimates⁴⁵ (3.7 ± 1.2 °C) obtained with other climate model simulations conducted as part of the Coupled Model Intercomparison Project, phase 6. However, this value is within the lower range of recent estimates compiled by the Intergovernmental Panel on Climate Change sixth assessment report⁴⁶ of Working Group 1, which identifies a very likely ECS range of 2–5 °C. To obtain a more realistic response to past long-wave radiative forcings in our palaeo-climate model simulation and to implicitly capture radiative effects of other CO₂-correlated forcings⁴⁷ from dust, vegetation, N₂O or CH₄, we therefore scaled the range of the applied CO₂ forcing¹⁵ by a factor of 1.5. The resulting effective ECS, which includes non-CO₂ greenhouse gas forcings, was in our case approximately 3.8 °C. Our result is in reasonable agreement with the Coupled Model Intercomparison Project phase 6 estimate and previous palaeo-climate estimates^18,19 of 3.2 °C, which were obtained from reconstructions of the global mean surface temperature and radiative forcings covering the last 784,000 years. Amplification of the CO₂ forcing in CESM1.2 led to a realistic representation of the amplitude of global mean, tropical and Antarctic temperature changes (Extended Data Figs. 1b and 2a, b) and to a simulated temperature range between Last Glacial Maximum and Late Holocene conditions of approximately 5.9 °C. This result is in close agreement with recent palaeo-proxy-based estimates²⁰. Similar to previous long-term transient climate model simulations conducted with Earth system models of intermediate complexity^7,48, the CESM1.2 simulations use an orbital acceleration factor of 5, which means that the 2-million-year orbital history is squeezed into 400,000 model years in CESM. The complete model trajectory is based on 21 individual chunks that were run in parallel, with each covering at least one interglacial–glacial cycle (Supplementary Table 2). Moreover, each chunk overlaps with the next chunk so that the issue of initial conditions and spin-up time can be evaluated properly. The final climate trajectory is obtained by combining the individual chunks and by using sliding linear interpolation in the chunk-overlap periods. The model simulation has been evaluated against numerous palaeo-proxy-based data (Fig. 1 and Extended Data Fig. 1). Unlike other Earth system models^49,50,51, the 2Ma simulation conducted with CESM1.2 does not generate strong internal millennial-scale variability such as that shown by Dansgaard–Oeschger cycles. The CESM1.2 data are provided on the climate data server of the Institute for Basic Science (IBS) Center for Climate Physics at https://climatedata.ibs.re.kr.

Topographic downscaling

The T31 spectral resolution of the 2Ma CESM1.2 simulation (approximately 3.75° × 3.75° horizontal resolution) is too coarse to properly capture important topographic barriers, which may have affected the dispersal and distribution of archaic humans. We applied simple downscaling of the simulated monthly surface temperatures T_s(t) onto a 1° × 1° horizontal grid by accounting for the difference in height Δh(t) between the ETOPO5 topographic dataset and the orographic forcing of the 2Ma experiment. The lapse rate-corrected temperatures were then calculated as T^*_s(t) = T_s(t) − gΔh(t), where g represents a constant average lapse rate of g = 6 °C per 1,000 m. The simulated rainfall p(t) was downscaled onto the high-resolution topography by accounting for temperature-dependent moisture availability through the Clausius–Clapeyron equation as p^*(t) = p(t)e^[17.625T*_s^/(T*_s^{+ 243.04) − 17.625}T_s^/(T_s^+ 243.04)].

A posteriori calculation of NPP

2Ma uses fixed plant functional types but a prognostic leaf area index. Therefore, we calculated the NPP a posteriori (Extended Data Figs. 5 and 9) using a simple empirical relationship among temperature, precipitation and tree fraction. The topographically downscaled temperature T^*_s (in degrees Celsius) and precipitation p^* (in millimetres per year) of the 2Ma simulation were used at every grid point to calculate the tree fraction⁵² as τ = 0.95{1 − e^[−β(T*_s ^− T_m^)]}p^*α/(p^*α + f), with the additional term f = be^[γ(T*_s ^− T_m^)], and the parameters β = 0.45, α = 3, b = 2.6 × 10⁶, γ = 0.155 and T_m = −15 °C; τ is capped between 0 and 1. Subsequently, the downscaled NPP can be calculated from an empirical model⁵³ as N* = {6,116[1 − e^{(−0.0000605}p^*)](1 − τ) + τ min(FP, FT)}f(CO₂), where the minimum (min) is taken over the mathematical terms FP = 0.551p^*1.055/e^(0.000306p^*) and FT = 2,540/[1 + e^{(1.584 − 0.0622}T*_s⁾]. The function f(CO₂) = [1 + 0.4ln(CO₂/280)/ln(2)] captures the bulk effect of CO₂ fertilization of plants⁵⁴ in the same way as the CLIMBER Earth system of intermediate complexity, and its time evolution is obtained from the transient CO₂ forcing of CESM1.2.

Extended dataset of archaeological and fossil hominin data

The SDM for the Homo genus was derived from a recent compilation of archaeological and fossil data²¹. The original data compilation published in 2020 (ref. ²¹) included 2,754 radiometric age estimates for fossil hominin occurrences, each accompanied by the confidence interval around the estimate, the fossil site name and the archaeological layer within the site (where available) from which the dated sample was derived, the geographical coordinates of the site and the possible attribution to one or more than one Homo species. Confident attributions to a single species generated a core record, whereas instances with multiple attributions formed an extended record. Six different species were recognized: H. habilis, H. ergaster, H. erectus, H. heidelbergensis, H. neanderthalensis and H. sapiens. The updated record, as shown in Supplementary Table 1, contains 3,245 data entries restricted to the temporal age interval of 2 Ma–30 ka; those from Australia and the Americas were excluded. Further, we combined H. habilis and H. ergaster into a single African Oldowan toolmaker species. Each occurrence is attributed to a given species depending on the presence of unambiguous anatomical remains, either singly or in connection to a specific lithic tool industry. This helped to guide identification if this was not otherwise possible from the bones and teeth alone (398 entries, 12%), the age limits of the individual species or the stone tool industry. For example, an occurrence in Africa older than the first appearance of H. sapiens at Jebel Irhoud⁵⁵ yet younger than the first appearance of H. heidelbergensis at Melka Kunture⁵⁶ is attributed to H. heidelbergensis. Moreover, French Mousterian stone tools have been unambiguously assigned to H. neanderthalensis, whereas Aurignacian tools were attributed to H. sapiens. When these criteria were applied, the core record included 94.5% of the attributions, 48.5% of which refer to H. neanderthalensis and 37.5% of which refer to H. sapiens. Where neither of these criteria was met (in the original compilation, the SDM acknowledges attribution uncertainty), we accounted for this by testing the stability of our results with respect to different versions of the SDM (Extended Data Fig. 10). For example, transitional industries (for example, the Levantine Mousterian or Lincombian–Ranisian–Jerzmanowician industries) received multiple attributions because they fit either H. sapiens or H. neanderthalensis in terms of toolmaker identity^57,58. A detailed explanation of this approach is provided as supplementary material for ref. ²¹ (https://ars.els-cdn.com/content/image/1-s2.0-S2590332220304760-mmc1.pdf).

A second source of uncertainty stems from dating. Although approximately 50% of the entries refer to the ¹⁴C method (>90% of which are based on accelerator mass spectrometry), other dating methods such as electron spin resonance (14% of the sample), thermoluminescence (12%) and optically stimulated luminescence (12%) are less precise than radiocarbon dating. Nonetheless, multiple datings are present for individual fossil sites, even within a single stratigraphic layer at the site. To account for uncertainties in species attribution and age, we ran our analyses according to the four different approaches described below.

1. 1.

   Multiple dates, tier 1. Only the core record, which excludes entries with uncertain species attributions, and all age estimates available for each archaeological layer are used. Multiple age estimates per layer are possible, and the age uncertainty for each is included in our analysis. This subdivision includes 3,060 data entries. Although the main analysis in our study is based on this case, we need to consider possible sampling biases due to the higher weights given to archaeological layers with multiple dates (Figs. 1–4 and Extended Data Fig. 10).

2. 2.

   Multiple dates, tier 2. The extended record, in which ambiguous species attributions are treated by randomly choosing among the possible candidate species, is used along with multiple age estimates (including uncertainties) per layer. This subdivision includes 3,245 (all) data entries (Extended Data Fig. 10).

3. 3.

   Single date, tier 1. Multiple age estimates for a single archaeological layer are combined in this approach to provide a minimum and maximum age for the layer. Each archaeological layer has only one entry, thereby eliminating possible sampling biases in the estimation of our CEM. This subdivision includes 1,567 data entries (Extended Data Fig. 10).

4. 4.

   Single date, tier 2. Age estimates for archaeological layers are treated as those in the single date, tier 1 category except that the extended record rather than the core record is used. This subdivision includes 1,652 data entries (Extended Data Fig. 10).

We acknowledge that our species subdivisions may be controversial and that these do not necessarily require constancy of morphology, habitat and behaviour. However, even though some species attributions such as H. heidelbergensis could be questioned, we remain confident that the majority of the record presents little challenge considering that 86% of the core data belong to the well-defined, widely accepted H. neanderthalensis or H. sapiens record and tool-making traditions. Thus, even though some species attributions might be considered invalid, widely accepted constraints are used. Clearly, to the best of current knowledge, 500,000-year-old remains in Africa can belong toneither H. sapiens nor H. habilis⁵⁹, irrespective of whether the name H. heidelbergensis is considered appropriate. To further reduce uncertainties, we tested the robustness of our main findings with four alternative scenarios (Extended Data Fig. 10) for species attribution and dating and excluding uncertain and poorly dated species (for example, Homo floresiensis, Homo naledi, Homo bodoensis, Homo longi and Denisovans), which are restricted to too few fossil sites for which no climatic variability can possibly be ascertained or do not currently include any other locality or remains in their definition. The final species assignments used in our study should be interpreted here as plausible working hypotheses.

Mahalanobis CEM

To derive the CEM (Extended Data Fig. 5) that best characterizes the habitable conditions for hominins, we chose four key climatic variables: annual mean temperature and precipitation (T*_am and P*_am, respectively), annual minimum precipitation (P*_min) and terrestrial NPP (N*). Obtained as 1,000-year downscaled averages (1° × 1° horizontal resolution), these variables, which relate to physiological constraints for hominin survival and the availability of food resources, are combined as a four-dimensional climate environment vector C(t) = (T*_am, P*_am, P*_min, N*) with 2,000 values in time (t) corresponding to 1,000-year (200-year) orbital (model) means from the 2Ma simulation. The fossil and archaeological data entries for the five individual hominin groups are described in the previous section. Although our main analysis focuses on the multiple date, tier 1 case (Methods and Supplementary Table 1), the robustness of our results was tested against other ways of treating species and age model uncertainties (Extended Data Fig. 10). The data entries are represented by their longitude λ_j,i and latitude φ_j,i coordinates and the respective average age t_j,i and age uncertainties Δt_j,i with i = 1, …, 5 corresponding to the five hominin groups. We defined the fossil state vector as zⁱ = (λ_1,i, φ_1,i, t_1,i, Δt_1,i, …, λ_n,i, φ_n,i, t_n,i, Δt_n,i) with n_i representing the total number of fossils in each group during the past 2 million years. We then built the matrix D from the four-dimensional climate data subsampled at the fossil data sites and the corresponding nearest ages. Age uncertainties were considered through a Monte Carlo sampling method, which expanded the length of the overall data vector. We obtained Dⁱ (4 × N_i matrix for each i = 1, …, 5) for each hominin group as Dⁱ = (T*_am(zⁱ), P*_am(zⁱ), P*_min(zⁱ), N*(zⁱ)). We then calculated the Mahalanobis squared distance model²³ for each group using ζ_i²(Dⁱ) = (Dⁱ − <Dⁱ>)^TS⁻¹ (Dⁱ − <Dⁱ>), where <…> represents the sample mean value and S⁻¹ is the inverse co-variance matrix obtained from the data. The Mahalabonis squared distances ζ_i were then transformed into a cumulative chi-squared distribution χ²_CDF in the four-dimensional climate data space C. When using 4 degrees of freedom²³, the corresponding probability H(C) = 1 − χ²_CDF(C) represents the likelihood of finding a fossil for a specific quadruplet within the four-dimensional climate data space in the HSM. We interpreted H as a probability, which we refer to as habitat suitability. Given the temporal evolution of C for every grid point of the downscaled 1° × 1° data over the last 2 million years, we were able to calculate the spatiotemporal habitat suitability for each downscaled grid point (x,y,t) in the model as H(x,y,t) = H(T*_am(x,y,t), P*_am(x,y,t), P*_min(x,y,t), N*(x,y,t)). The stability of the HSM was tested by using different dimensionalities and combinations of climate parameters such as annual mean and seasonal range of temperature and precipitation and annual mean and minimum values of temperature and precipitation. The key conclusions of our study remained essentially unchanged. Moreover, we tested the stability of our results against the omission of hominin sites with ambivalent original species attributions (multiple date, tier 2) and different treatment of archaeological ages (single date, tiers 1 and 2). The calculated H(x,y,t) was qualitatively very similar for the four different cases (Extended Data Fig. 10). Therefore, our main conclusions remain robust with respect to uncertainties in species attribution and archaeological layer dating.

Random climate trajectory

To address the question of whether the actual climate trajectory influenced the distribution of fossil and archaeological data, we developed a CEM and HSM in which fossil data were assigned to randomly chosen climate states from the CESM1.2 simulation under the constraint that the climate range selected must overlap with the total fossil age range of the respective species. We randomized the time variability of the four-dimensional climate data vector (annual mean temperature, annual mean precipitation, minimum precipitation and NPP) while keeping the co-variability among the climate vector components, as well as the mean state, invariant. The original HSM (H), which is based on the real trajectory of CESM1.2, and the model (H_scr) that we trained from a scrambled trajectory were then compared. The time-averaged differences between the models for H. sapiens, H. neanderthalensis and H. heidelbergensis were then interpreted as an indication of whether the realistic climate evolution influenced the observed hominin distributions in space and time relative to a system that maintains its orbital climate co-variance and mean state (Extended Data Fig. 8a–c) but does not consider the exact time evolution of glacial–interglacial and orbital cycles. The time-averaged difference between H(x,y,t) in the original HSM and H_scr(x,y,t) in the HSM derived from time-randomized climate data was then tested at each grid point using a paired t-test.

Reporting summary

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