Projecting Ancient Ancestry in Modern-Day Arabians and Iranians: A Key Role of the Past Exposed Arabo-Persian Gulf on Human Migrations

Abstract

The Arabian Peninsula is strategic for investigations centered on the early structuring of modern humans in the wake of the out-of-Africa migration. Despite its poor climatic conditions for the recovery of ancient human DNA evidence, the availability of both genomic data from neighboring ancient specimens and informative statistical tools allow modeling the ancestry of local modern populations. We applied this approach to a data set of 741,000 variants screened in 291 Arabians and 78 Iranians, and obtained insightful evidence. The west-east axis was a strong forcer of population structure in the Peninsula, and, more importantly, there were clear continuums throughout time linking western Arabia with the Levant, and eastern Arabia with Iran and the Caucasus. Eastern Arabians also displayed the highest levels of the basal Eurasian lineage of all tested modern-day populations, a signal that was maintained even after correcting for a possible bias due to a recent sub-Saharan African input in their genomes. Not surprisingly, eastern Arabians were also the ones with highest similarity with Iberomaurusians, who were, so far, the best proxy for the basal Eurasians amongst the known ancient specimens. The basal Eurasian lineage is the signature of ancient non-Africans who diverged from the common European-eastern Asian pool before 50,000 years ago, prior to the later interbred with Neanderthals. Our results appear to indicate that the exposed basin of the Arabo-Persian Gulf was the possible home of basal Eurasians, a scenario to be further investigated by searching ancient Arabian human specimens.

Arabian Peninsula, Iran, basal Eurasian lineage, ancient and archaic ancestry, out of Africa migration, main human population groups stratification

Introduction

Significance

Our work is an excellent example that the availability of informative ancient DNA samples from key geographic regions allows to directly test alternative scenarios for past founder contributions to the genomic pools of extant populations. We applied this approach to Arabian Peninsula (AP) and Iran populations. We detected a strong degree of isolation between western and eastern AP populations, with two clear continuums throughout time: one linking western Arabia and the Levant; and the other showing that eastern Arabia was a migration path with Iran and Caucasus. The eastern continuum was probably older, given the stronger signal of the basal Eurasian lineage and Iberomaurusians in that side of the AP. The basal Eurasian lineage is the signature of ancient non-Africans who diverged from the common European-eastern Asian pool before 50,000 years ago, prior to the interbreeding with Neanderthals. Our results indicate that the exposed basin of the Arabo-Persian Gulf must be considered amongst the possible homes of basal Eurasians, a scenario to be further investigated by continuing the pursuit for ancient Arabian human specimens.

Coalescence analyses on modern-day genomic data have been fundamental for inferring major events of the human past (Li et al. 2008; Soares et al. 2012; Hellenthal et al. 2014). By 2010, the successful technical improvements on ancient DNA (aDNA) cataloging contributed essential insights into the overall picture (Skoglund and Mathieson 2018). In fact, the new possibility of surveying diversity that became extinct is revealing an extremely rich picture involving the interplay between several population groups. New insights have been obtained on the out-of-Africa (OOA) migration (Stringer and Andrews 1988; Stringer 2014; Llorente et al. 2015; Schuenemann et al. 2017), namely, the divergence between the non-African populations (Lipson and Reich 2017; Lazaridis 2018; Yang and Fu 2018), and their earliest movements into Eurasia and America (Fu et al. 2014, 2016; Lazaridis et al. 2014; Raghavan et al. 2014; Rasmussen et al. 2014; Posth et al. 2018). The aDNA catalog already contains hundreds of whole-genome shotgun data and thousands of targeted capture data (single nucleotide polymorphism [SNP] array-like) from specimens that date as far back as 45 thousand years ago (ka). Unfortunately, there is still a strong bias towards European, Siberian, Native American, and Near Eastern specimens, so large gaps remain in the reconstruction of Asian and African genomic histories (Skoglund and Mathieson 2018; Pereira et al. 2021). A difficult-to-overcome factor responsible for this bias is the environmental influence in aDNA preservation, with arid deserts (AP, Sahara, and others) and humid tropical forests (large parts of Africa, Asia, and America) amongst the worst possible conditions (Hambrecht and Rockman 2017). Nevertheless, even for these regions, newly developed statistical tools are allowing to evaluate the contribution of neighboring ancient founders to the genomes of current populations (Leslie et al. 2015; Montinaro et al. 2015; Raveane et al. 2019).

AP played a major role or, at least, was in the path of the successful OOA migration at around 70–60 ka, and important genetic and archeology findings have been uncovered in the region (Rose et al. 2011; Fernandes et al. 2012; Petraglia et al. 2012; Groucutt et al. 2018; Scerri et al. 2018). Pleistocene Arabia was exposed to several climate change episodes that impacted modern human occupation (Petraglia et al. 2020). During wet periods, populations would expand from the coast to internal regions following river valleys, while contracting to refugia during arid periods. Climatic and archeological evidences indicate the existence of three main AP refugia (Rose 2010): the Red Sea coastal plain; the Dhofar Mountains and adjacent littoral zone in Yemen and Oman; and the exposed basin of the Arabo-Persian Gulf. Genomic studies focusing on extant AP populations (Alshamali et al. 2009; Cerný et al. 2009, 2011; Fernandes et al. 2012, 2015, 2019; Musilová et al. 2011; Rodriguez-Flores et al. 2016; Vyas et al. 2017) revealed a significant heterogeneity between western and eastern AP populations. This heterogeneity consists in a gradual cline most probably due to western AP (Saudi Arabia and Yemen) forming a somewhat continuum with eastern sub-Saharan Africa and the Levant, whereas eastern AP (Oman and Unite Arab Emirates or UAE) received inputs from Iranian and South Asian populations. Recently, we demonstrated (Fernandes et al. 2019) that the axis also matches distinctive positive selection signals for variants conferring malaria protection (higher in the west) and lactose tolerance (local adaptation in the western AP; adaptive admixture of the European/South Asian-derived allele in the eastern AP). However, analyses of basal mitochondrial DNA haplogroup N (Fernandes et al. 2012) suggested that the successful OOA migration had a more significant impact in the then exposed basin of the Arabo-Persian Gulf.

Despite human ancient remains being extremely scarce in AP, with archeological sites containing mainly tools and artifacts (Groucutt et al. 2018; Petraglia et al. 2012, 2020; Rose 2010; Rose et al. 2011), aDNA information from neighboring regions and from important periods can provide important clues to reconstruct the genetic history of its past inhabitants. A main issue concerns the ancestral group of successful OOA migrants, who have been previously designated as the bearers of the basal Eurasian lineage (Fu et al. 2015). The basal Eurasian bearers diverged from the common European-eastern Asian pool prior to 50 ka (Lazaridis et al. 2014; Yang and Fu 2018), hence they are necessarily an outgroup in relation to all known deep ancient hunter–gatherer (HG) Eurasian specimens, such as the 45 ka Ust’-Ishim (Fu et al. 2014) and the 39 ka Oase 1. The geographical location of this population group is still a matter of discussion, with higher proportions of the basal Eurasian lineage found in ancient samples from the Caucasus Mountains (Satsurblia and Kotias specimens; 13–9 ka; Jones et al. 2015), the Levant and Iran (14–3.5 ka; Lazaridis et al. 2016). To a lesser extent, the lineage is still present in modern populations from western Eurasia and southwest Asia (Lazaridis et al. 2016). It has been previously speculated that isolated basal Eurasian lineage descendants in the central Zagros Mountains of Iran, who were the first goatherds, spread afterwards into the Eurasian steppe (Broushaki et al. 2016; Lazaridis et al. 2016). Also the ancient Iberomaurusian specimens could be in part descendants of basal Eurasians (van de Loosdrecht et al. 2018), as they have a shared genetic affinity with early Holocene Near Easterners (best proxy are the Natufians) and one-third of input from sub-Saharan Africans (a mixture between western and eastern Africans). Thus, so far, southwestern Asia is still the most probable homeland for the hypothetical basal Eurasian population, although a more precise location is still difficult to pinpoint. Another element not to neglect is the little, if any, Neanderthal input in the basal Eurasian lineage (Lazaridis et al. 2016). This supports a Neanderthal-modern human admixture event after the basal Eurasian lineage divergence, but prior to the separation of the European and Asian pools. The later-on admixture of the basal Eurasian descendants with western Eurasians and southwest Asians possibly explains their observed lower Neanderthal input (Rodriguez-Flores et al. 2016; Vyas and Mulligan 2019; Almarri et al. 2021) than in present-day East Asians (no basal Eurasian input). Again, this new information for the Neanderthal-modern human interbreeding event can help to refine the geographic source of the basal Eurasian ancestors. The presence of modern humans in Israeli archeological sites, such as at Skhul and Qafzeh (Grün and Stringer 1991; Hershkovitz et al. 2015), and Neanderthals at Tabun (Simpson et al. 1998) and Kebara Cave (Trinkaus et al. 1991; Valladas et al. 1987) dating between 50 and 120 ka provides evidence that modern humans and Neanderthals overlapped geographically and temporally in the Levant (Hershkovitz et al. 2015). No such evidence was found in the southern parts of southwest Asia, although absence of evidence is not evidence of absence.

In this work, we merged and analyzed ancient genomes from key neighbors with the genomic profiles of modern-day Arabian populations in order to test alternative scenarios explaining their origins. Our data set consisted in 741,000 variants screened in 291 Arabians and 78 Iranians that we previously analyzed in detail in terms of population structure and selection events (Fernandes et al. 2019). Our work differs from others due to the good geographical resolution across the entire AP.

Results

Projection of Ancient Ancestry on Current Day Arabian and Iranian Diversities

To understand how the Arabian and neighboring populations are related to groups that occupied the region since the transition from the Paleolithic till the Bronze Age, we performed a PCA (principal components analysis) projecting ancient samples onto PCs estimated on modern allele-frequency variation (fig. 1). Modern Saudi Arabian and Yemeni samples clustered tightly, overlapping with the three Natufian samples, and were close to the Levant Pre-Pottery Neolithic B and C (PPNB and PPNC) and Levant Bronze Age samples. Most modern UAE, Oman and the Iranian samples mingled with Iran-Late Neolithic and -Neolithic, and Caucasus-Paleolithic, -Early Bronze Age, and -Copper Age. A few modern Arabian and Iranian samples were more dispersed and mingled with North African-Iberomaurusian and North African-Early Neolithic, in accordance with some sub-Saharan input in their genomes (supplementary fig. S1, Supplementary Material online displays a PCA with sub-Saharan and East Asian samples included; supplementary fig. S2, Supplementary Material online displays only the modern samples in colors for an easier visualization). This dispersion also shows the existence of some structure within each population. The ancient Anatolian, Balkan, Steppe, and European samples were projected distantly from the Arabian and Iranian modern samples.

Fig. 1.

PCA of present-day samples from the AP (western AP in blue and eastern AP in green), Iran (red) and neighboring populations (gray, to avoid visual clutter) with projected ancient samples in different colors and formats according to the legend. Abbreviations: E, Early; M, Middle; L, Late; HG, hunter–gatherer; N, Neolithic; C, Chalcolithic; BA, Bronze Age; IA, Iron Age; UP, upper Paleolithic.

The population structure pattern just described for the PCA was largely corroborated by the ADMIXTURE analysis (fig. 2A and B), for the optimal K (K = 9) estimated through cross-validation (CV) of the logistic regression. The dark blue component was frequent in modern AP populations (highest values in Saudi Arabia and Yemen) and the Near East, and in the ancient Levant Natufians and Levant Bronze Age samples. In contrast, a substantial proportion of modern samples from east AP and Iran shared a high frequency of the light purple component with Zagros-Neolithic and Anatolia/Caucasus across all periods; this component has been found modal in Levant and Caucasus populations. The orange component, which was very frequent in modern northern European samples, was also dominant in the HG samples from Europe and the Balkans and in the steppe populations. The light green component, the second most common in modern Europeans, especially from the south, and also present in modern Near easterns, Anatolians, and Iranians, was observed in all ancient samples (especially in Anatolian- and Southern European-Neolithic) except in the Zagros-Neolithic ancient sample. In the Zagros-Neolithic ancient sample, as well as in the ancient steppe populations, a considerable amount of the dark green component was detected. This component, typical of modern South Asian populations, was also observed in modern samples from Iran, UAE, and Oman. The impact of the African ancestry, represented by the pink and red components, was especially high in the ancient South Eastern African samples and has been found with important frequencies in ancient North African individuals (especially in Iberomaurusian samples). The same components were also found in part of modern samples from UAE and Oman, together with a few individuals from Yemen.

Fig. 2.

ADMIXTURE model-based clustering analysis of present-day populations (A) and projected ancient samples (B) for the optimal K (K = 9) estimated through CV of the logistic regression. Each individual is represented by a vertical (100%) stacked column of genetic components proportions shown in color. India* corresponds to India Gujarati, India** to India Andhra Pradesh, China* to China Xishuangbanna and China** to China Beijing. UP* corresponds to the Iberomaurusian samples, EPI* to the Natufian samples, N* to the Zagros samples, HG* to samples from Russia and IA* to the steppe samples. Abbreviations: E, Early; M, Middle; L, Late; HG, hunter–gatherer; N, Neolithic; C, Chalcolithic; BA, Bronze Age; IA, Iron Age; UP, upper Paleolithic; EPI, Epipaleolithic.

In order to better ascertain the impact of the sub-Saharan African input in the AP genomes (individuals with values > 40% sub-Saharan input had been previously removed from the analysis, to avoid bias by very recent migrants), we checked the population structure through a haplotype-based method (fineSTRUCTURE; fig. 3A and B and supplementary fig. S3, Supplementary Material online). A few clusters supported the differentiation between the west-east AP axis: cluster F is typical of west AP, representing 2/3 of Saudi and Yemeni populations, with less than 1/3 frequency in UAE and Oman; cluster D is characteristic of east AP, also shared with Iran, and practically inexistent in west AP; cluster A, was highly frequent in east AP (nearing 50%), and more restricted in west AP (around 20%). Cluster C, typical from the Levant showed quite limited presence in Saudi, UAE, Oman, and Iran; this cluster is related with cluster B only observed in Druze people (an isolated group in Levant). These ancestry backgrounds for each cluster are easily visualized when the ADMIXTURE plot is reorganized by cluster (fig. 3A). The mean sub-Saharan input values in the modern AP and Iranian samples distributed by these clusters were the following: 22% in A; 5% in C; 6% in D; 0.3% in E; and 3% in F. Definitely, cluster A still has a considerable input from sub-Saharan Africa, but all the others have limited values around 5%. We further tried to ascertain the lengths of the sub-Saharan blocks in clusters A and D, using a local ancestry algorithm (RFMix; supplementary fig. S4, Supplementary Material online; Maples et al. 2013). The RFMix software employs a linkage disequilibrium (LD) model between markers to infer ancestry for each segment of the genome between a mixture of putative ancestral panels of haplotypes. Supplementary fig. S4 and S5, Supplementary Material online display the median sizes (in centimorgan [cM]) of the sub-Saharan African blocks for each Arabian and Iranian individuals affiliated in clusters A and D. As recombination occurs along time, it will decrease the size of the blocks, so that the time of admixture is inversely correlated with the size of the blocks. The results indicated a clear tendency for shorter blocks in both clusters, compatible with an older age of sub-Saharan African admixture.

Fig. 3.

Population substructure in AP and neighboring regions (see supplementary fig. S3, Supplementary Material online for detailed tree). (A) FineSTRUCTURE tree with population clusters highlighted and respective ADMIXTURE profiles (the colors in this scheme are identical to fig. 2). The outgroup G comprises European populations. Cluster B is made up of only Druze individuals, who were not represented in (B). (B) Distribution of the clusters per AP country, in Iran, the Levant (summing up Palestinians, Jordanians, Syrians, and Lebanese individuals), Caucasus (Armenians, Georgians, and Lezgins), and Anatolia (Turkey).

Modeling the Ancestry Composition of AP

First, we modeled the ancestry composition of AP per country from ancient source populations (qpAdm algorithm). The most complex scenario that was feasible for the AP and Iranian populations tested here consisted in four sources (fig. 4A and supplementary table S1, Supplementary Material online): Caucasus HG; Iberomaurusian; Natufian; and Zagros farmer. Two-thirds of the west AP populations have ancestry shared with the Natufian, typical of the Levant, whereas east AP and Iran had a predominance of background shared with Caucasus HG and Zagros farmers. The similarity with the Iberomaurusian was higher in eastern than in western AP, but already residual in Iran. All the other feasible four-source scenarios (supplementary table S2, Supplementary Material online) consisted in a predominant Levant signal in the west (replacing either Caucasus or Zagros by Levant PPNB or Levant early bronze age), whereas the east had always a main Caucasus (sometimes Caucasus early bronze age replacing totally the Levant signal) and Zagros farmer similarity.

Fig. 4.

The ancestry composition of AP. (A) Mixture proportions on AP and Iran populations inferred by qpAdm as a combination of four ancient sources (this combination was the only one feasible in all the five tested populations). (B) CP/NNLS results for the same set of four ancient sources identified in the qpAdm analysis for each country. (C) CP/NNLS results for the same set of four ancient sources identified in the qpAdm analysis for each cluster observed in AP. (D) Basal Eurasian lineage contribution inferred through the statistic f4(Test, Han; Ust’-Ishim, Mbuti), where Ust’-Ishim and Mbuti were used as proxies for nonbasal and basal Eurasian lineages, respectively. The bullet identifies significantly negative values (Z-score < −3), whereas the cross identifies nonsignificant values, light red indicates modern populations and dark red marks ancient specimens. (E) Correlation between the basal Eurasian lineage estimates (inferred by qpAdm) and the Neanderthal introgression (measured through the f4(Test, Mbuti; Altai, Denisovan) test).

In order to have an independent estimation for admixture proportions of ancestries, we used another algorithm (CP/NNLS—CHROMOPAINTER analysis with the adaptation of the Nonnegative Least Squares [NNLS]) in the five modern populations (supplementary fig. S5, Supplementary Material online) divided by sources belonging to different archeological periods (Paleolithic, Neolithic, and Bronze Age). These results reflected the similarity within the western AP group with a continuum with Levant (Natufian, Levant PPNB, and Levant early Bronze Age), whereas the eastern AP-Iran group displayed more Caucasian/South Asian influence (Caucasus HG, Zagros Neolithic, and Steppe Bronze Age). As the four sources identified in qpAdm seem to also be highlighted in the CP/NNLS results per period, we performed an overall CP/NNLS analysis by testing Caucasus HG, Iberomaurusian, Natufian; and Zagros farmer as sources (fig. 4B and supplementary table S3, Supplementary Material online). As shown in fig. 4B, relative proportions of each of these sources are quite parallel between qpAdm and CP/NNLS inferences (linear correlation; r² = 0.8632), with a slight exception for Yemen, with a lower Zagros input and a higher Caucasus HG signal.

Given the high frequency in east AP of cluster A that has a significant sub-Saharan African input, we also calculated the CP/NNLS with the four sources (Caucasus HG; Iberomaurusian; Natufian; and Zagros farmer) per cluster (fig. 4C and supplementary table S4, Supplementary Material online). Only three clusters had input from Iberomaurusian: this input was overwhelming, and probably biased by a recent influence, in cluster A; it was very restricted in the east AP/Iranian-cluster D and more so in the Levant-cluster C. Curiously, the western AP-cluster F had significantly more influence from the Levant-emerged Epipaleolithic culture Natufian than the Levant–cluster C, Druze-cluster B, and Anatolian/Caucasus-cluster E. Zagros influence was restricted to the east AP/Iranian-cluster D, west AP-cluster F and cluster A.

Basal Eurasian Lineage

We followed Lazaridis et al. (2016) in applying the f4 (Test, Han; Ust’-Ishim, Mbuti) statistics to test for Basal Eurasian ancestry, by measuring the excess of allele sharing of Ust’-Ishim (ancient proxy for nonbasal ancestry) with a variety of Test populations compared with Han as a baseline. Figure 4D displays significant negative results in most of the tested modern populations from southwest Asia, compatible with the presence of the basal Eurasian lineage in these populations. The significant results were more extreme in eastern AP followed by western AP, the Near Eastern populations and Iran, Caucasus, Anatolia, and finally Europe. The ancient Iberomaurusian sample included in the analysis was the one with the highest input of the basal Eurasian lineage, even higher than the ones observed in modern eastern AP populations. Not surprisingly, cluster A still bearing a mean 20% sub-Saharan African admixture, was in a close position to the Iberomaurusian sample, with clusters D (eastern AP-Iran) and F (western AP) occupying a position similar to the other Levant and Iranian samples. The other ancient samples displayed similar f4 statistic values to local modern-day populations. The European Louschbour HG specimen, as expected, showed no input from the basal Eurasian lineage. To test if the signal is Basal Eurasian ancestry as opposed to unmodeled African ancestry, we retested f4 by changing the Mbuti proxy by: 1) Tswana from Botswana; 2) Biaka from central African Republic; 3) Khoisan from Namibia; 4) Yoruba from Nigeria; and 5) Khoisan from South Africa. The f4 statistic values were concordant in all these tests (supplementary fig. S6, Supplementary Material online).

In accordance with previous reports for ancient samples (Lazaridis et al. 2016; Vyas and Mulligan 2019), the proportion of the basal Eurasian lineage in the modern samples was negatively correlated with the Neanderthal input (r² = 0.92; P = 7.03 × 10⁻¹¹; fig. 4E; supplementary tables S5 and S6, Supplementary Material online). Concordantly with the above results, the modern samples with the highest amount of the basal Eurasian lineage were the east AP populations (∼45%), followed by the west AP and Iran (∼38%), Levant (∼32%; except Druze which displayed a lower 28%) and Caucasus (20–25%), whereas Europeans had the lowest values (<20%). In this analysis, the basal Eurasian lineage proportions between clusters A and D were quite similar (49% and 45%, respectively), indicating that the bias introduced by a probable recent sub-Saharan input in AP is very low. These proportions of the basal Eurasian lineage in the clusters more abundant in eastern AP were higher than the value observed for the dominant cluster in west AP (cluster F; 35%). The ancient specimens follow the regression line (the two outliers have higher missing genotypes), with the Iberomaurusian specimen displaying the highest basal Eurasian and the lowest Neanderthal proportions, as expected from the previous results. The qpAdm-based inference for the basal Eurasian input indicates comparable amounts in modern AP and ancient samples from Caucasus and Zagros, with a lower input of the Natufian.

As the basal Eurasian lineage input in AP could be due to admixture from later and already admixed sources, f4 tests were conducted in the form of f4(Test, Han; Ust’-Ishim, X), where X was an admixed ancient sample. Our results (supplementary fig. S7, Supplementary Material online) revealed a preponderant signature of the Natufian in west AP that was nevertheless low impacting in eastern AP (in fact, it had a higher influence in southern European populations than here). Interestingly, the test with the Iberomaurusian highlighted the similarity with the Levant and western AP, concordant with a main Near Eastern background of that specimen (van de Loosdrecht et al. 2018), followed by eastern AP (probably the basal Eurasian similarity), and then Iberia and Italy, an interesting result for the discussion of North Africa-South Europe connections (Hellenthal et al. 2014; Arauna et al. 2019). When the Basal signature is included in the evaluation, the Iberomaurusian specimen clusters with eastern AP, when left out it clusters with western AP. The Caucasus and Zagros specimens, themselves displaying several signals of admixture involving the periods under analysis (Jones et al. 2015), had a less clear pattern, so that founder and receptor become already too complex to disentangle.

Conclusions

In this work, we have shown that the availability of informative aDNA samples from neighboring key geographic regions allows to directly test alternative scenarios for past founder contributions to the genomic pools of extant populations (fig. 5).

Fig. 5.

Schematic representation of the hypothesis discussed in this article. (A) The OOA group of migrants, who followed a northern route (represented by the blue dashed arrow) and/or a southern route (represented by the red dashed arrow), diverged in two groups. One gave rise into the basal Eurasian group, hypothetically located in the Gulf region (peach highlighted region). The other interbred with Neanderthals, hypothetically in the Levant zone (green highlighted region), and eventually originating the Eurasian and Eastern Asian populations (originally without basal Eurasian ancestry). (B) Main signatures of population relationships in the period after 15 ka, for which informative aDNA of AP neighboring regions are beginning to be available (as indicated in the figure, with calibrated years before present [calBP]; geographical locations are approximate). There is a continuum between the Levant and western AP (purple highlight), in parallel to the continuum between Caucasus, Zagros, and eastern AP (yellow highlight). Gene flow between the two continua was probably stronger in the north (represented by the blue bidirectional arrow) than in the south, rendering basal Eurasian ancestry to be detected in Natufian and Iberomaurusian specimens.

It is now evident the continuum along time linking eastern AP with Iran and Caucasus, as reflected by the outcomes of the various analyses performed here. We extended the geographic scope (southwards) influenced by the two separated HG to farmer transitions, one in the southern Levant and the other in the southern Caucasus–Iran highlands, reported by others (Gallego-Llorente et al. 2016; Lazaridis et al. 2016). Our results reveal that the eastern continuum was probably older, given the stronger signal of the basal Eurasian lineage found in eastern AP. This designation of basal Eurasian focuses on the subsequent events that gave rise to the European and Asian population groups, but this lineage could be equally well described as the successful OOA migrant group highlighting the first successful settlement of a non-African region. This second designation calls the attention for a possible bias in identifying this signature due to a recent sub-Saharan input. By performing haplotype-based clustering and identifying clusters where the recent sub-Saharan input is negligible, we confirmed the basal Eurasian lineage signature in cluster D (eastern AP and Iran) followed by F (western AP).

The Iberomaurusian specimens are the ones bearing a very high level of basal Eurasian lineage because of two factors: 1) they derive from the main southwestern Asian parental population that migrated back to North Africa 20 ka; 2) they admixed with sub-Saharan Africans after their settlement in North Africa. Therefore, the Iberomaurusian specimens are probably the best ancient proxy for the basal Eurasian population group.

The highest levels of basal Eurasian lineage found in eastern AP reinforce previous genetic and archeological evidence of early human settlement in the exposed basin of the Arabo-Persian Gulf (Alshamali et al. 2009; Cerný et al. 2011; Al-Abri et al. 2012; Fernandes et al. 2012, 2015; Rodriguez-Flores et al. 2016; Vyas et al. 2017). These results, by relocating the geographic origin of the basal Eurasian population group from the Near East to the eastern part of the AP, reconcile better with the Neanderthal-modern human interbreeding event having probably occurred in the Levant. Lazaridis et al. (2016), who observed that the aDNA specimens with higher basal Eurasian lineage were from the Levant, suggested that either basal Eurasians thoroughly admixed into the Near East necessarily before the date of the analyzed Natufians but after the Neanderthal admixture, or the ancestors of basal Eurasians lived in the Near East, but did not participate in the Neanderthal admixture. Our results seem to fit the first scenario—basal Eurasians from eastern Arabia mixing with Near Easterns before the Natufian period—although they do not exclude other hypotheses. Thus, the exposed basin of the Arabo-Persian Gulf is a possible homeland of basal Eurasians, with an easy corridor linking the current Hormuz Strait to the Zagros and Caucasus steppes (the east continuum). Recent archeological evidence (Heydari-Guran and Ghasidian 2020) from the Zagros region indicates that this region was passable in the Middle and Upper Paleolithic and included intermountain plains connected to each other by valleys, associated with permanent water and raw material sources. If the basal Eurasians were located in the Arabo-Persian Gulf, a sister population would have mingled at a northern location, around Levant, with the autochthonous Neanderthals, and were the ancestors of the current European and Asians pools. The basal Eurasians and the Neanderthal admixed group were genetically close, so they most likely descended from the same African migrant group that had split somewhere earlier. The split might have occurred in southwestern Asia after the OOA migration (through either the northern or the southern route; Lahr and Foley 1994), or, alternatively, in Africa. In the latter scenario, the subset that gave rise to the basal Eurasian branch probably followed the southern route taking refugium in the exposed basin of the Arabo-Persian Gulf, whereas the direct ancestors of Europeans and Asians followed the northern route, mixed with Neanderthals, and hence moved forward, further splitting towards Europe and Asia. Current evidence does not allow us to disentangle between the two scenarios, which highlights the urgency of finding and analyzing ancient human specimens in eastern AP/Zagros region. These specimens do not need to be as old as 50 ka, 20–10 ka specimens should expectably bear a higher basal Eurasian ancestry and might represent excellent proxies for the demographic, admixture and dispersal events that occurred on the awake of OOA and around the split time of the main population groups existing nowadays.

Materials and Methods

Data Sets

The “modern data set” included 291 individuals from AP (100 Saudi Arabia, 61 Yemen, 59 Oman and 71 UAE) and 78 from Iranian populations previously genotyped with the Illumina Human 741,000 SNPs Omni Express Bead Chip (OmniExpress; Fernandes et al. 2019) and deposited in EGA (Archive with accession number: EGAS00001003335). These individuals are nonrelated and resided in Dubai at the time of sampling (between years 1991 and 2000), representing a widespread collection of birth origins from Saudi Arabia, Yemen, Oman, the UAE, and Iran. We further confirmed that they were nonrecent migrants from sub-Saharan Africa from the genomic data. This data set was merged with 1,998 worldwide relevant available samples (supplementary table S7, Supplementary Material online) through PLINK 1.9 (Chang et al. 2015). Minimum values of 2% were set for missing genotype per marker and missing genotype per sample, leading to a final set of 291,595 SNPs and 2,367 samples. Genotypes from the “modern data set” were phased with SHAPEITv2.r79044 (Delaneau et al. 2011) using the 1,000 Genomes phased data (Auton et al. 2015) as a reference panel and using the HapMap b37 genetic map (Frazer et al. 2007). For algorithms relying on independent markers, SNPs from the modern data set were pruned for LD, by removing any SNP that had an r² > 0.2 with another SNP, within a 50-SNPs sliding window with steps of 5 SNPs. After pruning, the “modern data set pruned” included 116,084 autosomal SNPs and 2,367 individuals. The ancient data set included 304 samples from several recent studies (downloaded from David Reich Lab website; https://reich.hms.harvard.edu/datasets; downloaded on April 11, 2019) as detailed in supplementary table S8, Supplementary Material online. The two data sets were merged in Ancient-Modern Data set (“AMD”) and checked for removal of SNPs with more than 50% of missing genotype (geno 0.5) and samples with less than 10,000 SNPs, resulting in 293,034 SNPs (or 57,338 SNPs when pruned) and 2,642 samples (275 ancient and 2,367 modern samples).

Population Structure and Differentiation

To visualize the genetic similarities of ancient and modern samples, PCA was performed in smartpca implemented in EIGENSOFT (Patterson et al. 2006; Price et al. 2006) on the AMD. The ancient samples were projected on the modern data set using the parameters “lsqproject” and “shrinkmode” options turned on.

To analyze the genetic relationships between modern and ancient samples, the ADMIXTURE software (Alexander et al. 2009) was applied on the “modern data set pruned” (but limited to the 57,338 SNPs), using maximum likelihood for components (K) from 2 to 20, with the optimal K estimated through CV of the logistic regression. The ancient samples were projected on the previously inferred ancestral allele frequencies.

Population structure was also evaluated on the “modern data set” based on haplotype information through CHROMOPAINTERv2 (CP) and fineSTRUCTURE v2.07 (Lawson et al. 2012), in which the modern samples were used both as donors and receivers (all vs. all method). The Ne is the “recombination scaling constant” and is directly related to the effective population size. It is used by CHROMOPAINTER to convert the values of the genetic distance between SNPs (taken from the genetic map) to the population-scaled values for these distances required by the algorithm. The θ is the “per site mutation rate parameter” and is used by the CHROMOPAINTER Hidden Markov Model (HMM) to allow for imperfect copying between haplotypes. The N_e and θ parameters were inferred performing 10 iterations on 6 randomly selected chromosomes [1, 6, 10, 15, 18, and 22] and around 10% of the total samples from each population making a total of 239 samples as in Busby et al. (2015) using the expectation–maximization (EM) option of CHROMOPAINTERv2.

RFMix software (Maples et al. 2013), was used in the modern AP and Iranian samples to infer ancestry for each segment of the genome between a mixture of putative ancestral panels of haplotypes. 1,000 Genomes populations were used as parental: Great Britain, representing European ancestry; Yoruba from Nigeria, as the sub-Saharan African ancestry; and Indian Telugu for the South Asian ancestry. Then the median sizes (in cM) of the sub-Saharan African blocks were estimated per AP and Iranian individuals affiliated in clusters A and D (supplementary fig. S4, Supplementary Material online).

Modeling the Ancestry Composition of AP

QpAdm tool from the ADMIXTOOLS software package (Patterson et al. 2012), which relates a set of “left” populations (the population of interest and potential sources of ancestry) to a set of “right” populations (diverse outgroups, supplementary table S9, Supplementary Material online), was used to test for the number of streams of ancestry from “right” to “left” and estimate mixture proportions. All the possible combinations of N = (2...5) and one set of right/left outgroups were evaluated, as done in Lazaridis et al. (2016). A P-value threshold of 0.01 was set and only feasible mixture proportions were reported.

The N_e and θ values obtained in the CP analysis of the modern data set were used in an “unlinked” CP analysis with the adaptation of the NNLS (CP/NNLS) function as in Leslie et al. (2015) and Montinaro et al. (2015), by adding ancient samples informative for this geographic region (and having the least number of missing genotypes). CP/NNLS estimates the proportions of the genetic contributions from ancient periods to the modern clusters (Raveane et al. 2019). For every tested model (per period—Paleolithic, Neolithic, Bronze Age; overall per country and per cluster—limiting to the four sources identified in the qpAdm results), the standard errors were estimated through a weighted jackknife bootstrap, by removing one chromosome at a time, and averaging the values taking into account the total number of markers analyzed for iteration (Busing et al. 1999). In order to evaluate the fitness of the CP/NNLS estimations, the sum of the squared residuals were inferred and the ones with the lowest values were selected as best estimations for the tested models (Raveane et al. 2019).

Impact of the Basal Eurasian Ancestry

To address the impact of the basal Eurasian lineage in AP (Lazaridis et al. 2014; Petr et al. 2019), the statistic-f4 implemented in the ADMIXTOOLS package (Patterson et al. 2012) was applied to the modern samples organized per countries and per cluster. The ∼45 ka Ust’-Ishim Upper Paleolithic Siberian (Fu et al. 2014), who is a proxy for nonbasal Eurasian ancestry (Lazaridis et al. 2016), was considered in the test f4(Modern population Test, Han; Ust’-Ishim, Mbuti), as performed in Fu et al. (2016). Significant negative values (Z < −3) are evidence for the test populations having ancestry from the basal Eurasian lineage.

Additionally, the qpAdm tool was also used to search for traces of the basal Eurasian lineage in modern AP populations. The tested model was in the form (Test, Mota, an early HG-represented by NWRussia—I0061), and one set of right outgroups (supplementary table S9, Supplementary Material online) as in Lazaridis et al. (2016). Mota is an ancient eastern African (4.5 ka; Llorente et al. 2015) who, similarly to the basal Eurasians, is basal to other non-Africans, hence indirectly providing an estimate of the mixture proportions of the basal Eurasian lineage. The 7 ka Karelia sample (Haak et al. 2015) was used as early HG. To relate the basal Eurasian lineage and Neanderthal proportions, the Altai Neanderthal, Denisovan, and Mbuti samples were extracted and merged with the modern phased data set (276,035 SNPs). We then performed a f4-statistics analysis in the form f4(Test, Mbuti; Altai, Denisova) as in Lazaridis et al. (2016) using the ratio implemented in the ADMIXTOOLS package qpDstat (Patterson et al. 2012), with the f4 mode = yes, indicating that f4 statistics not D-stats are computed. The two proportions were correlated as in Lazaridis et al. (2016), through a linear regression using the lm function in R.

Additional f4 tests were further conducted in the form of f4(Test, Han; Ust’-Ishim, X), in which X represented ancient samples meant to rule out the possibility of basal Eurasian ancestry deriving from recent admixture from a nonbasal Eurasian population. These ancient samples were: the Caucasus HG Kotias (9 ka) or Satsurblia (13 ka) (Jones et al. 2015); ∼13 ka Natufian farmer representative of the Levant input (Lazaridis et al. 2016); ∼9 ka Zagros farmer (sample ID WC1.SG from Broushaki et al. 2016); ∼15 ka North Africa Iberomaurusian (sample TAF011 from van de Loosdrecht et al. 2018).

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Acknowledgments

This work was financed by FEDER-Fundo Europeu de Desenvolvimento Regional funds through COMPETE 2020-Operacional Programme for Competitiveness and Internationalization (POCI), Portugal 2020, by Portuguese funds through FCT-Fundação para a Ciência e a Tecnologia, Ministério da Ciência, Tecnologia e Inovação in the framework of the project “Biomedical anthropological study in Arabian Peninsula based on high-throughput genomics” (POCI-01-0145-FEDER-016609), the Italian Ministry of Education, University and Research project Dipartimenti di Eccellenza Program (2018–2022)—Department of Biology and Biotechnology “L. Spallanzani,” University of Pavia (to A.T.). V.F. has a postdoc grant through FCT (SFRH/BPD/114927/2016). i3S is financed by FEDER-COMPETE 2020, Portugal 2020 and by Portuguese funds through FCT in the framework of the project “Institute for Research and Innovation in Health Sciences” (POCI-01-0145-FEDER-007274). Authors would like to thank Dr Francesco Bertolini for facilitating the research of A.R. in the last stage of the article preparation.

Data Availability

The data underlying this article were previously published in Fernandes et al. (2019) and deposited in EGA Archive at https://ega-archive.org and can be accessed with accession number: EGAS00001003335.