Climate effects on archaic human habitats and species successions

It has long been believed that climate shifts during the last 2 million years had a pivotal role in the evolution of our genus Homo1–3. However, given the limited number of representative palaeo-climate datasets from regions of anthropological interest, it has remained challenging to quantify this linkage. Here, we use an unprecedented transient Pleistocene coupled general circulation model simulation in combination with an extensive compilation of fossil and archaeological records to study the spatiotemporal habitat suitability for five hominin species over the past 2 million years. We show that astronomically forced changes in temperature, rainfall and terrestrial net primary production had a major impact on the observed distributions of these species. During the Early Pleistocene, hominins settled primarily in environments with weak orbital-scale climate variability. This behaviour changed substantially after the mid-Pleistocene transition, when archaic humans became global wanderers who adapted to a wide range of spatial climatic gradients. Analysis of the simulated hominin habitat overlap from approximately 300–400 thousand years ago further suggests that antiphased climate disruptions in southern Africa and Eurasia contributed to the evolutionary transformation of Homo heidelbergensis populations into Homo sapiens and Neanderthals, respectively. Our robust numerical simulations of climate-induced habitat changes provide a framework to test hypotheses on our human origin.

It has long been believed that climate shifts during the last 2 million years had a pivotal role in the evolution of our genus Homo [1][2][3] .However, given the limited number of representative palaeo-climate datasets from regions of anthropological interest, it has remained challenging to quantify this linkage.Here, we use an unprecedented transient Pleistocene coupled general circulation model simulation in combination with an extensive compilation of fossil and archaeological records to study the spatiotemporal habitat suitability for five hominin species over the past 2 million years.We show that astronomically forced changes in temperature, rainfall and terrestrial net primary production had a major impact on the observed distributions of these species.During the Early Pleistocene, hominins settled primarily in environments with weak orbital-scale climate variability.This behaviour changed substantially after the mid-Pleistocene transition, when archaic humans became global wanderers who adapted to a wide range of spatial climatic gradients.Analysis of the simulated hominin habitat overlap from approximately 300-400 thousand years ago further suggests that antiphased climate disruptions in southern Africa and Eurasia contributed to the evolutionary transformation of Homo heidelbergensis populations into Homo sapiens and Neanderthals, respectively.Our robust numerical simulations of climate-induced habitat changes provide a framework to test hypotheses on our human origin.
During the past 5 million years (Ma), a gradual transition in climate conditions has occurred from the warmer and wetter Pliocene (5.3-2.6 Ma) to the colder and drier Pleistocene (2.6-0.011Ma).During this time, tropical savannahs and open grasslands expanded in central-eastern Africa 4 , which, according to the savannah hypothesis 5 and variants thereof 6 , contributed to the early evolution of our human ancestors.Milankovitch cycles in solar insolation and climate (Extended Data Figs.1-3), particularly the eccentricity-modulated precessional cycle (Extended Data Fig. 1a), further created multiple human migration corridors from sub-Saharan Africa into northern Africa, the Arabian Peninsula and Eurasia [7][8][9][10] .The existence of these corridors is well supported by fossil, archaeological 9 and genetic 11 evidence.A possible effect of astronomical forcings on early hominin evolution has been suggested in the context of the variability selection hypothesis 3,12,13 , which posits that early hominin evolution, selection and speciation were influenced by alternating periods of high and low variability in climate and resources.
To better quantify the impact of spatially heterogenous orbital-scale climate variability 14 (Extended Data Fig. 4) on human evolutionary transitions, we conducted an unprecedented transient global coupled general circulation model (CGCM) simulation covering the global climate history of the last 2 Ma (henceforth referred to as the 2Ma simulation).2Ma is based on the state-of-the-art Community Earth System Model version 1.2 (CESM1.2) in 3.75° × 3.75° horizontal resolution forced with ice-sheet distribution and elevation as well as CO 2 evolution obtained from another transient intermediate-complexity model simulation 15 and astronomical insolation changes 16 (Methods).2Ma, which uses an orbital acceleration 17 factor of 5, reproduces key palaeo-climate records such as tropical sea surface temperatures, Antarctic temperatures, the eastern African hydroclimate and the East Asian summer monsoon in close agreement with palaeo-reconstructions (Extended Data Figs. 1  and 2), which supports the realism of our CGCM-based simulation.Glacial-interglacial variability is characterized by a global mean temperature amplitude of approximately 2-3 °C (5-6 °C) during the Early (Late) Pleistocene (Extended Data Fig. 2a), which is consistent with estimates from an Earth system model of intermediate complexity 15 and palaeo-climate data constraints [18][19][20] .
To quantify the relationship between climate and the presence of hominin species, we built a climate envelope model (CEM; Methods and Extended Data Fig. 5).This CEM was derived from an extended version of a previously published species database (SDB) 21,22 composed of geochronologically constrained hominin fossils and archaeological layers containing lithic industries (Fig. 1a-e and Methods) and topographically downscaled (1° × 1° grid) 1,000-year averaged data of climate Article variables from 2Ma, which are relevant for human survival.These factors included annual mean precipitation, temperature, yearly minimum precipitation and net primary productivity (NPP; Methods).The 3,245 data entries of the extended SDB (Supplementary Table 1) contain information about location, age, age uncertainty and hypothesized species, selected among early African Homo (combining Homo habilis and Homo ergaster as one group), Eurasian Homo erectus, Homo heidelbergensis, Homo neanderthalensis and Homo sapiens.The spatiotemporal climate fields of the 2Ma simulation (see Extended Data Figs. 3 and 4 for select locations) were extracted for the species-presence locations and ages in the SDB and were then statistically aggregated as a CEM.Subsequently, using the Mahalanobis metric 23,24 and the spatiotemporal climate evolution in 2Ma, we derived a habitat suitability model (HSM; Methods) for each species, which quantifies the probability of finding fossil and/or archaeological evidence of the species at a given time and geographical location.
The key goals of our study were (1) to address how past climate changes have affected archaic human habitats; (2) to test whether the current fossil and archaeological records (location and age of each hominin species) have been affected by the orbital-scale evolution of our climate system; (3) to identify common climate envelopes and therefore potential contact zones of hominin groups; and (4) to identify linkages between regional climate shifts and evolutionary diversification.

Time-averaged habitats
To illustrate the connection between climate and the temporal and geographical extent of hominin species, we focused on habitat suitability calculated from the CEM.The simulated time-averaged maps of hominin habitat suitability (Fig. 1a-e) exhibit several interesting features.In particular, the suitable habitat for early African Homo (Fig. 1e) is composed of relatively narrow corridors that begin in southern Africa and run northward throughout the rift valley, straddle the Intertropical Convergence Zone and cut across southern Africa in a northwestsoutheast direction.Such a limited range and high spatiotemporal heterogeneity of habitat suitability are consistent with high levels of environmental specialization and sensitivity to regional environmental perturbations, such as eccentricity-modulated precessional cycles (Fig. 1f-i).Even though we included only Eurasian fossils and artefacts for H. erectus in the HSM, the predicted global habitat suitability of this species is far more extensive than that of any other hominin species analysed here (Fig. 1d).This is consistent with the concept that H. erectus was, on an evolutionary timescale, a flexible generalist who roamed Earth for more than 1 Ma and inhabited a wide range of different environmental conditions (Extended Data Fig. 6).Even though H. erectus and early African Homo fossil records are treated as geographically disjunct (Fig. 1d, e), there is still regional overlap in their climate  Homo heidelbergensis  1).The time-averaged habitat suitability (blue to white shading) covering the period of respective hominin presence can be interpreted in terms of probability (Methods), with values ranging from 0 (habitat unsuitable) to 1 (habitat extremely suitable).Coloured circles represent the locations of fossils and/or archaeological artefacts associated with the five hominin groups.f-i, Time series for precession (blue) and eccentricity (f) and simulated regional habitat suitability at selected sites of archaeological interest for H. habilis and H. ergaster (treated jointly as early African Homo), H. heidelbergensis and H. sapiens (g-i).The centre locations of a 4° × 4° average area include Jebel Irhoud (34° N, 4° W), a region near Lake Turkana (0° N, 34° E) and Blombos cave (34° S, 21° E).
envelopes inside Africa (Fig. 2b and Extended Data Fig. 7d), which is consistent with a deeper ancestral linkage between these two groups 25 .For H. heidelbergensis, we observed a time-averaged habitat suitability pattern that was qualitatively similar to that of H. neanderthalensis (Fig. 1b, c).By comparing the climate niches of H. sapiens (Fig. 1a) with those of other hominin species, we determined that our own species was best equipped to cope with dry conditions (Fig. 1g and Extended Data Fig. 6c).This extended climatic tolerance of H. sapiens was introduced into the CEM by a group of fossils and archaeological artefacts located in northeastern Africa, the Arabian Peninsula and the Levant (Fig. 1a and Supplementary Table 1).This tolerance of dry conditions greatly enhanced the mobility of H. sapiens, which may have further facilitated the documented multiple-wave dispersals into Eurasia across the Sinai passage or Bab-el Mandeb strait into the Levant (Extended Data Fig. 6c) and the Arabian Peninsula 9 , respectively.

Climate impacts on species distributions
The temporal evolution of our HSM exhibits pronounced Milankovitch cycles (Figs.1g-i and 3, Extended Data Fig. 6 and Supplementary Videos 1-5).Tropical regions are characterized mostly by precessional cycles, which are modulated by eccentricity cycles of 80-120 thousand years (ka) and 405 ka (Fig. 1f), whereas extratropical locations show a stronger component of 80-120 ka due to CO 2 and ice-sheet forcings (Extended Data Fig. 6b, d).Notably, regional climate changes and the resulting habitat changes were driven not only by the interplay of local forcings but also by remote effects such as eastern equatorial Pacific temperature changes, as suggested recently 26 by a synthesis of African hydroclimate proxy records and tropical sea surface temperature reconstructions.
To further test whether orbital-scale climate variability influenced the observed spatiotemporal distribution of hominin species, we recalculated the CEM for each species using the fossil and archaeological data in combination with a time-scrambled trajectory of the CESM1.2simulation (Methods).The resulting new CEM is different from the original one in that it assigns different, randomized temporal climate states to the fossil and archaeological data while maintaining the overall regional co-variability of the climate components and the long-term mean state.By comparing the long-term mean difference in the habitat suitability projections of the null-hypothesis model with the original one, we could then ascertain whether Milankovitch cycles influenced the distribution of fossils and archaeological sites on a regional level.The results for H. sapiens, H. neanderthalensis and H. heidelbergensis (Extended Data Fig. 8a-c) showed statistically significant differences

Article
(P < 0.05, paired t-test) in the calculated habitat suitability, with values of more than 0.05 in magnitude attained when comparing the unshuffled and shuffled models over parts of Asia, Europe and Africa.This documents that the orbital-scale trajectory had an important role in determining where and when hominin species lived.

Species successions
To identify locations where potential succession or speciation of hominin groups may have taken place, we calculated the species overlap as the co-variance of habitat suitability between the different hominin groups (Fig. 2 and Extended Data Fig. 7).We assumed that species that interacted with or emerged from each other probably shared similar regional climate envelopes, at least during their transition time.
For H. neanderthalensis and H. heidelbergensis, the highest values of niche overlap were found in Europe (Fig. 2b and Extended Data Fig. 7b), which also hosts archaeological artefacts and fossils from both species 27 (Fig. 3b) and has been regarded as the 'birthplace' of Neanderthals 28 .By comparing the zonal mean overlap for H. sapiens and H. heidelbergensis with their respective fossil and archaeological sites, we identified two key areas with climatic conditions that were suitable for joint occupancy outside Europe: central-eastern Africa and southern Africa (Fig. 3).In addition to habitat overlap (Fig. 2), we calculated the regional habitat similarity as an indicator for potential evolutionary transitions such as baseline evolution or speciation events (Extended Data Fig. 7a).A more detailed analysis into the simulated regions of orbitally varying species overlap indicated two pronounced periods of reduced habitat suitability in southern Africa for H. heidelbergensis at 415-360 ka and 340-310 ka (Fig. 3d).These prolonged eras of climatic stress were further characterized by low probabilities for fossil and archaeological records in subequatorial Africa.Subsequently, from 310 ka to 200 ka, high values of habitat suitability correlated with the first evidence of H. sapiens in southern Africa in terms of both fossils and archaeological artefacts [29][30][31] (Figs.2a and 3d) as well as presence of the earliest mitochondrial DNA lineage (L0) of southern African origin 32 .The disappearance of H. heidelbergensis from Africa could potentially be explained by progressive evolution of H. heidelbergensis into H. sapiens, which would be consistent with the presence of their respective fossils and archaeological artefacts at about 200-300 ka (Supplementary Table 1) and their similar values for regional habitat  suitability (Fig. 3d).By contrast, a larger habitat discrepancy between H. heidelbergensis and H. sapiens (Fig. 3c) occurred in central Africa, indicating that gradual species transition or diversification is less likely to have occurred in this region than in southern Africa, at least from a climate envelope perspective.Another major climate disruption in southern Africa around 210-200 ka (Fig. 3d) during the austral summer perihelion (Fig. 3a) could have increased the regional environmental stress on H. sapiens, leading to dispersal and, subsequently, genetic diversification.This timing is consistent with the presence of the first known mutation event that occurred in our reconstructed common mitochondrial ancestry 32 , even though considerable uncertainties in dating and methodology still exist 33 .Overall, our analysis suggests that the emergence of H. sapiens and the gradual disappearance of H. heidelbergensis in southern Africa coincided with long-term climatic anomalies during Marine Isotope Stages 11 and 9.

Speciation and dispersal
We combined a transient Pleistocene climate model simulation with an extensive compilation of hominin fossils and archaeological artefacts to study the environmental context of hominin evolution.On the basis of the resulting HSM and palaeogenetic evidence 34,35 , we propose the following scenario (Fig. 4): about 850-600 ka, H. heidelbergensis, which may have originated from H. ergaster in eastern Africa (Extended Data Fig. 7e), split into southern and northern African branches, the latter of which included northern African and Eurasian populations.
The intensified dispersal into off-equatorial regions may have occurred during periods of high eccentricity around 680 ka and 580 ka, which increased habitat suitability in otherwise unhospitable regions (Fig.  3b), gradually evolved into H. neanderthalensis between 400 ka and 300 ka.Side branches to northwestern Africa, back-propagation, multiple dispersals 37 , interbreeding 38 and subsequent speciation 39 may have further complicated the picture.
Recent studies have suggested that the sequence of hominin speciation events and the long-term positive trend in brain size may have been linked to past climatic shifts in Africa 40 .Our analysis supports the notion of strong Milankovitch cycles in early hominin habitat suitability in central Africa (Fig. 1h).Moreover, during the early Pleistocene (2-1 Ma), early African Homo populations occupied two main habitats: one in central-eastern Africa and the other in southern Africa (Fig. 2b).On average, these groups preferred geographical regions that were characterized by relatively stable NPP values of 200-380 gC m −2 per  2a).The proposed divergence at 400-300 ka of H. heidelbergensis into Neanderthals in Europe is consistent with recent genetic estimates 34 .This scenario is also in agreement with Neanderthal whole-genome data 44 that suggest a population split between Neanderthal-Denisovan and modern human lineages between 550 and 765 ka and a divergence between Neanderthals and Denisovans around 445-473 ka.Possible eccentricity-modulated windows for early non-coastal north-south migrations occurred around 700 ka and 600 ka during periods of high eccentricity, according to the calculated HSM (see inset time series for 4° × 4° averages centred near 21° N, 31° E and 20° S, 31° E).
Article year (Fig. 2c).Within Africa, early African Homo populations adapted mostly to local orbital-scale variations in climate and NPP (Fig. 1g-i and Extended Data Figs.3, 4 and 9), as reflected also in their habitat suitability.After the mid-Pleistocene transition and with the emergence of H. heidelbergensis between approximately 885 and 865 ka, the dynamics again changed remarkably.H. heidelbergensis began to migrate into Eurasia and other regions, encountering along their journey a much wider spatial range of NPP, from extremely low values of 20 gC m −2 per year to values exceeding 600 gC m −2 per year (Fig. 2c).These migrating groups crossed large spatial gradients in climate and NPP that far exceeded the temporal ranges in NPP experienced by their more stationary Early Pleistocene predecessors.This transition to global wanderers about 0.8-0.6Ma must have required H. heidelbergensis to acquire new adaptation skills, which in turn also strengthened their ability to further expand their geographical range, thereby providing a strong positive migration-climate adaptation feedback.Our analysis clearly shows that H. erectus had already undergone such a transition from regional dweller to early global wanderer before 1.8 Ma (Fig. 2c).Together with the H. heidelbergensis evidence, this indicates that dispersals from Africa always involved an adaptive shift, either biological or cultural, to wider climate envelopes.Therefore, to understand hominin evolution during the Pleistocene, the full spatial and temporal complexity of the climate signal and the corresponding habitat suitability must be considered.

Discussion
The main conclusions of our analysis are robust with respect to the existing uncertainties in species attribution, particularly for the period from 1 to 0.3 Ma, and the dating of archaeological layers, as demonstrated by key HSM calculations with four different scenarios that accounted for these factors (Methods and Extended Data Fig. 10).Although our study is based on species-stratified fossil and archaeological input data, our calculation of species overlap as HSM co-variability allowed us to treat potential species transitions and successions in human evolutionary history quantitatively and to identify their spatiotemporal characteristics.To the best of our knowledge, such research has not been reported thus far.The HSM captures regionally distributed patchworks of habitable areas in agreement with a general multiregional perspective 41 (Figs. 1 and 4).According to our CEM, southern and eastern Africa as well as the region north of the Intertropical Convergence Zone emerge as potential long-term refugia for various types of archaic humans.As the climate changed on orbital timescales, these refugia shifted geographically, creating population patterns with greater complexity.Further analysis of the pan-African connectivity of refugia in our HSM dataset, as shown in the inset in Fig. 4, will increase understanding of hominin dispersal, interbreeding and cladogenetic transitions as well as potential cultural exchanges.
In summary, we demonstrated that astronomically forced climate shifts were a key factor in driving hominin species distributions 42 and dispersal and were probably important for diversification 43 .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-04600-9.
Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

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 15 , ice sheets 15 and astronomical insolation conditions 16 .CESM1.2 has a relatively low standard equilibrium climate sensitivity (ECS) of 2.4 °C per CO 2 doubling, which lies outside the likely range of estimates 45 (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 46 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 2 -correlated forcings 47 from dust, vegetation, N 2 O or CH 4 , we therefore scaled the range of the applied CO 2 forcing 15 by a factor of 1.5.The resulting effective ECS, which includes non-CO 2 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 2 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 20 .Similar to previous long-term transient climate model simulations conducted with Earth system models of intermediate complexity 7,48 , the CESM1.2simulations 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 .

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 21 .The original data compilation published in 2020 (ref. 21) 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 55 yet younger than the first appearance of H. heidelbergensis at Melka Kunture 56 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. 21(https://ars.els-cdn.com/content/image/1-s2.0-S259033222 0304760-mmc1.pdf).
A second source of uncertainty stems from dating.Although approximately 50% of the entries refer to the 14 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.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. 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. 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. 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 to neither H. sapiens nor H. habilis 59 , 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 i = (λ 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 i (4 × N i matrix for each i = 1, …, 5) for each hominin group as D i = (T* am (z i ), P* am (z i ), P* min (z i ), N*(z i )).We then calculated the Mahalanobis squared distance model 23 for each group using where <…> represents the sample mean value and S −1 is the inverse co-variance matrix obtained from the data.The Mahalabonis squared distances ζ i were then transformed into a cumulative chi-squared distribution χ 2 CDF in the four-dimensional climate data space C. When using 4 degrees of freedom 23 , the corresponding probability H(C) = 1 − χ 2 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.2simulation 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
and orbital forcing from astronomical estimates 16 .The simulated surface temperatures and precipitation from 2Ma on a ~3.75-degree horizontal grid are then downscaled to a 1x1 degree horizontal grid by including lapse-rate corrected topographic features.Net Primary Production is calculated from empirical parameterizations of the downscaled temperature and precipitation fields (Extended Data Fig. 8d-f).Using an extended fossil and archaeological database (Supplementary Table 1) in combination with the downscaled annual mean temperatures, annual mean rainfall, minimum rainfall and net primary productivity, a statistical Mahalanobis distance-based climate envelope model (CEM) is derived.The model is then forced for every land grid point on a 1x1 degree grid with the temporal evolution of the downscaled climate variables to obtain the temporal evolution of the habitat suitability (HSM) for each of the 5 hominin groups and at every grid point.The impact of resolution on key features in simulated net primary production is further illustrated in Extended Data Fig. 8d-f.Extended Data Fig. 6 | Orbital effects on regional habitability in Eurasia.Timeseries for precession (blue) and eccentricity (a) and simulated regional habitat suitability at selected sites of archaeological interest for Homo habilis, Homo ergaster, Homo erectus, Homo heidelbergensis, Homo neanderthalensis, atmosphere resolution, obtained from empirical NPP model using 1000-year late Holocene average of total precipitation and surface temperatures simulated by the CESM1.22Ma experiment; middle: same as left, but using altitude-corrections for temperature and precipitation as downscaling onto a 1x1 degree target grid, showing the emergence of key topographic features in Africa.This resolution was chosen in the calculations of the climate envelope model; right: for illustration, same as middle but for a 0.25x0.25 degree horizontal target grid.The qualitative gain in terms of regional details from T31 to 1x1 degree resolution outweighs the additional gain going from 1x1 degree resolution to a 0.25x0.25 degree grid.
Extended Data Fig. 10 | Robustness with respect to age and species assignment uncertainties.Upper panels, same as for Fig. 1c, but for the 4 different cases single date-tier 1,2 and multi date-tier 1,2 (Methods).Multi-date tier 1 is the standard case used in our analysis.Lower panels, same as Fig. 3b-d, but calculated also for the 4 cases of assessing species and age uncertainties.

Fig. 1 |
Fig. 1 | Hominin species habitat suitability.a-e, African-Eurasian species distribution calculated from a Mahalanobis distance model using fourdimensional climate envelope data of topographically downscaled temperature, precipitation and NPP changes simulated by 2Ma (Methods) and the locations and ages of fossil and archaeological sites (Supplementary Table1).The time-averaged habitat suitability (blue to white shading) covering the period of respective hominin presence can be interpreted in terms of probability (Methods), with values ranging from 0 (habitat unsuitable) to

Fig. 2 |
Fig. 2 | Habitat overlap, succession and geographical distribution of fossils.a, Green shading represents a Hovmoeller (time-latitude) diagram of the zonal mean of the spatial scalar product of habitat suitability for H. sapiens and H. heidelbergensis.Circles represent the corresponding average age and latitudinal distributions of fossils and archaeological artefacts.High values of habitat overlap coinciding with joint presence of fossils indicate possible locations of hominin interaction, diversification and possibly speciation.b, Same as a but for H. neanderthalensis and H. heidelbergensis (right side) and H. erectus and H. ergaster-habilis.'Out-of-Africa' migration periods are marked as OOA.Potential regions for gradual diversification and transformation are indicated by dual-coloured boxes.c, NPP (Methods) for each fossil and archaeological site (coloured circles), selected by averaging the 2Ma NPP data in a 6° × 6° vicinity around the fossil sites and for their respective fossil ages.The size of the circle represents the great-circle (Haversine) distance to a grid point in central-eastern Africa (4° N, 36° E), with larger circles indicating closer proximity to this location.

Fig. 3 |
Fig. 3 | Hominin species succession in Europe and southern Africa.a, Eccentricity (orange) and precession (blue) indices from Fig. 1f.b, Habitat suitability calculated from the CEM for H. heidelbergensis (dark red curves) and H. neanderthalensis (red curves) in Europe (4° × 4° average centred around 45° N, 6° E). c, Same as b but for H. heidelbergensis (dark red) and H. sapiens (orange) in eastern-central Africa (4° × 4° average centred near 5° S, 36° E). d, Same as c but for southern Africa (4° × 4° average centred near 24° S, 24° E).The shaded curves represent probability estimates of the occurrence of respective fossil and/or archaeological data obtained from the ages and age

Fig. 4 |
Fig. 4 | Schematics of H. heidelbergensis succession.On the basis of fossil ages, we propose a split of H. heidelbergensis into northern and southern branches (blue shading, habitat suitability) around 850-650 ka.The gradual transition at 300-200 ka of H. heidelbergensis into H. sapiens in southern Africa is supported by fossil and archaeological data in this region and habitat overlap estimates (Fig.2a).The proposed divergence at 400-300 ka of H. heidelbergensis into Neanderthals in Europe is consistent with recent genetic estimates34 .This scenario is also in agreement with Neanderthal

and
Homo sapiens.The centre locations of a 5°x5° averaging area are selected as: b) Neanderthal (51°N, 7°E); c) Shkul cave (33°N, 35°E); d) Denisova Cave (51°N, 85°E).Extended Data Fig. 7 | Potential hominin overlap regions.a) Shading indicates the relative occurrence frequency (relative to the entire 2-million-year simulation history) of when the habitat suitability for Homo sapiens and Homo heidelbergensis at a grid point both exceed 0.3 and exhibit a difference of less than 0.1.This constraint is considered here as a measure of habitat similarity during overlap times.Circles show the respective fossil and archaeological sites from the overlap periods.Cyan boxes highlight areas, for which there are mixed fossils available and high habitat similarity.These areas are potential areas for speciation events or gradual evolutionary transformations.b) same as a), but for Homo neanderthalensis and Homo heidelbergensis; c) same as a) but for Homo erectus and Homo heidelbergensis; d) same as a) but for Homo erectus and Homo ergaster/habilis; e) same as a) but for Homo ergaster/habilis and Homo heidelbergensis.The size of the circles scales with the age of the fossils.Extended Data Fig. 8 | Climate effects on time-mean habitat suitability of Homo sapiens, Homo neanderthalensis, Homo heidelbergensis.a-c) difference between time-averaged habitat suitability estimated from original climate envelope model H and climate envelope model H scr calculated from vector-randomized climate variables (Methods) for 3 hominin groups.This type of time randomization (Methods) maintains the time-mean and climate co-variance among the 4-dimensional climate input vector.Results show the effect of the original orbital-scale climate trajectory (relative to a randomized climate trajectory) on hominin habitat suitability.Only grid points with p < 0.05 are shown in colors based on a paired t-test using H(x,y,t) and H scr (x,y,t).d-f) Illustration of late-Holocene altitude-corrected downscaling of net primary production: Left: Simulated NPP (gC/m 2 /year) for original T31 36  insets).The southern branch experienced considerable climatic stress in southern Africa during Marine Isotope Stages 11 and 9, which could have accelerated either a gradual or a cladogenetic transition into H. sapiens36.The Eurasian populations of the northern branch further bifurcated around 430 ka, possibly giving rise to Denisovans, which populated parts of central and eastern Asia.Inside central Europe, H. heidelbergensis, which experienced strong local climatic stress due to eccentricity-modulated ice-age cycles (Fig.