Where did homo sapiens first settle after migrating from Africa and colonising the world?

Article https://doi.org/10.1038/s41467-024-46161-7

The Persian plateau served as hub for Homo sapiens after the main out of Africa dispersal

Received: 5 August 2023 Accepted: 16 February 2024

Check for updates

Leonardo Vallini 1 , Carlo Zampieri1, Mohamed Javad Shoaee2,
Eugenio Bortolini 3, Giulia Marciani 3,4, Serena Aneli 5, Telmo Pievani1, Stefano Benazzi3, Alberto Barausse1,6, Massimo Mezzavilla 1,
Michael D. Petraglia 7,8,9 & Luca Pagani 1,10

A combination of evidence, based on genetic, fossil and archaeological find- ings, indicates that Homo sapiens spread out of Africa between ~70-60 thou- sand years ago (kya). However, it appears that once outside of Africa, human populations did not expand across all of Eurasia until ~45 kya. The geographic whereabouts of these early settlers in the timeframe between ~70-60 to 45 kya has been difficult to reconcile. Here we combine genetic evidence and palaeoecological models to infer the geographic location that acted as the Hub for our species during the early phases of colonisation of Eurasia. Leveraging on available genomic evidence we show that populations from the Persian Plateau carry an ancestry component that closely matches the population that settled the Hub outside Africa. With the paleoclimatic data available to date, we built ecological models showing that the Persian Plateau was suitable for human occupation and that it could sustain a larger population compared to other West Asian regions, strengthening this claim.

A growing body of evidence indicates that the colonisation of Eurasia by Homo sapiens was not a simple process, as fossil and archaeological findings support a model of multiple migrations Out of Africa from the late Middle Pleistocene and across the Late Pleistocene1–6. Traces of these early dispersals are also evidenced in the genome of our Nean- derthal relatives, which illustrate interbreeding events as humans moved into Eurasia7–10. Early dispersals of our species were likely accompanied by population contractions and extinctions, though succeeded by a subsequent, large-scale wave at ~70–60 kya11–15, from which all modern-day non-Africans descend16,17.

The geographically widespread and stable colonisation of Eurasia appears to have occurred at ~45kya through multiple population expansions associated with a variety of stone tool technologies18,19. Earlier incursions into Europe have been recorded20–23, however, they

failed to leave a significant contribution to later populations. A chronological gap of ~20ky between the Out of Africa migration (~70–60 kya) and the stable colonisation (~45 kya) of West and East Eurasia can be identified, for which the geographic location and genetic features of this population are poorly known. On the basis of genetic and archaeological evidence, it has been suggested that the Eurasian population that formed the first stable deme outside Africa after ~70–60 kya can be characterised as a Hub population18, from which multiple population waves emanated to colonise Eurasia, which would have had distinct chronological, genetic and cultural char- acteristics. It has also been surmised that the Hub population cannot be seen as simply the stem from which East and West Eurasians diverged. Instead, this was a more complex scenario, encompassing multiple expansions and local extinctions18. Previous studies, however,

1Department of Biology, University of Padova, Padova, Italy. 2Department of Archaeology, Max Planck Institute for Geoanthropology, Jena, Germany. 3Department of Cultural Heritage, University of Bologna, Bologna, Italy. 4Research Unit Prehistory and Anthropology, Department of Physical Sciences, Earth and Environment, University of Siena, Siena, Italy. 5Department of Public Health Sciences and Pediatrics, University of Turin, Turin, Italy. 6Department of Industrial Engineering, University of Padova, Padova, Italy. 7Human Origins Program, Smithsonian Institution, Washington, DC 20560, USA. 8School of Social Science, The University of Queensland, Brisbane, QLD, Australia. 9Australian Research Centre for Human Evolution, Griffith University, Brisbane, QLD, Australia. 10Institute of Genomics, University of Tartu, Tartu, Estonia. e-mail: leo.vallini.lv@gmail.com; luca.pagani@unipd.it

Nature Communications | (2024)15:1882

1

     

1234567890():,;

1234567890():,;

Article

https://doi.org/10.1038/s41467-024-46161-7

Fig. 1 | Relationship and legacy of the West and East Eurasian Core populations. Summary of the major population events (A) and schematic representation of our reference space and expected position of admixed and unadmixed populations (B – note that this panel is a rotation of the blue and red portion of A); derived allele sharing of each ancient or modern individual/population with Kostenki14 and Tianyuan (C); East Asians in red, Oceanians in orange, Native Americans in pink,

have failed to delve into the potential geographic location of this Hub population24, the overall scarcity of fossil evidence of Homo sapiens between 60 and 45 kya anywhere across Eurasia.

The aforementioned scenario was grounded in evidence stem- ming from ancient genomes from West and Central Eurasia25,26 and China27, indicating that the ancestors of present-day East Eurasians emerged from the Hub at ~45 kya (Fig. 1A, red branch). These emergent groups subsequently colonised most of Eurasia and Oceania, though these populations became largely extinct and were assimilated in West Eurasia28 by a more recent expansion that took place by ~38kya (Fig. 1A, blue branch). The first of these two expansions, whose asso- ciated ancestry we name here the East Eurasian Core (EEC), left des- cendants in Bacho Kiro, Tianyuan, and most present-day East Asians and Oceanians. The second expansion, which we name the West Eur- asian Core (WEC), left descendants in Kostenki14, Sunghir, and sub- sequent West Eurasians, and in the genome of palaeolithic Siberians29. Crucially, the Hub population accumulated some drift together with the WEC in the millennia that elapsed between the EEC and the WEC expansions (Fig. 1A, grey area). Despite its key role during the peopling

South Asians in yellow, Northern South Asians in green, West Eurasians in blue, Levantines in cornflowerblue, ancient samples in black (N stands for Neolithic, WHG for Western Hunter Gatherers). A graph with all populations analysed, and individual population names are in Supplementary Fig. 1, Source Data in Supple- mentary Data 3.

of Eurasia, the geographic location and the genetic characteristics of the Hub population, remain obscure24. The outlined scenario is com- plicated by the need to account for the Basal Eurasian population (Fig. 1A, green), a group30 that split from other Eurasians soon after the main Out of Africa expansion, hence also before the split between East and West Eurasians. This population was isolated from other Eurasians and later on, starting from at least ~25 kya31,32, admixed with popula- tions from the Middle East. Their ancestry was subsequently carried by the population expansions associated with the Neolithic revolution to all of West Eurasia.

Given the current impossibility of directly inferring the homeland of the Hub population from fossil remains, here we combine available genetic evidence (including both ancient and present-day genomes) and palaeoecological models to infer the geographic region that acted as a Hub for the ancestors of all present-day non-Africans during the initial colonisation of Eurasia. With our work, we show that populations from the Persian Plateau carry an ancestry component that closely matches the population that settled the Hub outside Africa, therefore pointing to the Persian Plateau as suitable for human occupation

Nature Communications | (2024)15:1882

2

Article

https://doi.org/10.1038/s41467-024-46161-7

throughout 60–40 kya, indirectly shedding light on the early interac- tions and admixture of our species with Neanderthals33 and the rela- tionships between the main Eurasian and the elusive Basal Eurasian human population30 as well as informing on where future archae- ological investigations should be focused.

Results

Rationale

The characteristics of the Hub population were outlined from a genetic perspective, including available data from extant and ancient popula- tions. In consideration of the complex relationship between the Hub population and the EEC and WEC expansions, we aimed to retrieve the genetic profile of the population (Fig. 1A, grey) that remained in the Hub location after the separation from EEC (Fig. 1A, red) which would have shared minimal drift with the WEC (Fig. 1A, blue). In light of patterns of shared drift, we hypothesized that, within a Eurasian landscape, the relic Hub population should be found among those showing a West Eurasian signature, albeit with the least shared drift within representatives of the WEC expansion.

We build a framework by describing ancient and extant popula- tions in terms of their derived allele sharing (DAS) with the two oldest unadmixed representatives of the EEC and WEC waves. The two reference samples that we considered are Tianyuan27, a 40 ky old individual from East Asia, and Kostenki1434, a 38 ky old individual from Western Russia. DAS with Tianyuan defines the red, vertical axis (in T units) in Fig. 1B, while DAS with Kostenki14 defines the blue, horizontal axis (in K units). In the absence of admixture and other confounders, most ancient and extant Eurasian genomes should fall either on the red or on the blue axes in a position proportional to the evolutionary time spent together with Tianyuan or with Konstenki14 (R and B in Fig. 1B). Under these ideal circumstances, a population falling at the intersec- tion of the two axes (black dot in Fig. 1B) would mimic the genetic features of the Hub population at the time when the EEC wave departed, while a population falling just right of this point, along the blue axis, would show the legacy of the Hub population after the EEC expansion.

There are, however, at least three confounders that complicate this search: archaic admixtures, interaction between EEC and WEC waves after they left the Hub location (the dashed, blue/red line in Fig. 1A), and admixture with Basal Eurasians (the dashed, blue/green line in Fig. 1A). The genetic contribution from archaic hominins33,35 decreases as a function of time due to purifying selection36,37 and is likely to bias the affinities between ancient samples through archaic sharing. This can be overcome by removing sites where Denisovans or Neanderthals share the same derived alleles as the other human populations considered here. On the other hand, groups that experi- enced additional gene flow from archaic humans (Oceanians and Bacho Kiro)25 will experience a reduction in their K and T coordinates. The other two confounders are illustrated in Fig. 1A, B and are the result of admixture events that occurred in the last 40 ky. On one hand, an admixture between descendants of the EEC (R, for red, in Fig. 1B) and the WEC (B, for blue, in Fig. 1B) expansions would give birth to an admixed population (BR in Fig. 1B) occupying a position intermediate between the two source populations. On the other hand, the reported admixture of certain groups with the so-called Basal Eurasian population30 starting at least from 25 kya31,32 would decrease both their T and K coordinates. Basal Eurasians are a group that is thought to have split from other Eurasians before the divergence between the EEC and the WEC and perhaps even before the establishment of the Hub population. This is, therefore, prior to the intersection of the blue and red axes in Fig. 1B, and is represented as a green dot in the bottom left part of Fig. 1B. Assuming that each analysed genome (Supplementary Data 1–2) was putatively affected by all these confounders, its resulting observed (Obs, Fig. 1B) point should therefore be brought to the B position through a number of corrections before assessing its putative

role as a Hub relic. We also note that a fourth confounder may be represented by a recent interaction with Sub-Saharan genetic components38.

DAS with Kostenki/Tianyuan

When examining the positions of modern and ancient populations (Supplementary Data 1–2, Supplementary Fig. 1A) within the reference space defined by the DAS with Kostenki14 and Tianyuan (Fig. 1C, Supplementary Fig. 1B, Supplementary Data 3), it may be observed that East Asians and Pre-Ice Age Europeans resemble the positions held by R and B in Fig. 1B. African populations align along the bisector but before the intersection of the red and blue lines representing KT (the K coordinate of Tianyuan) and TK (the T coordinate of Kostenki14), ordered by how early they diverged from the other populations (Khoe and-San are at the bottom, Supplementary Data 3). Oceanians show a greater affinity for Tianyuan compared to Kostenki14, with these populations being predominantly part of the EEC wave15,18,39. These populations are, however, shifted towards the bottom left corner of Fig. 1C by their admixture with Denisovans35, which reduces both K and T values (see Methods for more details on the role of archaic alleles in our calculations). A similar effect can be observed in the Initial Upper Palaeolithic (IUP) Bacho Kiro individuals, which have been shown to have a Neanderthal ancestor in their recent genealogy25. As unadmixed South-East Eurasians, the populations from the Andaman Islands are aligned along the red axis at a lower T than East Asians; this is com- patible with the shorter time spent together with the ancestors of Tianyuan during the EEC movement away from the Hub location, and with the earliest evidence of human occupation in South Asia at ~48–45 kya40.

The central portion of the top right quadrant of Fig. 1C is occu- pied, as expected, by groups resulting from the admixture between the EEC and the WEC waves anytime after they expanded from the Hub location. South Asians follow a cline connecting Andamanese and West Eurasians; this is expected given the reported interaction between Ancestral South Indian (ASI) and Ancestral North Indian (ANI) populations41,42. The position of Native Americans suggests a primarily East Asian ancestry, with a smaller contribution from palaeolithic West Eurasian populations43,44. Similarly, Mal’ta and Yana fall in an inter- mediate position between the two axes, the result of a palaeolithic admixture between EEC and WEC groups18. The 45 ky old individual from Ust’ Ishim26 falls close to the intersection of the red and blue axes, where we expect the 45 kya Hub population to be positioned. This location is anticipated as the lineage of this individual forms almost a trifurcation with Tianyuan and Kostenki14. Zlatý kůň lies further behind, in accordance with its lineage being unequivocally basal to the split between EEC and WEC populations18,21. West Eurasians, North Western South Asians, and Levantines occupy the area below the bisector, compatible with an admixture between EEC and WEC, or below the blue axis, further complicated by the presence of Basal Eurasian or African components in these populations. Since interac- tion with Sub-Saharan genetic components38 will have an effect similar to the interaction with Basal Eurasians, we excluded populations showing evidence of gene flow from Africa (Methods, Supplementary Fig. 2, Supplementary Data 4, 5).

Accounting for DAS confounders

Starting from the K and T coordinates estimated from each ancient or modern genome in Fig. 1C, we assume each assessed genome is potentially affected by Basal Eurasian and/or by admixture between the EEC and the WEC. We exploited the fact that the coordinates of an admixed population correspond to the mean of the two source populations weighted by each source’s contribution (we empirically validated this equation by creating synthetic mixed individuals, see Supplementary Data 6–8, Supplementary Fig. 3). For each tested population we assumed the empirically retrieved K and T coordinates

Nature Communications | (2024)15:1882

3

Article

https://doi.org/10.1038/s41467-024-46161-7

to be KObs and TObs. On the basis of this, we computed the K coordinate of the source WEC population (KB in Fig. 1B) following the scheme in Fig. 1B (see Material and Methods). Given that the goal of our analyses was to determine the population with the smallest corresponding KB coordinate (the one that could fall closest to the Red/Blue origin in Fig. 1), we assessed whether the analytical framework that we devel- oped could retrieve a correct ranking of proximity to the Hub popu- lation. Using msprime45, we performed coalescent simulations to obtain the distinct WEC source populations under different demo- graphic scenarios46,47 (Methods, Supplementary Code 1–2, Supple- mentary Fig. 4). In particular, we simulated WEC populations with different allele sharing with Kostenki14 (i.e., mimicking WEC popula- tions with different distances from the Hub population) and mixed them with the EEC and Basal Eurasians (Supplementary Data 9). We found that our approach retrieves the correct ranking along the K coordinate in the majority of cases, and with an accuracy of >0.9 in all cases where the admixed populations are at least 50% WEC, and the mixing WEC sources have at least 3 ky of differential allele sharing with Kostenki14 (Supplementary Fig. 5A). Decreasing the differences between the mixing WEC sources (hence decreasing the differences in their K coordinates) results in higher WEC fractions in the admixed population to maintain a comparable level of accuracy (Supplemen- tary Fig. 5B, C). The lowest accuracy in our simulated results is obtained when the compared populations are the result of the mixture of the WEC and EEC sources (i.e., without any contribution from Basal Eurasians), a composition that appears to be virtually absent in mod- ern Eurasian populations.

Having confirmed the validity of our approach we tested the existing data. We found that after accounting for East and Basal Eur- asian confounders, the populations that harbour the WEC component closer to the Hub population (grayscale gradient of population points in Fig. 2A, Supplementary Data 11) are the ones whose West Eurasian ancestry is related to the hunter gatherers and early farmers from Iran48. This is a genetic ancestry commonly referred to as the Iran Neolithic30 or the East Meta49, here named Iran HG for clarity (Sup- plementary Data 11). The Iran HG ancestry is widespread not only in modern-day Iran but also across ancient and modern samples from the Caucasus (in particular in the Mesolithic hunter gatherers of that region) and in the northwestern part of South Asia50. Along the blue axis of genetic similarity to Kostenki14, these populations come before modern and ancient groups from the Levant and, in turn, before groups from Europe and other areas associated with the Anatolian

Neolithic expansion49,51–53. The furthermost groups along this axis are post- and pre-LGM European hunter gatherers, which is expected owing to their genetic proximity to Kostenki14.

We report the affinity to the Hub population of each analysed group showing at least a 75% WEC genetic fraction (the strictest threshold we identified in our validation), after correcting for the aforementioned confounders, in grayscale on the points on a geo- graphic map (Fig. 2A). Assuming at least partial population continuity in the past 40 ky, we identify focal areas54 for the Hub. These locations are defined as parts of the map where the available West Eurasian genetic component is closest to the Hub population and where such genetic proximity decreases as a function of geographic distance from the said location (Fig. 2A). For each location of the map, with at least one sampled population within 2000 km, we computed the shortest overland distance between all other sampled populations, and we then estimated the Pearson correlation with the corresponding values of −KB. A negative r value (light colours in Fig. 2A) is expected in the locations from where the K coordinates can only increase. The focal area for the Hub location, the lightest shade in Fig. 2A, falls in the region between the southern shores of the Caspian Sea and the Persian Gulf (then not submerged). The results are qualitatively the same when changing the inclusion criteria in terms of the fraction of the West Eurasian genetic component in the admixed individuals (Supplemen- tary Fig. 6A) or the age of the samples (Supplementary Fig. 6B). We performed the same focal area analysis with the inferred proportion of Basal Eurasian (Fig. 2B, Supplementary Fig. 6C, D), and found that the most likely entry point for such an ancestry is separated from the putative Hub location. This appears to be linked to the Levant, sug- gesting either that the Levant, the Arabian Peninsula or North Africa was a potential location for this elusive population.

In integrating the genetic results within a spatially explicit model, should be noted that post Neolithic expansions might have con- tributed to the spread of a Hub-like component beyond its homeland; for example, towards northern South Asia, along with the expansion of the so-called Iranian Neolithic genetic components48,50. In addition, other population movements might have diluted its presence in the Hub location with the arrival of other WEC components with a lower Hub affinity (e.g., via the Eastward spread of Anatolian Neolithic components)50. Therefore, putative legacies of the Hub may be found over a large area, stretching from the Southern Caucasus to northern South Asia, though this may not have always been the case. In the Caucasus, pre-LGM hunter gatherers were more closely related to early

Fig. 2 | Focal areas for the Hub and Basal Eurasian genetic components. Focal area for the Hub (A, white to light yellow hues show the most likely Hub location from a genetic perspective) and for Basal Eurasian ancestry (B, dark hues show higher proportion of a Basal Eurasian component) based on at least 75% WEC

ancestry. The grayscale in individual population points is proportional to the inferred proximity to the Hub (A) or the proportion of Basal Eurasian ancestry (B). kya stands for thousands of years ago. Source Data in Supplementary Data 11.

Nature Communications | (2024)15:1882

4

Article

https://doi.org/10.1038/s41467-024-46161-7

agriculturalists from western Anatolia31,32 than to the Mesolithic hunter gatherers (CHG, carrying an ancestry strictly related to Iran HGs). This suggests an expansion of populations from the Hub population to the Caucasus between 25 and 13 kya. This would, therefore exclude the Caucasus as a location for the Hub unless a more complex scenario, such as a double population replacement, is postulated. The presence of a Western Eurasian component in northern South Asia has tradi- tionally been explained as the result of the eastward expansion of Iranian farmers48. A recent study, however, reported the presence of this ancestry in a ~4500 year old sample from the Indus Valley, and inferred that it split from Iranian farmers before the advent of agri- culture, suggesting that the WEC genetic component may predate the Iranian Neolithic expansion55. Nevertheless, as the case of the Caucasus has shown, genetic continuity before the advent of agriculture might not necessarily mean that it dates back to the timeframe of interest. While we can not exclude it, a long term presence of a population Hub in South Asia is at odds with the existence of an indisputably EEC genetic component referred to as ASI (or AASI) that made up the majority of the pre-Neolithic genetic landscape50.

We inferred the focal area from which Basal Eurasian ancestry expanded into Eurasia, an event reported to have taken place before the Neolithic, but after the WEC expansion. Figure 2B points to northern Africa56, the Arabian peninsula or the Levant as the most likely homeland for Basal Eurasians, and its expansion across West Asia follows a gradient of radial distance from such areas. Given the inferred focal area for Basal Eurasians, and the notion that such a population was genetically separated from the Hub population until at least 38 kya, it follows that the Hub location must have been physically separated from the Basal Eurasian homeland. For this reason, the putative Hub location should be geographically distinct from the location of Basal Eurasians, a criterion that is met in the Persian Plateau.

Palaeoecological modelling

The emerging picture for the putative Hub location was based on genetic data alone and reliant on assumptions of population con- tinuity, or at least the absence of major population replacements across the area. To test the plausibility of our continuity assumption, we utilised palaeoecological models as additional and independent layers of information to assess whether such a geographic region was suitable for human occupation at any given time. Recently, palaeocli- matic reconstructions that include the time interval between 70 and 30 thousand years ago have been published57,58. We took advantage of this data to build a species distribution model following the procedure described by Rodríguez and colleagues59 (see Methods, Supplemen- tary Figs. 7–11) and to reconstruct the areas with environmental con- ditions suitable for human occupation throughout that period (Fig. 3A, Supplementary Fig. 12). We then combined these results with the estimates of Net Primary Productivity (NPP) for each region and the relation between NPP and hunter gatherer population density59. This was used to estimate the maximum sustainable human population (carrying capacity) through time in different geographic regions (Fig. 3B, geographic regions shown in Supplementary Fig. 13).

Our palaeoecological model shows that the putative Hub location identified in Fig. 2A (and framed in Fig. 3A) could have supported human occupation throughout most of the time periods between 70 and 30 kya. Our model also shows that the region could have poten- tially sustained a population size much higher than other Western Asian regions during the same time intervals. This inference can be made even without taking into account the presence of Mesopotamian rivers and the rich hydrological network of the Persian Plateau60, which likely extended and interconnected areas potentially inhabitable beyond the limits of those deemed suitable by our model. It is inter- esting to note that the Mesopotamia/Persian Plateau becomes patchier in terms of habitat suitability between 60 and 50 kya (Supplementary Fig. 12); and then reconnects again after 50 kya, a trend also reflected

Fig. 3 | Palaeoecological inference of areas suitable for human occupation.

Number of time intervals between 70 and 30 kya (thousands of years ago, one every 2000 years, 21 in total) in which each area is predicted as suitable for human occupation by our model (A). Maximum sustainable human population in each geographic region over time (B). The geographic regions are shown in Supple- mentary Fig. 13. The black frame in A shows the area predicted to be the Hub location from genetic evidence. Source Data in Supplementary Data 14.

by the increase in carrying capacity of the region (Fig. 3B). It is possible that this may have provided an ecological trigger for the IUP expansion that occurred sometime before 45 kya18,61.

Due to the scant availability of direct palaeoclimate records from the Persian Plateau, the data we used for our palaeoecological mod- elling has been validated locally on a single datapoint from Iran. In spite of this and thanks to the higher resolution of the validating data provided by other areas with similar environmental conditions, the habitat inferred by our model can be considered reliable as indepen- dently validated by other approaches. In recent work, researchers utilised paleohydrological mapping along with existing paleoenviron- mental proxies and paleoclimate modelling to reconstruct the histor- ical climate of the Persian Plateau, aiming to explore climatically influenced pathways for hominin dispersals during MIS 5 and MIS 360. The paleoclimate model for MIS 3 indicates substantial increases in moisture for both the Zagros Mountains (70–30 kya) and the northern Persian Plateau (50–40kya) during specific periods within MIS 3. Notably, these conditions could have potentially supported hominin habitation in these areas, a feature also picked up by our model. This aligns with the findings of available proxy records and the spatial dis- tribution of archaeological sites60,62–67. When combining this evidence, it appears clear how the ecological variability of the Persian Plateau picked up in our model corresponds with a period of increased aridity during MIS 4 and a later recovery at the onset of MIS 3, which brought more favourable environmental conditions, though not as ideal as those of MIS 560. It is relevant to note that even at its lowest carrying

 

Nature Communications | (2024)15:1882

5

Article

https://doi.org/10.1038/s41467-024-46161-7

Fig. 4 | Combination of palaeoecological and genetic analyses. In light yellow, within the black frame, are geographic locations that are putative Hub focal areas and predicted habitable areas. The areas are compiled on the basis of at least 90% of

capacity prediction, the value of the Persian Plateau is almost equal to the highest values in other regions of the Middle East. The Persian Plateau is also higher in comparison to other periods, providing a clue for its competitive advantage over surrounding areas to serve as a Hub location. The palaeoclimatic results, therefore, confirm the suitability of the Persian Plateau for human occupation, validating and refining the picture emerging from the genetic data (Fig. 4, light yellow area). Furthermore, the presence of a viable area located on both shores of the Red Sea and stretching across the Mediterranean Sea would seem to offer a suitable habitat for the Basal Eurasian populations, partly disconnected from the putative Hub location (Fig. 4, green area).

Fossil and archaeological evidence

The hypothesised links of the EEC wave with the IUP, and of the WEC wave with the Upper Palaeolithic (UP) are informative about the expected material culture to be found in the putative Hub location18. According to this scenario, IUP and UP technologies are expected to appear after ~45 kya and ~40 kya, respectively. On the other hand, if the Hub location hosted the population from which even earlier, failed attempts to colonise Eurasia stemmed, such as the assemblage asso- ciated with the Zlaty Kun expansion18,21, cultures perhaps more basal than IUP and UP are also to be expected in the area. With this in mind, and in light of the scarcity of hominin remains in the Persian Plateau dated before 40 kya, a brief examination of the fossil and archae- ological record is warranted.

Archaeologists working in the Persian Plateau have begun to pay considerable attention to hominin dispersals in the Late Pleistocene60,68,69. A recent climatic and palaeohydrological study has indicated that the region experienced ameliorated conditions in MIS 5 (130–71 kya) and MIS 3 (57–29 kya), consistent with the widespread location of Middle Palaeolithic (MP) sites in the Zagros Mountains and in lowland riverine and lake settings. MP sites in the Zagros Mountains, dating to between ~77–40 kya, are sometimes in association with rare Neanderthal fossils68,70. No fossils of Homo sapiens have yet been recovered in association with the MP in the Persian Plateau, though associations are found in Arabia at ~85 kya3 and in the Levant in MIS 571 and at ~55kyr72. The latter finding from Manot Cave is dated to 55 ± 5 kya based on the calcitic patina that covered the bone, hence indicating a minimum age3. If the actual date was older this individual

the time intervals inspected by our palaeoclimatic analysis or those located along major rivers. In green are the habitable areas that might have hosted the Basal Eurasian population.

would be connected with the former sites, which are too old to represent the ancestors of modern-day non-Africans and are instead associated with an older dispersal OoA. On the other hand, if the dating of Manot reflects its actual age (hence around the time of the OoA), this finding is not at odds with our inference of the Hub being in the Persian Plateau, as the site lies on one of the possible routes leading to it from Africa13 If this fossil was instead to be considered a member of the Hub population, based on the lack of direct genetic information and given the lack of support for subsequent Levantine individuals to be direct descendants of the Hub population, one needs to postulate a later complete replacement of the hypothetical Levantine Hub population with a less basal lineage, which we here deem to be a non-parsimonious scenario. Archaeologists working in Iran have considered that MP Levallois technologies may also be the product of modern humans60,73. And in fact, the MP assemblages of Mirak, dated to ~55 kya, and situ- ated on the Dasht-I Kavir ancient lakes and wetland systems in the Iranian Central Plateau60,74 have been suggested to represent Homo sapiens. This assertion is based on a 3D geometric-morphometric comparison of lithic assemblages, which show the consistency of MP and UP tool reduction strategies at Mirak and their similarity with Levantine assemblages manufactured by modern humans60,75.

Standardised blade and bladelets industries of the Persian Plateau have been referred to as the Baradostian, the Zagros Aurignancian, and the IUP68,76,77. These UP techno-complexes appear to make a relatively rapid appearance across the Persian Plateau78, in support of a wide- spread replacement of earlier populations which may have included the Neanderthals and earlier Homo sapiens populations68.

Discussion

Here we characterised the genetic and geographic features of the Eurasian population Hub, where the ancestors of all present-day non- Africans lived between the early phases of the OoA expansion (~70–60 kya) and the broader colonisation of Eurasia (~45 kya). Our study was not designed to reconstruct earlier population expansions from the Hub20,21 as well as previous Out of Africa events, which may have left traces in the fossil and archaeological records but contributed little or none of the genetic makeup of present-day individuals79.

Our results showed that the genetic component closest to the Hub population is represented in ancient and modern populations in

Nature Communications | (2024)15:1882

6

Article

https://doi.org/10.1038/s41467-024-46161-7

the Persian Plateau. Such a component, after mixing with Basal and East Eurasian ancestries, resurfaced in the palaeogenetic record, pre- viously referred to as the Iranian Neolithic, the Iranian Hunter Gath- erer’ or the East Meta49.

This genetic perspective was further refined when combined with palaeoecological evidence, pointing to the Persian Plateau belt sur- rounding the Central Iranian Desert, and encompassing the Southern Caspian shores, the Zagros Mountains, and the Persian Gulf and Mesopotamia, as the most likely Hub location (Fig. 4). Although no genetic data is currently available, we cannot exclude the northeastern portion of the Arabian peninsula as the southern edge of the Hub80 since the area was putatively linked across the then drier Persian Gulf. The ecological separation offered by the harsh climatic conditions of the inner parts of present-day southern Levant and western Arabia seems to provide an explanation for the long-lasting disconnect between the Hub population and Basal Eurasians30, which enabled the two human groups to build their characteristic drift components over many thousand years (their potential locations are shown in green and light yellow colour in Fig. 4).

Information from archaeological evidence advocates for a long presence of modern humans in the region, dated to at least 44 kya when considering the earliest IUP and UP record68,76,77 or to as early as 55 kya when including certain MP sites that may relate to modern humans60,75. Furthermore, the attested presence of Neanderthals in the Zagros mountains until at least 40 kya81, supports the idea that the ancestors of all living Eurasians spent ~20 ky in the Hub location and admixed with our archaic relatives, potentially over a prolonged per- iod of time82, ultimately differentiating into the populations that eventually led to the colonisation of Eurasia, Oceania, and the Amer- icas. At around 45 kya, a population expansion associated with the IUP61,83 emanated from the Hub18, spreading the EEC ancestry and colonising most of Eurasia25,84 This expansion left descendants only among present-day East Eurasians27 and Oceanians, but largely faded and was assimilated in Europe28 by a subsequent expansion, carried out by WEC groups, taking place ~38 kya18.

While earlier incursions into Europe have been documented20,21, they left no trace in the gene pool of present-day and ancient populations. The reasons behind this significant temporal gap between the main OoA event and the first durable expansions deeper into Eurasia can only be hypothesised; but the early stages of the colonisation of Eurasia surely provided with several challenges. A certain amount of time would have been necessary to demo- graphically recover from the OoA bottleneck; and new environ- mental stressors had to be dealt with either with biological adaptation (see Supplementary Note 1; Supplementary Fig. 14, Sup- plementary Data)85 or by developing technological innovations, as testified by the archaeological record. Finally, ecological stressors, exemplified by changes in habitat connectivity and the interaction with local archaic groups might have also played an important role. To this extent, the time elapsed as a single Hub population might have served as an incubator for the development of cultural inno- vations later observed almost simultaneously at opposite ends of the Pleistocene Eurasian world. Among these features the most notable is the presence of rock art at ~40 kya in Sulawesi, Indonesia86, com- patible in age with the oldest European art (41–35 kya)87,88; as well as the innovative usage of projectile weapons, recorded in Europe89–91, the Levant92, and South Asia93.

In conclusion, our multidisciplinary effort shed light on the millennia that separated the Out of Africa expansion and the differ- entiation of Eurasians into Europeans, East Asians, and Oceanians. Our study also pointed to the Persian Plateau as the most likely candidate for a Hub location of Homo sapiens populations, in con- sideration of genetic, paleontological, and archaeological evidence, thereby illustrating that this is a key region for future archaeological investigations.

Methods

Dataset

We downloaded individual genotypes13,21,25–27,29,30,43,48,50,94–127 for 1,233,013 sites from Allen Ancient_DNA_resource_(version_v.50.0_ https://reich.hms.harvard.edu/allen-ancient-dna-resource-aadr- downloadable-genotypes-present-day-and-ancient-dna-data) and con- verted it dataset to plink format128 using admixtools129 convertf; we then merged it with the Zlatý kůň21 and Bacho Kiro25 genomes (pro- cessed from bam file), Gumuz13 genomes and with individuals sequenced in94. We selected good quality aDNA samples relevant to our analyses and included all modern-day populations, limiting the number of individuals to 20 (randomly selected) if more than that was available. The full set of individuals analysed is reported in Supple- mentary Data 1, and a per-population summary is in Supplemen- tary Data 2.

We set the reference allele to the ancestral identified in the 1000 genomes project phase 3 release123 using plink --reference-allele and kept only autosomes while removing the alleles for which information on the ancestral allele was missing. We also challenged our approach by artificially swapping a fraction of the ancestral/derived allele assignments and found that results are robust even when 50% of alleles are swapped. We used a custom script to identify the derived allele of archaic humans (Neanderthals: Vindija33, Altai110, and Chagyrskaya111; Denisovans: Altai Denisova35 and the genome of a first-generation hybrid between the two species113) and removed them from the data- set. This left us with 833070 SNPs.

Derived allele sharing

For each individual or population, using a custom script, we computed the DAS, which is the proportion of derived alleles that each individual shares with a given reference sample) with Kostenki1434, Tianyuan27, and Ust’Ishim26, only considering SNPs with missing call rates not exceeding 25% (using plink --geno 0.25 option).

Removing individuals with African admixture

To identify the populations showing admixture with Africans we ran an unsupervised admixture130 analysis on the modern samples of our dataset from k = 2 to k = 10, choose k = 7 (Supplementary Fig. 2, Sup- plementary Data 4) since it is the one with the lowest cross-validation value and then excluded all populations that show a proportion higher than 0.05 of the ancestry that is maximised in Sub-Saharan African populations (blue in Supplementary Fig. 2). With this step, on top of African populations, we excluded BedouinA and BedouinB, Palesti- nians, Jordanians, and Yemenite_Yew.

To test whether any of the ancient populations had a contribution from either East or West Africans, we confronted them with French (that were shown by the previous analysis to not have African admix- ture) by computing the test D(X, French, Yoruba or Gumuz, Chimp) and excluded those with a significantly positive value (Supplemen- tary Data 5).

Artificially admixed populations
In order to create the chimeric genomes, we first selected three source populations and four individuals from each; we chose Han, Sunghir, and Gumuz as these populations somewhat mimic the R, B, and BEA points in the plot of Fig. 1A. The individuals selected are shown in Supplementary Data 6. We then extracted individual chromosomes for each and merged them as shown in Supplementary Data 7, then computed the proportion of ancestry from each source based on the number of SNPs each chromosome has.

We finally computed the DAS with Kostenki14 and Tianyuan for each source and admixed population (Supplementary Fig. 3, Supple- mentary Data 8) both empirically (on chimeric genomes) and analyti- cally (by averaging the K and T values of the source populations weighted by their contribution) and compared them. The ratio of the

Nature Communications | (2024)15:1882

7

Article

https://doi.org/10.1038/s41467-024-46161-7

two values never deviated from 1 more than 0.008 and very often much less than that (Supplementary Data 8).

Identification of the WEC source population for each sampled population
From this starting point, we made some reasonable assumptions to retrieve the position of the B source population (Fig. 1A) for each Obs (Fig. 1A) population.

The first step is determining the Basal Eurasian contribution in each population. Since Basal Eurasians diverged before the split between EEC and WEC30, and Ust Ishim forms a near trifurcation with the two above-mentioned groups, all Eurasian populations should have a similar DAS with Ust Ishim (U), unless they received some contribution from a more basal population such as Basal Eurasians or an African population (populations displaying a sizeable African con- tribution, however, were identified and removed from the analyses, see Methods, Supplementary Fig. 2, Supplementary Data 4, 5). We chose as baseline U the mean (Umean) of the U of Kostenki14 (UK) and Tianyuan (UT), and since Basal Eurasians have so far not been sampled, we used the U, K, and T of Gumuz as a conservative proxy for it in Eq. (1) and (2). Gumuz is an East African population that split from Eurasians just before the OOA and, unlike most other East Africans, lack any contribution from Eurasian populations13. In this step, if the proportion of Basal Eurasian results is negative, which can be the case for some populations that have a U slightly higher than Umean, we set it to zero.

genetic component. We then used the coordinates of Han as source if the test was significantly negative, the ones of Onge if it was positive and their average if it did not significantly deviate from zero (Supple- mentary Data 10); we did not use Irula_S and their coordinates directly in our DAS analyses because these individuals were genotyped on a different chip with fewer SNPs than the 1240 K; however ASI has been shown to be a trifurcation with East Asian and Andamanese populations42. With this final piece of information, we can compute the proportion of EEC admixture in the hypothetical BRX point (note that this is different from the proportion in population X if admixture with Basal Eurasians is also present). In this step, if pEECBRx results are negative or higher than 1, we set it to zero and one, respectively.




pEECBR = TBR 
 TWEC = TEEC 
 TWEC ð3Þ

Where:
pEECBR = the proportion of EEC ancestry in the hypothetical BR

population
TWEC = T coordinate of the source WEC population (B in Fig. 1A) of

the BR point
TEEC = T coordinate of the source EEC population (R in Fig. 1A) of

the BR point
Finally, we compute the coordinates of the unadmixed WEC

source (B in Fig. 1A) of each admixed population. 



KB = KBR 
KEEC pEECBR = 1
pEECBR ð4Þ Additionally, to get the total proportion of EEC and WEC ancestry




pBEAX = UX 
 Umean = UBEA 
 Umean

ð1Þ

Where:
pBEAX = proportion of Basal Eurasian ancestry in population X UX = derived allele sharing with Ust’ Ishim of population X
UBEA = derived allele sharing with Ust’ Ishim of the Basal Eurasian

population
Umean = average of the U of Kostenki14 and Tianyuan
Having estimated the proportion of Basal Eurasian we can then

find the position a population would occupy if it did not have it (BR in Fig. 1A).

in population X, one can compute: pEECX=pEECBR 1
pBEAX







KBR = KX 
KBEA pBEAX = 1
pBEAX 



TBR = TX 
TBEA pBEAX = 1
pBEAX

ð2:1Þ ð2:2Þ

Where:
pWECX = the proportion of WEC ancestry population X

Coalescent simulations

Where:
KX, TX = DAS with Kostenki14, Tianyuan of population X
KBEA, TBEA = DAS with Kostenki14 and Tianyuan of the Basal Eur-

asian population
KBR, TBR = K and T coordinates of the hypothetical BR population

corresponding to population X, i.e., the coordinates it would have if it lacked Basal Eurasian ancestry

We can then determine whether population X has any contribu- tion from an EEC population, knowing that if this was 0 then its TBR coordinate should be equal to TWEC; if not, then we need to know the T coordinates of the two sources (B and R in Fig. 1A) in order to deter- mine the proportion of EEC ancestry. We approximate TWEC with the average of the T coordinates of the available WEC populations older than 30 ky (Sunghir116, Vestonice1696, Muierii2109, and Krems95).

For TEEC, we first need to identify the source EEC population (R in Fig. 1A). Unadmixed EEC populations can be grouped into East Asians (red in Fig. 1B and with higher T values) and Andamanese/ASI (brown in Fig. 1B, with lower T values); to determine which of these two groups was contributing to each population (or if they both did contribute), we used the test D(Irula_S, Han, X, Mbuti). Irula_S are the masked ASI genomes obtained through ancestry deconvolution by Yelmen and colleagues42, hence a good proxy for an unadmixed South Asian

We simulate the OOA 60 kya, with Basal Eurasians (BEA in Supple- mentary Fig. 4) splitting soon after (57.5 kya) and the split between EEC and WEC, with the former leaving the Hub18, 46 kya (allowing the time for them to reach Ust’Ishim and Bacho Kiro by ~45 kya). Both Tianyuan and Kostenki14 split off from their respective main branches (EEC and WEC) 34 generations (~1000 years) before their death. We simulated two different West Eurasian populations: WEC and WEC2, with WEC2 staying in the Hub longer than WEC (and Kostenki14), and hence closer to it from a genetic point of view. We then have each of these populations acting as a source for admixture events with Basal Eur- asians (BEA) and East Eurasians in different proportions (Supplemen- tary Data 9). We simulated 10 genome chunks for each population and averaged the DAS values obtained from the different chunks of KOS, TIA, and UST populations and then applied the procedure, assump- tions, and approximations included, described in the previous para- graph, to identify the B source for each admixed population where WEC is the source and compared it with all other admixed populations where WEC2 is the source and vice versa to see whether, regardless of the ancestry composition of a population, we could retrieve the cor- rect ranking in terms of vicinity to the Hub. We repeated the process in three different scenarios where WEC2 (the source closer to the Hub) shares 3000, 2000, and 1000 years less than WEC with Kostenki14 (i.e., either WEC and WEC2 left the Hub at the same time, but WEC2 bran- ched off earlier or WEC2 remained in the Hub when WEC and

 1 
 pBEAX pEECX = the proportion of EEC ancestry population X

ð5:1Þ ð5:2Þ

pWECX =

1 
 pEECBR



Nature Communications | (2024)15:1882

8

Article

https://doi.org/10.1038/s41467-024-46161-7

Kostenki14 left), and repeated the process over 500 simulations for each of the three scenarios.

An in-house python3 (reported in Supplementary Code 1) script was used to simulate genomes under the demographic scenario shown in Supplementary Fig. 4A, and the ancestry composition of the admixed populations is shown in Supplementary Data 9. Simulations rely on the msprime library45 and Ne parameters are taken from the model of Gravel and colleagues46. The simulated genomic chunks are written as a vcf file, which is then converted to plink128. When com- paring the relative position of the inferred sources of WEC ancestry in each admixed population we assigned value 1 if the ranking was correct and 0 if it wasn’t, then averaged the matrix over 500 simulations. The resulting matrices are shown in Supplementary Fig. 5A–F.

To test whether our method remained robust under different demographic scenarios, we tested its validity again under the topology proposed by Kamm and colleagues47, but we slightly modified it by incorporating an initial bottleneck when WEC and EEC populations diverge/leave the Hub and a subsequent exponential growth post Out of Africa with parameters inspired by those inferred by Gravel and colleagues46; the model is graphically represented in Supplementary Fig. 4B and the code to generate the simulated genomes is reported in Supplementary Code 2. Finally, since we relied heavily on aDNA data that is largely capture based, we assessed the impact of ascertainment bias using our simulated data: we repeated the analysis for the simu- lations of the scenario analysed in Supplementary Fig. 5C, restricting the analysis to the SNPs with an minor allele frequency higher than 0.05 in modern Eurasian populations, since the ascertainment process is more likely to identify common variants in largely studied popula- tions. The results do not differ significantly (wilcoxon test p = 0.9513, mantel statistic r = 0.9998 –p = 0.001 | maximum absolute deviation between two corresponding cells in the two matrices = 0.016); the results are shown in Supplementary Fig. 5G.

Palaeoecological model

Following the method described in Rodriguez et al.59 we built a species distribution model (SDM) with presence and background data based on the MaxEnt algorithm131 to identify environmentally suitable areas for Homo sapiens in Eurasia (10°E to 140°W; 0° to 75°N) during the period 70 kya to 30 kya.

We used the georeferenced database of hominin remains and artefacts used in Timmermann et al.132, which is a derived version of the one published in Raia et al.133, to fit the SDM. We selected only entries referring to strong evidence (i.e., tier 1, single date in their annotation) of the presence of Homo sapiens, dated within the target period (70–30 kya) and area, with dating uncertainty (i.e., difference between maximum and minimum age) lower than 4500 ky.

Estimates for 18 bioclimatic and topographic variables (including temperature and precipitation indices, and elevation (Supplementary Data 12) in the target period with high spatial (0.5° × 0.5°) and temporal (2kya timesteps) resolution, for the entire area of interest were accessed and extracted using the R package pastclim57,134. To avoid biases in fitting the species distribution model using redundant pre- dictors, the 18 bioclimatic and topographic variables were tested for collinearity and multicollinearity on the whole studied area (883260 spatial points). Pairwise Pearson’s correlation coefficients (r; Supplementary Fig. 7) were calculated using the R package Hmisc135, and the predictors strongly correlated with others (|r| > 0.9) were excluded from further modelling. The more general predictor was retained when strongly correlated ones occurred (e.g., Annual pre- cipitation over Precipitation of wettest month). Moreover, multi- collinearity R2, Tolerance, and Variance Inflation Factor (VIF), all three describing the linear dependence of one predictor on multiple others, were computed using the R package fuzzySim136. We iteratively excluded the variable with the highest VIF until all remaining variables matched the condition VIF < 5 as done by Rodriguez and colleagues59

(Supplementary Data 13). The candidate variables selected through this process were elevation, minimum annual temperature (BIO5), temperature annual range (BIO7), mean temperature of the wettest trimester (BIO8), precipitation seasonality (BIO15), and precipitation of the driest (BIO17), of the warmest (BIO18) and of the coldest (BIO19) trimesters.

Palaeoclimatic data and georeferenced observations were trans- formed to the Sphere Two-Point Equidistant projection (ESRI:53031).

Given the unbalanced distribution of sampling effort in the target area, we tested the effect of both sampling size (expressed as a number of considered spatial points) and sampling strategy in selecting the background points required by the MaxEnt algorithm. Absolute sam- ple error at six sample sizes (n = 1, 10, 100, 1000, 10000, 100000) was computed for each of the selected predictors, and the procedure was repeated 1000 times at each sample size. Absolute sample error dis- tribution showed to be stable when n > 1000, so we decided to use n = 10,000 (Supplementary Fig. 8a). The effect of (i) uniformly random sampling over all the time periods, (ii) uniformly random sampling stratified by time period (i.e., weighted by the number of archae- ological sites per time period), and (iii) effort-weighted random sam- pling was tested by comparing predictors distributions at the background points (Supplementary Fig. 8b). For this test we used a sample size of n = 10,000. The sampling effort used to weigh back- ground sampling for the third tested strategy was estimated on the collective records of Homo sapiens from Timmermann et al.132 dataset. We selected only records of tier 1 with single date in their annotation with average age smaller than 70 ky. Record points were converted into a continuous surface using a 2D kernel density estimation through the R package ks137. We then used the resulting probability surface to randomly sample background points. No considerable differences can be observed between uniform and stratified strategies, while an effect can be observed when using the effort-weighted strategy. Then, we decided to use the effort-weighted strategy to sample background points and avoid biases due to the non-uniform distribution of pre- sence samples over the study area138.

A total of 54 different configurations of the SDM resulted from the combination of six enhancing functions applied to the feature space, namely L, LQ, LQH, LQP, LQHP, LQHPT (with L = linear, Q = quadratic, H = hinge, P = product, and T = threshold), and nine regularisation multipliers (i.e., 0.2, 0.4, 0.6, 0.8, 1.0, 1.5, 2, 3, and 4). All the SDM configurations were calibrated, and their performance was measured, using the same initial covariates, sample data, and background points.

All SDM configurations were fitted through cross-validation, and then their performance was evaluated considering three indicators: the relative Akaike Information Criterion corrected for small sample sizes (ΔAICc), the average Emission Rate of the 10%-percentile of presence points (OR.10), and the Area Under receiver-operator-Curve (AUC) of the test subsample. None of the SDM configurations tested performed coherently in all three indicators. Models enhanced only by the linear function (L) outperformed in terms of omission rate but performed poorly in ΔAICc and AUC. Opposite results were obtained for the models enhanced by the most complex set of features (LQHPT). Thus, the model with a feature space enhanced by linear, quadratic, and product functions (LQP), and with a regularisation multiplier of 0.2 was selected for further analysis since it represented the best balance among the three indicators (ΔAICc = 978, OR.10 = 0.188, AUC = 0.665; Supplementary Fig. 9).

Within the eight selected predictors, temperature annual range (BIO7), mean temperature of the wettest trimester (BIO8), precipita- tion seasonality (BIO15), and minimum annual temperature (BIO5) showed to be the most informative ones in estimating human niche (Supplementary Fig. 10).

Then, we predicted habitat suitability for H. sapiens in the target area and time periods (Supplementary Fig. 11), producing continuous probability surfaces which were then binarised based on a threshold of

Nature Communications | (2024)15:1882

9

Article

https://doi.org/10.1038/s41467-024-46161-7

suitability equal to 0.238 corresponding to the 5%-percentile of pre- dicted values at observed sites (Supplementary Fig. 12).

Finally, for each geographic region (boundaries shown in Sup- plementary Fig. 13), we computed the maximum sustainable human population using the estimates of NPP and the regression equations linking NPP to Hunters population density published by Rodriguez and colleagues59 (Fig. 3B, Supplementary Data 14).

Reporting summary

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

Data availability

No new data was produced for this study. Genomic data for modern and ancient individuals is available from the Allen-Ancient DNA Resource_(https://reich.hms.harvard.edu/allen-ancient-dna-resource- aadr-downloadable-genotypes-present-day-and-ancient-dna-data). Gumuz genomic data is available from the European Genome- Phenome Archive under accession number EGAS00001000238. Alti- bathymetric geographic maps of the ETOPO1 dataset139,140 and the inferred sea levels in the past from the “Global 1 Ma Temperature, Sea Level, and Ice Volume Reconstructions” dataset141 are available from the NOAA Paleoclimatology Programme. Paleoclimatic Data57 can be accessed using the R package Pastclim58. The data to plot Fig. 1C is reported in Supplementary Data 3. The data to plot Fig. 2 and Sup- plementary Fig. 6 is reported in Supplementary Data 11 (the coordi- nates of the populations are reported in Supplementary Data 2). The data to plot Fig. 3B is reported in Supplementary Data 14. Figure 4 combines data from Figs. 2 and 3A, the rivers position was obtained from Natural Earth (free vector and raster map data natur- alearthdata.com). The map in Supplementary Fig. 1A is made with Natural Earth, the coordinates of the samples shown are reported in Supplementary Data 2. The data to plot Supplementary Fig. 1B is reported in Supplementary Data 3. The data to plot Supplementary Fig. 2 is reported in Supplementary Data 4. The data to plot Supple- mentary Fig. 3 is reported in Supplementary Data 8. The data to plot Supplementary Fig. 14 is reported in Supplementary Data 15. Supple- mentary Fig. 4A, B can be plotted using the code in Supplementary Codes 1 and 2. The code to plot Fig. 3 and Supplementary Figs. 8–12 is available142.

Code availability

Supplementary Codes 1 and 2 (the code to run the coalescent simu- lations) and the custom script used to identify the derived variants in each individual, as well as the code to carry out palaeoecological modeling available on Zenodo (https://zenodo.org/records/ 10571649)142.

References

1. Grün, R. et al. U-series and ESR analyses of bones and teeth relating to the human burials from Skhul. J. Hum. Evol. 49, 316–334 (2005).

2. Hershkovitz, I. et al. The earliest modern humans outside Africa. Science 359, 456–459 (2018).

3. Groucutt, H. S. et al. Homo sapiens in Arabia by 85,000 years ago. Nat. Ecol. Evol. 2, 800–809 (2018).

4. Harvati, K. et al. Apidima Cave fossils provide earliest evidence of Homo sapiens in Eurasia. Nature 571, 500–504 (2019).

5. Groucutt, H. S. et al. Multiple hominin dispersals into Southwest Asia over the past 400,000 years. Nature 597, 376–380 (2021).

6. Freidline, S. E. et al. Early presence of Homo sapiens in Southeast Asia by 86-68 kyr at Tam Pà Ling, Northern Laos. Nat. Commun. 14, 3193 (2023).

7. Kuhlwilm, M. et al. Ancient gene flow from early modern humans into Eastern Neanderthals. Nature 530, 429–433 (2016).

8. Posth, C. et al. Deeply divergent archaic mitochondrial genome provides lower time boundary for African gene flow into Nean- derthals. Nat. Commun. 8, 16046 (2017).

9. Petr, M. et al. The evolutionary history of Neanderthal and Deni- sovan Y chromosomes. Science 369, 1653–1656 (2020).

10. Peyrégne, S., Kelso, J., Peter, B. M. & Pääbo, S. The evolutionary history of human spindle genes includes back-and-forth gene flow with Neandertals. Elife 11, e75464 (2022).

11. Soares, P. et al. The Expansion of mtDNA Haplogroup L3 within and out of Africa. Mol. Biol. Evol. 29, 915–927 (2012).

12. Schiffels, S. & Durbin, R. Inferring human population size and separation history from multiple genome sequences. Nat. Genet. 46, 919–925 (2014).

13. Pagani, L. et al. Tracing the route of modern humans out of Africa by using 225 human genome sequences from Ethiopians and Egyptians. Am. J. Hum. Genet. 96, 986–991 (2015).

14. Malaspinas, A.-S. et al. A genomic history of aboriginal Australia. Nature 538, 207–214 (2016).

15. Mallick, S. et al. The Simons Genome Diversity Project: 300 genomes from 142 diverse populations. Nature 538, 201–206 (2016).

16. Posth, C. et al. Pleistocene mitochondrial genomes suggest a single major dispersal of non-Africans and a late glacial population turnover in Europe. Curr. Biol. 26, 827–833 (2016).

17. Bergström, A., Stringer, C., Hajdinjak, M., Scerri, E. M. L. & Sko- glund, P. Origins of modern human ancestry. Nature 590, 229–237 (2021).

18. Vallini, L. et al. Genetics and material culture support repeated expansions into paleolithic eurasia from a population hub out of Africa. Genome Biol. Evol. 14, evac045 (2022).

19. Slimak, L. The three waves: rethinking the Structure of the first Upper Paleolithic in Western Eurasia. bioRxiv https://doi.org/10. 1101/2022.10.28.514208. (2022).

20. Slimak, L. et al. Modern human incursion into Neanderthal terri- tories 54,000 years ago at Mandrin, France. Sci. Adv. 8, eabj9496 (2022).

21. Prüfer, K. et al. A genome sequence from a modern human skull over 45,000 years old from Zlatý kůň in Czechia. Nat. Ecol. Evol. 5, 820–825 (2021).

22. Marciani, G. et al. Lithic techno-complexes in Italy from 50 to 39 thousand years BP: an overview of lithic technological changes across the middle-upper Palaeolithic boundary. Quat. Int. 551, 123–149 (2020).

23. Benazzi, S. et al. Early dispersal of modern humans in Europe and implications for Neanderthal behaviour. Nature 479,
525–528 (2011).

24. Vallini, L. & Pagani, L. The future of the Eurasian past: highlighting plotholes and pillars of human population movements in the Late Pleistocene. J. Anthropol. Sci. 100, 231–241 (2022).

25. Hajdinjak, M. et al. Initial Upper Palaeolithic humans in Europe had recent Neanderthal ancestry. Nature 592, 253–257 (2021).

26. Fu, Q. et al. Genome sequence of a 45,000-year-old modern human from western Siberia. Nature 514, 445–449 (2014).

27. Yang, M. A. et al. 40,000-year-old individual from asia provides insight into early population structure in Eurasia. Curr. Biol. 27, 3202–3208.e9 (2017).

28. Posth, C. et al. Palaeogenomics of upper palaeolithic to neolithic European hunter-gatherers. Nature 615, 117–126 (2023).

29. Sikora, M. et al. The population history of northeastern Siberia since the Pleistocene. Nature 570, 182–188 (2019).

30. Lazaridis, I. et al. Genomic insights into the origin of farming in the ancient Near East. Nature 536, 419–424 (2016).

31. Lazaridis, I. et al. Paleolithic DNA from the Caucasus reveals core of West Eurasian ancestry. bioRxiv 423079 https://doi.org/10.1101/ 423079. (2018).

Nature Communications | (2024)15:1882

10

Article

https://doi.org/10.1038/s41467-024-46161-7

32. Allentoft, M. E. et al. Population genomics of stone age Eurasia. bioRxiv (2022). https://doi.org/10.1101/2022.05.04.490594. (2022).

33. Green, R. E. et al. A draft sequence of the Neandertal genome. Science 328, 710–722 (2010).

34. Seguin-Orlando, A. et al. Paleogenomics. Genomic structure in Europeans dating back at least 36,200 years. Science 346, 1113–1118 (2014).

35. Reich, D. et al. Genetic history of an archaic hominin group from Denisova Cave in Siberia. Nature 468, 1053–1060 (2010).

36. Sankararaman, S. et al. The genomic landscape of Neanderthal ancestry in present-day humans. Nature 507, 354–357 (2014).

37. Harris, K. & Nielsen, R. The genetic cost of Neanderthal intro- gression. Genetics 203, 881–891 (2016).

38. Moorjani, P. et al. The history of African gene flow into Southern Europeans, Levantines, and Jews. PLoS Genet. 7, e1001373 (2011).

39. Wall, J. D. Inferring human demographic histories of Non-African

    populations from patterns of allele sharing. Am. J. Hum. Genet.

    100, 766–772 (2017).

40. Wedage, O. et al. Microliths in the South Asian rainforest ~45-4 ka:

    mew insights from Fa-Hien Lena Cave, Sri Lanka. PLoS One 14,

    e0222606 (2019).

41. Reich, D., Thangaraj, K., Patterson, N., Price, A. L. & Singh, L.

    Reconstructing Indian population history. Nature 461,

    489–494 (2009).

42. Yelmen, B. et al. Ancestry-specific analyses reveal differential

    demographic histories and opposite selective pressures in mod-

    ern south Asian populations. Mol. Biol. Evol. 36, 1628–1642 (2019).

43. Raghavan, M. et al. Upper Palaeolithic Siberian genome reveals

    dual ancestry of Native Americans. Nature 505, 87–91 (2014).

44. Willerslev, E. & Meltzer, D. J. Peopling of the Americas as inferred

    from ancient genomics. Nature 594, 356–364 (2021).

45. Kelleher, J., Etheridge, A. M. & McVean, G. Efficient coalescent

    simulation and genealogical analysis for large sample sizes. PLoS

    Comput. Biol. 12, e1004842 (2016).

46. Gravel, S., Henn, B. M. & Gutenkunst, R. N. Demographic history

    and rare allele sharing among human populations. Proc. Natl.

    Acad. Sci. USA 108, 11983–8 (2011).

47. Kamm, J., Terhorst, J., Durbin, R. & Song, Y. S. Efficiently inferring

    the demographic history of many populations with allele count

    data. J. Am. Stat. Assoc. 115, 1472–1487 (2020).

48. Broushaki, F. et al. Early neolithic genomes from the eastern fertile

    crescent. Science 353, 499–503 (2016).

49. Marchi, N. et al. The genomic origins of the world’s first farmers.

    Cell 185, 1842–1859.e18 (2022).

50. Narasimhan, V. M. et al. The formation of human populations in

    South and Central Asia. Science 365, eaat7487 (2019).

51. Bramanti, B. et al. Genetic discontinuity between local hunter-

    gatherers and central Europe’s first farmers. Science 326,

    137–140 (2009).

52. Skoglund, P. et al. Origins and genetic legacy of neolithic farmers

    and hunter-gatherers in Europe. Science 336, 466–469 (2012).

53. Lazaridis, I. et al. Ancient human genomes suggest three ancestral

    populations for present-day Europeans. Nature 513,

    409–413 (2014).

54. Bortolini, E. et al. Inferring patterns of folktale diffusion using

    genomic data. Proc. Natl. Acad. Sci. USA 114, 9140–9145 (2017).

55. Shinde, V. et al. An ancient harappan genome lacks ancestry from

    steppe pastoralists or Iranian farmers. Cell 179, 729–735.e10

    (2019).

56. Pagani, L. & Crevecoeur, I. What is Africa? A human perspective.

    Human origins and dispersal, words, bones ... https://www. academia.edu/39001199/What_is_Africa_A_Human_ Perspective. (2019).

57. Beyer, R. M., Krapp, M. & Manica, A. High-resolution terrestrial climate, bioclimate and vegetation for the last 120,000 years. Sci. Data 7, 236 (2020).

58. Leonardi, M., Hallett, E. Y., Beyer, R., Krapp, M. & Manica, A. pastclim 1.2: an R package to easily access and use paleoclimatic reconstructions. Ecography https://doi.org/10.1111/ecog.
06481. (2023).

59. Rodríguez, J., Willmes, C., Sommer, C. & Mateos, A. Sustainable human population density in Western Europe between 560.000 and 360.000 years ago. Sci. Rep. 12, 6907 (2022).

60. Shoaee, M. J. et al. Defining paleoclimatic routes and opportu- nities for hominin dispersals across Iran. PLoS One 18, e0281872 (2023).

61. Zwyns, N. et al. The Northern Route for Human dispersal in Central and Northeast Asia: new evidence from the site of Tolbor-16, Mongolia. Sci. Rep. 9, 11759 (2019).

62. Djamali, M. et al. A late Pleistocene long pollen record from Lake Urmia, Nw Iran. Quat. Res. 69, 413–420 (2008).

63. Carolin, S. et al. Iranian speleothems: Investigating Quaternary climate variability in semi-arid Western Asia. in EPSC 2016–2923 (2016).

64. Mehterian, S. et al. Speleothem records of glacial/interglacial climate from Iran forewarn of future water availability in the interior of the Middle East. Quat. Sci. Rev. 164, 187–198 (2017).

65. Ghafarpour, A., Khormali, F., Meng, X., Tazikeh, H. & Stevens, T. Late Pleistocene climate and dust source from the Mobarakabad loess–paleosol sequence, northern foothills of the Alborz moun- tains, northern Iran. Front. Earth Sci. 9, https://www.frontiersin. org/articles/10.3389/feart.2021.795826/full. (2021).

66. Kehl, M. et al. Pleistocene dynamics of dust accumulation and soil formation in the southern Caspian Lowlands - new insights from the loess-paleosol sequence at Neka-Abelou, northern Iran. Quat. Sci. Rev. 253, 106774 (2021).

67. Torabi, M., Fattahi, M., Ghassemi, M. R., Amiri, A. & Saffar, M. M. Central Iranian playas: the first 3D near-surface model. J. Arid Environ. 208, 104847 (2023).

68. Shoaee, M. J., Vahdati Nasab, H. & Petraglia, M. D. The Paleolithic of the Iranian Plateau: Hominin occupation history and implica- tions for human dispersals across southern Asia. J. Anthropol. Archaeol. 62, 101292 (2021).

69. Ghasidian, E., Kafash, A., Kehl, M., Yousefi, M. & Heydari-Guran, S. Modelling Neanderthals’ dispersal routes from Caucasus towards east. PLoS One 18, e0281978 (2023).

70. Heydari, M., Guérin, G., Zeidi, M. & Conard, N. J. Bayesian lumi- nescence dating at Ghār-e Boof, Iran, provides a new chronology for Middle and Upper Paleolithic in the southern Zagros. J. Hum. Evol. 151, 102926 (2021).

71. Shea, J. J. Transitions or turnovers? Climatically-forced extinctions of Homo sapiens and Neanderthals in the east Mediterranean Levant. Quat. Sci. Rev. 27, 2253–2270 (2008).

72. Hershkovitz, I. et al. Levantine cranium from Manot Cave (Israel) foreshadows the first European modern humans. Nature 520, 216–219 (2015).

73. Heydari-Guran, S. & Ghasidian, E. Consistency of the ‘MIS 5 Humid Corridor Model’ for the Dispersal of Early Homo sapiens into the Iranian Plateau. In: Pearls, Politics and Pistachios: Essays in Anthropology and Memories on the Occasion of Susan Pollock’s 65th Birthday 219–238 (2021).

74. Nasab, H. V. et al. The open-air Paleolithic site of Mirak, northern edge of the Iranian Central Desert (Semnan, Iran): evidence of repeated human occupations during the late Pleistocene. C. R. Palevol. 18, 465–478 (2019).

75. Hashemi, S. M., Vahdati Nasab, H., Berillon, G. & Oryat, M. An investigation of the flake-based lithic tool morphology using 3D

Nature Communications | (2024)15:1882

11

Article

https://doi.org/10.1038/s41467-024-46161-7

geometric morphometrics: a case study from the Mirak Paleolithic

site. Iran. J. Archaeol. Sci. Rep. 37, 102948 (2021).

76. The middle and upper paleolithic archeology of the levant and

    beyond. (Springer Nature, Singapore).

77. Ghasidian, E. Rethinking the upper paleolithic of the Zagros

    mountains. PaleoAnthropology 2019, 240–310 (2019).

78. Ghasidian, E., Heydari-Guran, S. & Mirazón Lahr, M. Upper Paleo-

    lithic cultural diversity in the Iranian Zagros Mountains and the expansion of modern humans into Eurasia. J. Hum. Evol. 132, 101–118 (2019).

79. Pagani, L. et al. Genomic analyses inform on migration events during the peopling of Eurasia. Nature 538, 238–242 (2016).

80. Garba, R., Meredith-Williams, M., Neuhuber, S., Gier, S. & Usyk, V. Arid occupation of south-eastern Arabia: a new late pleistocene site at Wadi Asklat, south-central Oman—dating and paleoenvir- onmental reconstruction. https://meetingorganizer.copernicus. org/EGU23/EGU23-8109.html. (2023).

81. Heydari-Guran, S. et al. The discovery of an in situ Neanderthal remain in the Bawa Yawan Rockshelter, West-Central Zagros Mountains, Kermanshah. PLoS One 16, e0253708 (2021).

82. Iasi, L. N. M., Ringbauer, H. & Peter, B. M. An extended admixture pulse model reveals the limitations to human-neandertal intro- gression dating. Mol. Biol. Evol. 38, 5156–5174 (2021).

83. Zwyns, N. The initial upper paleolithic in central and east Asia: blade technology, cultural transmission, and implications for human dispersals. J. Paleolithic Archaeol. 4, 19 (2021).

84. Yang, S.-X. et al. Initial upper palaeolithic material culture by 45,000 years ago at Shiyu in northern China. Nat. Ecol. Evol. https://doi.org/10.1038/s41559-023-02294-4. (2024).

85. Tobler, R. et al. The role of genetic selection and climatic factors in the dispersal of anatomically modern humans out of Africa. Proc. Natl Acad. Sci. USA 120, e2213061120 (2023).

86. Aubert, M. et al. Pleistocene cave art from Sulawesi, Indonesia. Nature 514, 223–227 (2014).

87. Valladas, H. et al. Bilan des datations carbone 14 effectuées sur des charbons de bois de la grotte Chauvet. Bull. Soc. Bot. France 102, 109–113 (2005).

88. Pike, A. W. G. et al. U-series dating of Paleolithic art in 11 caves in Spain. Science 336, 1409–1413 (2012).

89. Teyssandier, N., Bon, F. & Bordes, J.-G. Within Projectile range: some thoughts on the appearance of the Aurignacian in Europe. J. Anthropol. Res. 66, 209–229 (2010).

90. Sano, K. et al. The earliest evidence for mechanically delivered projectile weapons in Europe. Nat. Ecol. Evol. 3, 1409–1414 (2019).

91. Metz, L., Lewis, J. E. & Slimak, L. Bow-and-arrow, technology of the first modern humans in Europe 54,000 years ago at Mandrin, France. Sci. Adv. 9, eadd4675 (2023).

92. Yaroshevich, A., Kaufman, D. & Marks, A. Weapons in transition: reappraisal of the origin of complex projectiles in the Levant based on the Boker Tachtit stratigraphic sequence. J. Archaeol. Sci. 131, 105381 (2021).

93. Langley, M. C. et al. Bows and arrows and complex symbolic displays 48,000 years ago in the South Asian tropics. Sci. Adv. 6, eaba3831 (2020).

94. Lazaridis, I. et al. The genetic history of the Southern Arc: a bridge between West Asia and Europe. Science 377, eabm4247 (2022).

95. Teschler-Nicola, M. et al. Ancient DNA reveals monozygotic newborn twins from the Upper Palaeolithic. Commun. Biol. 3, 1–11 (2020).

96. Fu, Q. et al. The genetic history of ice age Europe. Nature 534, 200–205 (2016).

97. Mao, X. et al. The deep population history of northern East Asia from the Late pleistocene to the holocene. Cell 184, 3256–3266.e13 (2021).

98. Prüfer, K. et al. A high-coverage Neandertal genome from Vindija Cave in Croatia. Science 358, 655–658 (2017).

99. Gallego Llorente, M. et al. Ancient Ethiopian genome reveals extensive Eurasian admixture throughout the African continent. Science 350, 820–822 (2015).

100. Jones, E. R. et al. Upper Palaeolithic genomes reveal deep roots of modern Eurasians. Nat. Commun. 6, 8912 (2015).

101. Harney, É. et al. Ancient DNA from Chalcolithic Israel reveals the role of population mixture in cultural transformation. Nat. Com- mun. 9, 3336 (2018).

102. Agranat-Tamir, L. et al. The genomic history of the bronze age Southern Levant. Cell 181, 1146–1157.e11 (2020).

103. Antonio, M. L. et al. Ancient rome: a genetic crossroads of Europe and the Mediterranean. Science 366, 708–714 (2019).

104. Feldman, M. et al. Late Pleistocene human genome suggests a local origin for the first farmers of central Anatolia. Nat. Commun. 10, 1218 (2019).

105. Haber, M. et al. Continuity and admixture in the last five millennia of levantine history from ancient canaanite and present-day Lebanese genome sequences. Am. J. Hum. Genet. 101,
274–282 (2017).

106. van de Loosdrecht, M. et al. Pleistocene North African genomes link near eastern and sub-Saharan African human populations. Science 360, 548–552 (2018).

107. Mathieson, I. et al. The genomic history of southeastern Europe. Nature 555, 197–203 (2018).

108. González-Fortes, G. et al. Paleogenomic evidence for multi- generational mixing between neolithic farmers and mesolithic Hunter-gatherers in the lower Danube basin. Curr. Biol. 27, 1801–1810.e10 (2017).

109. Svensson, E. et al. Genome of Peştera Muierii skull shows high diversity and low mutational load in pre-glacial Europe. Curr. Biol. 31, 2973–2983.e9 (2021).

110. Prüfer, K. et al. The complete genome sequence of a Neanderthal from the Altai Mountains. Nature 505, 43–49 (2014).

111. Mafessoni, F. et al. A high-coverage Neandertal genome from Chagyrskaya Cave. Proc. Natl Acad. Sci. USA 117,
15132–15136 (2020).

112. Meyer, M. et al. A high-coverage genome sequence from an archaic Denisovan individual. Science 338, 222–226 (2012).

113. Slon, V. et al. The genome of the offspring of a Neanderthal mother and a Denisovan father. Nature 561, 113–116 (2018).

114. Allentoft, M. E. et al. Population genomics of Bronze Age Eurasia. Nature 522, 167–172 (2015).

115. Mathieson, I. et al. Genome-wide patterns of selection in 230 ancient Eurasians. Nature 528, 499–503 (2015).

116. Sikora, M. et al. Ancient genomes show social and reproductive behavior of early upper paleolithic foragers. Science 358, 659–662 (2017).

117. Skourtanioti, E. et al. Genomic history of neolithic to bronze Age Anatolia, Northern Levant, and Southern Caucasus. Cell 181, 1158–1175.e28 (2020).

118. Kılınç, G. M. et al. The demographic development of the first farmers in Anatolia. Curr. Biol. 26, 2659–2666 (2016).

119. Gokhman, D. et al. Differential DNA methylation of vocal and facial anatomy genes in modern humans. Nat. Commun. 11, 1189 (2020).

120. Hofmanová, Z. et al. Early farmers from across Europe directly descended from Neolithic Aegeans. Proc. Natl Acad. Sci. USA 113, 6886–6891 (2016).

121. de Barros Damgaard, P. et al. The first horse herders and the impact of early Bronze Age steppe expansions into Asia. Science 360, eaar7711 (2018).

122. Bergström, A. et al. Insights into human genetic variation and population history from 929 diverse genomes. Science 367, eaay5012 (2020).

Nature Communications | (2024)15:1882

12

Article

https://doi.org/10.1038/s41467-024-46161-7

123. 1000 Genomes Project Consortium. et al. A global reference for human genetic variation. Nature 526, 68–74 (2015).

124. Mondal, M. et al. Genomic analysis of Andamanese provides insights into ancient human migration into Asia and adaptation. Nat. Genet. 48, 1066–1070 (2016).

125. Raghavan, M. et al. POPULATION GENETICS. Genomic evidence for the Pleistocene and recent population history of Native Americans. Science 349, aab3884 (2015).

126. Skoglund, P. et al. Genetic evidence for two founding populations of the Americas. Nature 525, 104–108 (2015).

127. Fan, S. et al. African evolutionary history inferred from whole genome sequence data of 44 indigenous African populations. Genome Biol. 20, 82 (2019).

128. Chang, C. C. et al. Second-generation PLINK: rising to the chal- lenge of larger and richer datasets. Gigascience 4, s13742–015–0047–8 (2015).

129. Patterson, N. et al. Ancient admixture in human history. Genetics 192, 1065–1093 (2012).

130. Alexander, D. H., Novembre, J. & Lange, K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 19, 1655–1664 (2009).

131. Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modeling of species geographic distributions. Ecol. Modell. 190, 231–259 (2006).

132. Timmermann, A. et al. Climate effects on archaic human habitats and species successions. Nature 604, 495–501 (2022).

133. Raia, P. et al. Past extinctions of homo species coincided with increased vulnerability to climatic change. One Earth 3, 480–490 (2020).

134. Leonardi, M., Hallett, E. Y., Beyer, R., Krapp, M. & Manica, A. pastclim: an R package to easily access and use paleoclimatic reconstructions. bioRxiv (2022). https://doi.org/10.1101/2022.05. 18.492456.

135. Harrell Jr, F., & Dupont, Ch. Hmisc Harrell miscellaneous. R package version 4.2-0. - references - scientific research publish- ing. https://www.scirp.org/(S(vtj3fa45qm1ean45%20vvffcz55))/ reference/referencespapers.aspx?referenceid=3118318. (2019).

136. Barbosa, A. M. fuzzySim: applying fuzzy logic to binary similarity indices in ecology. Methods Ecol. Evol. 6, 853–858 (2015).

137. Duong, T. Kernel Smoothing [R package ks version 1.14.0]. (2022).

138. Barber, R. A., Ball, S. G., Morris, R. K. A. & Gilbert, F. Target‐group backgrounds prove effective at correcting sampling bias in Max- ent models. Divers. Distrib. 28, 128–141 (2022).

139. NOAA National Geophysical Data Center. ETOPO1 1 arc-minute global relief model. (2009).

140. Amante, C. & Eakins, B. W. ETOPO1 global relief model converted to PanMap layer format. (2009) https://doi.org/10.1594/ PANGAEA.769615. (2009).

141. Bintanja, R., van de Wal, R. S. W. & Oerlemans, J. Modelled atmospheric temperatures and global sea levels over the past million years. Nature 437, 125–128 (2005).

142. Vallini, L. et al. The Persian Plateau served as Hub for Homo sapiens after the main out of Africa dispersal, leovallini/ValliniHub2024: ValliniHub2024. https://doi.org/10.5281/zenodo.10571649.

Acknowledgements

The authors would like to thank Dr. Burak Yelmen for sharing deconvo- luted haplotypes of Irula_S genomes. Computation carried out thanks to University of Padova Strategic Research Infrastructure Grant 2017: “CAPRI: Calcolo ad Alte Prestazioni per la Ricerca e l’Innovazione”. EB, GM, and SB are funded by ERC-724046-SUCCESS “The earliest migra- tion of Homo sapiens in Southern Europe. Understanding the biocultural processes that define our uniqueness” awarded to Stefano Benazzi. LP is funded by the Italian Ministry of University and Research (PRIN 2022B27XYM).

Author contributions

Conceived the study: LV, LP. Genomic data analysis: LV, MM, LP. Palaeoecological modelling: CZ, AB, LV. Provided archaeological con- text: MJS, EB, GM, MDP. Wrote the manuscript: LV, CZ, MJS, MM, MDP, LP. revised the manuscript and provided interpretation of results: EB, GM, SA, TP, SB.

Competing interests

The authors declare no competing interests

Additional information

Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41467-024-46161-7.

Correspondence and requests for materials should be addressed to Leonardo Vallini or Luca Pagani.

Peer review information Nature Communications thanks Pasquale Raia, Yassine Souilmi and the other, anonymous, reviewer(s) for their con- tribution to the peer review of this work. A peer review file is available.

Reprints and permissions information is available at http://www.nature.com/reprints

Publisher’s note Springer Nature remains neutral with regard to jur- isdictional 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 licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence 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 licence, visit http://creativecommons.org/ licenses/by/4.0/.

© The Author(s) 2024

Nature Communications | (2024)15:1882

13