Genomic transformation and social organization during the Copper Age-Bronze Age transition in southern Iberia.

The emerging Bronze Age (BA) of southeastern Iberia saw marked social changes. Late Copper Age (CA) settlements were abandoned in favor of hilltop sites, and collective graves were largely replaced by single or double burials with often distinctive grave goods indirectly reflecting a hierarchical social organization, as exemplified by the BA El Argar group. We explored this transition from a genomic viewpoint by tripling the amount of data available for this period. Concomitant with the rise of El Argar starting ~2200 cal BCE, we observe a complete turnover of Y-chromosome lineages along with the arrival of steppe-related ancestry. This pattern is consistent with a founder effect in male lineages, supported by our finding that males shared more relatives at sites than females. However, simple two-source models do not find support in some El Argar groups, suggesting additional genetic contributions from the Mediterranean that could predate the BA.

INTRODUCTION

During the last centuries of the third millennium BCE, the societies of Europe, the Near East, and Egypt experienced large-scale social and political upheavals. Settlement abandonment, depopulation, the disappearance of communication networks, and major political disruptions at the end of the Akkadian empire and the Old Kingdom in Egypt have often been interpreted in the light of a climatic crisis, known as the 4.2k event (–). Recently, the possibility of substantial population movements, causing social instability during the third millennium BCE, has been proposed as a further explanation for the changes observed at the end of the Copper Age (CA) in central and western Europe (–). Signs of social and economic turnover are particularly marked in the southern half of the Iberian Peninsula (), where the CA is associated with exceptional demographic growth, a diversity of monumental settlements and funerary structures, widespread copper metallurgy, and a sophisticated, large-scale manufacture and exchange of symbolic goods, among others [e.g., (, )]. Moreover, this period is characterized by a diversity of settlement types, including fortified sites, ditched enclosures, and so-called megasites, some of which exceeded 100 ha in size (e.g., Valencina de la Concepción and Marroquíes Bajos) and all of which were formed at around 3300 to 2800 BCE, therefore predating the Bell Beaker horizon.

This period is also associated with a major increase in interconnectedness and mobility. On the basis of available radiogenic (Sr) isotope studies, the percentage of southern Iberian individuals who were buried in locations other than where they grew up ranges between 8 and 74% (, ). Ivory from Africa and the Near East (–), amber from Sicily (), and ostrich eggshells from Africa () are indicative of transregional connections. However, evidence of a strong political centralization and economic inequality remains elusive or inconclusive (–).

Archaeogenetics has suggested that the remarkable development during the (south) Iberian CA was coupled with a strong population continuity attested since the Neolithic [e.g., (, , , –)]. However, the Late CA anthropological and archaeological records from the north and central Iberia show the first individuals carrying “steppe-related ancestry” by ~2400 calibrated (cal) BCE, which are often but not exclusively linked to Bell Beaker–associated artifacts (, ). In parallel, African ancestry was also reported in one individual, which suggests discrete movement/mobility of people ().

The beginning of the Bronze Age [Early Bronze Age (EBA)] in Iberia (2200 to 1550 cal BCE) marks a clear population turnover, suggested by both the omnipresence of steppe-related ancestry in all individuals directly postdating 2200 BCE. An even more notable shift can be observed in the frequency of Y-chromosome haplogroups in males, which are almost exclusively of the R1b-P312 type that was completely absent in Iberia before 2400 BCE (, , , ).

The turn from Late CA to the EBA at the end of the third millennium BCE saw the demise of fortified settlements such as Los Millares and ditched megasites such as Valencina and Perdigões in southern Iberia, while in the southeast it is concomitant with the rise of new hilltop occupations of smaller sizes (<0.5 ha) (, , ). Substantial and densely built hilltop settlements, distinguished by a specific intramural burial rite and characteristic ceramic and metal types, appear around 2200 BCE in the fertile tertiary basins framed by mountain ranges, running parallel to the coast of southeastern Iberia (). This area of circa (c.) 3500 km2 is considered to be the core of the El Argar “culture,” which is one of the most outstanding examples of an early complex society in European prehistory with evidence for social stratification (, , ). The origins of El Argar are still unclear, as there are no hybrid contexts where El Argar elements appear in Late CA settlements or vice versa. Although the early El Argar material record shares some traits with the Bell Beaker complex (), such as V-perforated buttons, Palmela-type points, or perforated stone plaques so-called “archers’ wristguards,” the characteristic Bell Beaker pottery is absent. Upon the discovery of the monumental fortification in the 5-ha hilltop settlement of La Bastida, dated to around 2200 BCE, a possible eastern Mediterranean contribution was reconsidered (). The intramural burials in large storage vessels (pithoi), the circulation of silver rings and bracelets, and the characteristic footed Argar cup have also been interpreted as signs of Aegean or Near Eastern contacts (), although all these features emerged during later phases of El Argar. The genealogy of different characteristic material traits of early El Argar, such as apsidal buildings, intramural burial, metal casting technology, and the halberds as distinguished weapons, is reminiscent of several social developments in southeastern, central, and western Europe, with possible links of still uncertain origin ().

Between 2000 and 1800/1750 BCE, El Argar expanded across a wider region of southeastern Iberia and entered the central Spanish Meseta. Characteristic El Argar items, such as halberds, are also present beyond this territory. The leading figures within El Argar society seem to have been a class of warriors armed with halberds and daggers. These weapons, occasionally associated with golden arm rings, also became the insignia of political domination in the male elite burials of central Europe around 2000 to 1800 BCE, a time when this specialized weapon had already gone out of use in the rest of Europe (). At the same time, a growing part of the population, particularly children, was buried in cists, artificial caves, pottery vessels, or pits inside the distinctive settlements of El Argar.

During its final phase (1800/1750 to 1550 BCE), El Argar reached an outstanding level of economic and social development. A vast amount of grinding tools, large workshops, and storage facilities found in the larger hilltop settlements (1 to 6 ha) imply that certain groups managed to control the flow of resources and the workforce of wider regions (). The establishment of a dominant and hereditary class, as well as increasing social asymmetries, become recognizable in the intramural burials, monumental architecture (, ), and the spatial relation between both (, , ).

Lull and colleagues (, ) have argued that social conflict and environmental degradation resulting from deforestation and extensive dry farming could have led to the abandonment or destruction of all El Argar settlements around 1550 BCE. However, the rise of a similar economic organization, architecture, and funerary record in inner Alicante, as highlighted by Cabezo Redondo, suggests that during the peak of El Argar at least some groups managed to establish themselves in a territory that was influenced by the Valencian BA, a neighboring “cultural group” located in southeastern Iberia.

In this study, we aimed to understand the importance of population dynamics in the collapse of the highly dynamic Iberian CA world, the rise and development of El Argar, as well as its relations between neighboring BA groups in western Europe and the Balearic Late BA (LBA). We also explore the genetic makeup of the BA groups in Iberia and the Balearic Islands in relation to the other western and central Mediterranean islands, such as Sardinia and Sicily. In total, we characterize the genomic profiles of 96 BA individuals associated with El Argar and contemporaneous societies, as well as 34 CA and 6 LBA individuals.

RESULTS

Here, we report genome-wide data in the form of 1.24 million ancestry-informative single-nucleotide polymorphisms (SNPs) for 136 southern Iberian individuals, covering a time span of 2000 years from the Late Neolithic (LN)/CA (3300 cal BCE) to the LBA (1200/1000 cal BCE) (Fig. 1, table S1.1, and text S1). Five of these individuals were also shotgun-sequenced to 1.5× to 5.2× coverage to reconstruct the complete genomes to add to the growing record of ancient western Eurasian genomes for future studies. Our new dataset includes various groups and types of sites such as sepulchral caves from the LN/CA (Cova d’en Pardo, n = 7; Cueva de las Lechuzas, n = 10), simple pit burials of pre–Bell Beaker CA megasites [Valencina (PP4-Montelirio), n = 11], a collective hypogeum with Late CA burials (Camino del Molino, n = 6), and EBA pit graves (Molinos de Papel, n = 3) to provide a diverse comparative dataset for the Argaric societies, which are at the center of our study. For the latter, we include the almost complete and archaeologically well-defined sites of La Almoloya (n = 67) and La Bastida (n = 10), alongside Cerro del Morrón (n = 3) and Lorca town [Madres Mercedarias Church (n = 1), Los Tintes (n = 1), and Zapatería (n = 1) sectors]. We also analyzed individuals of the so-called Valencian BA, including Cabezo Redondo (n = 1), Peñón de la Zorra (n = 1), Puntal de los Carniceros (n = 4), and La Horna (n = 3), individuals of the BA from Catalonia (Miquel Vives, n = 1), as well as LBA individuals from the Balearic Islands of Menorca (Es Forat de ses Aritges, n = 6) (Fig. 1, table S1.1, and text S1).

Fig. 1.

Geographic locations as well as cultural and chronological information of the studied sites.

Map of Iberia with sites mentioned in the main text (table S1.1). The area shaded in teal highlights the maximum extent of EBA El Argar. Chronological time scale for published and new individuals analyzed in this study (bottom panel). Directly radiocarbon-dated individuals are plotted according to their mean calibrated date (2-sigma range), and jitter option within their specific time range was applied for individuals that were dated by archaeological context (table S2.1). Random jitter was applied to the Y axis to prevent overplotting.

We initially screened 244 individuals for DNA preservation by shotgun sequencing ~5 million Illumina single end reads per partial uracil-DNA glycosylase (UDG)–treated DNA library (), followed by assessment of % endogenous human DNA (>0.1%), average read length between 40 and 75 base pairs (bp), and characteristic ancient DNA (aDNA) damage patterns >3% (tables S1.1 and S1.2). Libraries fulfilling these criteria were retained for hybridization capture of 1.24 million informative sites (1240k SNP panel) in the human genome (). Libraries were independently captured for mitochondrial genomes (tables S1.3 and S1.4 and text S2) [() modified following ()] and mappable regions of the Y chromosome (table S1.1 and text S3) (). Following standard aDNA processing pipelines (), we quantified contamination rates at the autosomal (in males) () and mitochondrial levels (both genetic sexes) (), showing low contamination estimates for the individuals in this study (<3% for nuclear and <1% for mitochondrial DNA; see Materials and Methods). By determining the genetic sex (), we observed two cases of sex chromosome aneuploidies, one XXY (Klinefelter) individual and one XXX (trisomy X) individual (table S1.1, text S4, and fig. S1). We obtained endogenous human aDNA contents ranging from 1.3 to 29.5% on the targeted 1240k SNPs, which correspond to 2004 to 831,086 SNP positions covered (tables S1.1 and S1.2). For population genetic analyses, we only retained individuals with more than 40,000 SNPs on the 1240k panel, excluding 18 low-coverage individuals in downstream analyses. We performed kinship analysis (–) for all newly typed individuals (fig. S2; tables S1.5, S1.6, and S1.7; and text S5), the result of which will be featured in a separate study. From each pair of first-degree relatives, we then excluded the individual with lower SNP coverage. In summary, for in-depth population genetic analysis, we merged genome-wide data from 122 newly typed individuals with a comparative set of previously published ancient and modern-day individuals, which are reported by the Reich Lab (https://reich.hms.harvard.edu/datasets; please see the Supplementary Materials for a detailed list of references), as well as genotypes from recently published studies (table S2.1) (–).

Genetic substructure in the Iberian CA

We first explored the genetic affinities of the newly typed CA individuals by performing a principal components analysis (PCA) with relevant ancient individuals projected onto PCs calculated from a set of modern-day West Eurasian populations genotyped with the Human Origins (HO) SNP panel (Fig. 2A) (). The new CA individuals from southern Iberia fall onto a position that partially overlaps with previous Middle Neolithic (MN), Middle/Late Neolithic (MLN), and CA (non-steppe) groups from Iberia but are slightly positively shifted in their coordinates for PC1 toward previously published Early Neolithic (EN) groups from Iberia and later groups such as Sardinia Chalcolithic, suggesting an equally small hunter-gatherer (HG) ancestry contribution in the CA individuals of southern Iberia. These results run counter to previous suggestions of steppe-related ancestry diffusion into southern Iberia in the early CA based on supposed formal similarities in one element of the Valencina material culture (). To formally test the differential contribution of the main HG source in European Neolithic farmers, called Western HG (WHG), which represents the predominant form of ancestry in prefarming western Europe (, ), we calculated f4-statistics of the form f4(SE/SW_Iberia_CA, N/NE/C_Iberia_CA; WHG, Mbuti) and obtained a significantly negative value for N_Iberia_CA and C_Iberia_CA (Fig. 2B and table S2.2), indicative of higher WHG ancestry in these regions, a signal that is also present in the preceding MLN period, as shown by f4-statistics of the form f4(SE/SW_Iberia_MLN, N_Iberia_MLN; WHG, Mbuti) (table S2.2).

Fig. 2.

Key population genetic analyses of CA groups.

(A) PCA of published and newly genotyped CA groups projected onto 2018 modern-day West Eurasians (gray dots). Ancient individuals projected and their correspondent labels are listed in table S2.1. (B) f4-statistics showing significant differences in terms of WHG ancestry in northern and southern CA groups (error bars indicate 3 standard errors) (table S2.2 also includes the results for MLN groups). (C) f4-statistics highlighting the higher affinity of southern CA Iberian groups to GoyetQ2 (error bars indicate 3 standard errors) (table S2.3). (D) Modeling Iberian MLN and CA without steppe-related ancestry with a three-way qpAdm admixture model using distal sources Anatolia_N, WHG, and GoyetQ2. Southern MLN and CA Iberian groups show an extra minor ancestry component represented by GoyetQ2 (error bars indicate 1 standard error; numbers in brackets are P values for qpAdm model) (table S2.4).

Leveraging insights from previous studies about the differential HG legacy in western Europe (, –), we explored the HG ancestry of geographically diverse CA Iberians in different ways. We first calculated f4-statistics of the form f4(GoyetQ2, WHG; test, Mbuti), where test represents all MLN and CA groups (Fig. 2C and table S2.3). Despite the resulting f4-values being significantly negative for all test groups (indicating WHG admixture), we observed a geographically related gradient with southern groups (MLN and CA) showing lower negative values, indicating a differential relationship with GoyetQ2, who acts as a proxy for Magdalenian-related ancestry (Fig. 2C). We confirmed this subtle signal using qpAdm outgroup-based ancestry modeling with Anatolia_Neolithic, WHG (Loschbour and KO1), EHG [Eastern European HGs; the predominant form of ancestry in pre-agropastoralist eastern Europe; ()], and GoyetQ2 as distal sources (table S2.4). Here, we found that the CA groups from southern Iberia differ with respect to the overall quantity of HG ancestry, but the latter is also qualitatively different. Estimates of Magdalenian-associated HG ancestry in southern Iberian CA individuals vary from 6.1 ± 1.3% to 7.2 ± 1.3%, which reflects the geographic structure of local HG groups (, ) and suggests a certain amount of genetic continuity since the Neolithic (Fig. 2D and table S2.4). In the case of SE_Iberia_CA, adding GoyetQ2 as a third source improves the model fit slightly from P = 2.53 × 10−6 to P = 0.008. Adding Iran_N or Jordan_PPNB as a fourth source (Anatolia, WHG, GoyetQ2, and Iran_N/Jordan_PPNB) further improves the model fit for SE_Iberia_CA by an order of magnitude, i.e., from P = 0.008 to P = 0.014 (with Iran_N) or P = 0.027 (with Jordan_PPNB) (table S2.4, text S7, and fig. S3). In turn, adding other populations as a fourth source does not improve the model fit (text S7). The unknown source contributing to SE_Iberia_CA groups is likely to have carried an excess of Levantine and/or Iran_N-like ancestries compared to the distal source Anatolia_Neolithic, as these together have been found admixed in Anatolian and Levantine CA groups (~6000 to 5000 BCE) (). This finding suggests a subtle contribution that was spread early along the Mediterranean or, alternatively, different sources of early farmer ancestry during the Neolithic transition with varying proportions of Levantine and/or Iran_N-like components when compared to Anatolia_N used here. Removing Anatolian HG (AHG) from the outgroups also improves the model fit (P = 0.045; table S2.4), indicating that the Neolithic ancestry is not well represented by using Anatolia_N as a distal proxy and might come from another farmer group more similar to AHG than Anatolia_N (text S7). More individuals from the Neolithic and CA across the Mediterranean will be needed to track this contribution more confidently.

The Magdalenian-related ancestry was not detectable in other contemporaneous Mediterranean populations (e.g., Sardinia_Neolithic, Sardinia_Chalcolithic, Sardinia_EBA, Sicily_EBA, and Italy_CA), which renders gene flow in the direction from southern Iberia to the central Mediterranean less plausible. However, applying the same qpAdm model (text S7, fig. S3, and table S2.5), we detect a previously unreported amount of Iran_N-like ancestry in central Mediterranean groups from Sardinia, ranging from 2.8 ± 1.2% in Sardinia_Chalcolithic to 5.8 ± 1% in Sardinia_Nuragic_BA. Adding Iran_N as a source to model Sicily_EBA improves the model fit but without reaching P values ≥0.05, making gene flow from the Italian Peninsula to Sardinia more likely (Italy_CA also shows Iran_N-like ancestry) than from Sicily (Sicily_EBA) (text S7, fig. S3, and table S2.5). Notably, when modeling Sicily_EBA, we obtain P values ≥0.05 by removing AHG from the outgroups in a three-source model (Anatolia_Neolithic, WHG, and Yamnaya_Samara) or by removing Morocco_Iberomaurusian in a four-source model (Anatolia_Neolithic, WHG, Yamnaya_Samara, and Iran_N), which points to genetic substructure in the Mediterranean before the arrival of steppe-related ancestry. However, this finding does not rule out limited genetic input from other Mediterranean populations to southeastern CA Iberians (text S7).

Genetic turnover in the southern Iberian BA and the rise of El Argar

When compared to the preceding CA groups, individuals associated with El Argar BA and SE_Iberia_BA form a dense cloud in PC space that is shifted toward populations with steppe-related ancestry from central Europe and positioned between individuals from the Iberian CA and those of the north Iberian BA (Fig. 3A). A similar genetic shift is noticeable in the ADMIXTURE () results (text S8 and fig. S4), where Iberian BA individuals show an additional component when compared to preceding CA individuals, and this is further supported by an f4-test of the form f4 (Anatolia_Neolithic, test; Yamnaya_Samara, Mbuti), in which test iterates through newly reported CA and BA individuals (Fig. 3B and table S2.6). Here, the shift to negative f4-values in EBA individuals indicates shared drift with Yamnaya_Samara individuals, which is absent in earlier CA individuals.

Fig. 3.

Evidence for steppe-related ancestry in EBA individuals from southeastern Iberia.

(A) West Eurasian PCA calculated with modern populations (gray dots) () on which relevant ancient CA, Bell Beaker and BA groups discussed in the text were projected (correspondent labels are listed in table S2.1).(B) f4-statistics showing the increased affinity to Yamnaya_Samara (absent in CA individuals), which implies the presence of steppe-related ancestry in EBA Iberians. This plot also shows the almost-complete turnover in Y-chromosome lineages in male individuals (color-filled squares) during the southeastern Iberian EBA (error bars indicate 3 standard errors) (table S2.6).

This observation suggests a substantial amount of steppe-related ancestry in El Argar BA individuals, which we tested formally and directly with f4-statistics of the form f4(Argar_Iberia_BA/SE_Iberia_BA, SE_Iberia_CA; Yamnaya_Samara, Mbuti) (fig. S5A and table S2.7). Significantly positive f4-values confirmed the presence of steppe-related ancestry in all BA individuals. We then tested for differences in affinity to steppe-related ancestry by contrasting northern versus southern BA individuals using f4(N/NE/C_Iberia_BA, Argar_Iberia_BA/SE_Iberia_BA; Yamnaya_Samara, Mbuti) (fig. S5B and table S2.8). The resulting f4-values confirmed a smaller amount of steppe-related ancestry in individuals from the Argaric sites La Almoloya and La Bastida compared to the rest of Iberia_BA groups, especially when compared to those from northern Iberia (fig. S5B and table S2.8), despite the complete turnover to lineage R1b-P312 (except for one subadult male in La Bastida) visible in the Y-chromosome record (Fig. 3B, table S2.6, and text S8). However, at the intrasite level, we observe no significant differences with respect to the amount of steppe-related ancestry between the early and late phase of La Almoloya and La Bastida based on PCA and formal f4-statistics (Fig. 3A and fig. S5), which suggests that the contribution is homogenized across the population.

The newly analyzed LBA individuals from the Balearic Islands (Aritges_LBA) showed less steppe-related ancestry than the published Late CA individual from Mallorca (previously named Mallorca_EBA), confirming a decrease of steppe-related ancestry over time in agreement with (). The amount of steppe-related ancestry in Mallorca_CA_Stp is similar to the first CA groups in Iberia (C_Iberia_CA_Stp and NW_Iberia_CA_Stp) (fig. S5A) and could have arrived at the same time, as there is no clear evidence of human occupation in the archipelago before their arrival (, ).

Olalde and colleagues () had proposed for the Iberian Peninsula that the first contribution of steppe-related ancestry was diluted during the BA by admixture with descendants of local CA groups, but increased again during the LBA–Iron Age during a second pulse. We do not detect this second pulse in the Balearic Islands, in agreement with the findings of (), albeit based on one LBA individual from Menorca (fig. S5, A and B).

The situation observed in southern Iberia offers a specific example of the demographic changes occurring during the transition from CA to EBA societies around 2200 BCE. The absence of steppe-related ancestry in our newly typed individuals from the Early CA at Valencina (c. 2900 to 2800 BCE) and the Late CA collective burial cave Camino del Molino, which is contemporaneous with the Beaker Horizon (albeit without Beaker pottery) and the first appearance of steppe-related ancestry in a double tomb with Bell Beaker–associated artifacts at Molinos de Papel, a site only 400 m away from the latter, provides the temporal framework for the arrival of steppe-related ancestry in southern Iberia. Taking the youngest date from Camino del Molino [Beta-261524, 3830 ± 40 years before present (BP)] and the oldest from Molinos de Papel (MAMS-11826, 3780 ± 30 BP) as reference points, we can show that both C14 dates overlap at the 95% level and date the arrival of steppe-related ancestry in southeastern Iberia to ~2200 cal BCE. At a cross-regional scale, none of the south Iberian CA individuals analyzed so far show evidence of any steppe-related ancestry, whereas the oldest individuals from the El Argar sites of La Bastida (BAS024 and BAS025) and La Almoloya (ALM019), directly dated to the 21st century BCE based on the radiocarbon evidence, show clear evidence of such ancestry. Consequently, the population mixture must have occurred in the preceding 150 years at the latest. Estimating the time of admixture of the putative sources using the program DATES (Distribution of Ancestry Tracts of Evolutionary Signals; https://github.com/priyamoorjani/DATES) produced compatible admixture dates in all new BA groups, ranging from ~2550 BCE ± 189 years to ~1642 BCE ± 133 years (table S2.9 and text S8). This could have occurred after the collapse of the CA social system in some aggregation megasites such as Valencina at around 2350 BCE () and more or less contemporary in others, such as Los Millares at around 2300/2200 BCE (). The demise of CA cultures ~2350 to 2200 cal BCE has been linked to the so-called 4.2–thousand year cal BP climatic event, as there is a temporal overlap between this event and the collapse of large megasites (, –), although direct, conclusive, paleoclimatic evidence for Iberia remains sparse (, ). The social and economic changes that lead into the BA, as well as the arrival of the new genetic ancestry, appear to have spread from North to South Iberia (, ) and thus might be linked. However, it remains unclear whether this expansion was opportunistic (i.e., followed the consequences of a potential climatic deterioration) or causal to the actual shifts observed in archaeology. Notably, the arrival of new genetic ancestry is not paralleled by the same social changes in all regions of Iberia. In the coastal regions of southeastern Iberia, the shift to single or double burial practices inside settlements occurred ~2200 BCE, coinciding with the beginning of El Argar. However, in inner Iberia, the associated grave goods and settlement pottery are more related to late Bell Beaker CA (leaf-shaped tanged point, V-perforated button, and incised Beaker pottery) than to the material culture of the early El Argar “core” territory (). In La Mancha, the southern part of the central Spanish plateau, the first individual with a confirmed steppe-ancestry component (), a male from tomb 4 of Castillejo del Bonete dated to ~2100 BCE, coincides with the foundation of monumental settlements, such as Las Motillas fortifications.

This shows that while the genetic contribution of steppe-related ancestry to Iberia was a long-term process starting around 2400 cal BCE in the northern and central regions (), from where it spread southward over ~300 years. At a local scale, this change might have occurred faster. A similar situation might have existed in central Portugal, where we still find individuals with no steppe-related ancestry in collective CA burials (Galería da Cisterna and Cova da Moura) around 2300 to 2200 cal BCE. However, after 2100 cal BCE, all individuals from all sites carry steppe-related ancestry, in line with R1b-P312 becoming the predominant Y-chromosomal lineage present not only in El Argar but also in the rest of BA Iberia.

Mediterranean and central European ancestries shaped the genetic profile of southeastern BA groups in Iberia

To explore the genetic turnover and the contribution of the local groups to the newly formed BA genetic profile in Iberia, we systematically tested a series of qpAdm models. We started by using the distal ancestry sources Anatolia_N, WHG, GoyetQ2, Yamnaya_Samara, and Iran_N to model the genetic ancestry components of Iberian BA groups (table S2.10 and fig. S6). We found that the local traces of GoyetQ2, a characteristic but variable component of southern Iberia CA individuals, were no longer detectable, suggesting a dissolution of geographic substructure in BA Iberia with respect to HG ancestry. We explain this by the spread of steppe-related ancestry from North to South () that also contributed northern and central Iberian ancestry to the South, diluting the subtle GoyetQ2 signal to a level beyond the limits of detection (text S8). By using the same qpAdm model, we also observed that Almoloya_Argar_Early, Almoloya_Argar_Late, SE_CabezoRedondo_BA, and Bastida_Argar cannot be modeled with Yamnaya_Samara as a single source but find better support with a combination of Iran_N and Yamnaya_Samara, however, without reaching P values ≥0.05 in Almoloya_Argar_Early and Late and SE_CabezoRedondo_BA (table S2.10 and fig. S6). These El Argar groups (Almoloya and Bastida) are also slightly shifted to the right on the PC1 axis, in the direction of Mediterranean BA groups with excess Iran_N-like ancestry, such as “Minoans,” who only carry Iran_N-like ancestry but not steppe-related ancestry, or “Mycenaeans,” who carry a mix of both (), and that has also been shown for some BA individuals from Sicily_MBA () and for Sardinians here (Fig. 3A).

To explore the reasons for the model rejection observed in southeastern BA groups, we tested several qpAdm models (text S8). We used two- and three-way competitive qpAdm models to test whether a third source with Iran_N-like ancestry was needed. In the two-way model (table S2.11), Germany_Bell_Beaker (the largest number of individuals representing BB) was used as a fixed proximal ancestry source together with a local CA source (either N_Iberia_CA, C_Iberia_CA, or SE_Iberia_CA). For target groups in which these two sources were rejected, we iteratively tested central and eastern Mediterranean populations as a third source, as some of them are known to carry Iran_N-like ancestry (Fig. 4A and table S2.11). Using C_Iberia_CA as a local source of ancestry and Germany_Bell_Beaker as a proxy for steppe ancestry, only models with Almoloya_Argar_Early and Late as a target were rejected, suggesting that a third component was required for these Iberian BA groups. We found that adding any Mediterranean population as a third source improves the model fit, albeit without reaching P values >0.05. However, we obtained P values >0.05 when using Iran_Chalcolithic as a third source, which suggests that a higher proportion of Iran_N-like ancestry is needed (text S8 and table S2.11). Notably, when we exchanged C_Iberia_CA with the local SE_Iberia_CA, the model did not find support for a larger number of Iberian BA groups (N_Iberia_BA, SE_Iberia_BA, Almoloya_Argar_Early and Late, and C_Iberia_CA_Stp), and adding a third source in this constellation did not improve the model fit either (Fig. 4A and table S2.11).

Fig. 4.

Modeling genetic ancestry in Iberian BA groups.

(A) Modeling Iberian BA individuals with steppe-related ancestry as a two- and three-way qpAdm admixture model using proximal sources C_Iberia_CA/SE_Iberia_CA, Germany_Bell_Beaker, and Iran_C (table S2.11). (B) Modeling Iberian BA individuals with steppe-related ancestry as a two- and three-way qpAdm admixture model using proximal sources C_Iberia_CA_Stp, C_Iberia_CA/SE_Iberia_CA, and Iran_N (table S2.1). P values for each group are given inside each column. Faded colors indicate rejected models when applying a P value cutoff of ≤0.05.

To rule out that the rejection of models involving Germany_Bell_Beaker + C_Iberia_CA in La Almoloya (Early and Late) is due to specific north/central European LN admixture in Germany_Bell_Beaker that was not present in the source of steppe-related ancestry in Iberia, we fixed the source of direct steppe-related ancestry to Yamnaya_Samara, as we expect the same distal steppe source contributing to all Iberian groups, and iterated through the proximal local CA sources, which in these models contribute 100% of the Chalcolithic (i.e., Neolithic farmer) ancestry. We added distal Iran_N (instead of Iran_Chalcolithic, which could contribute other confounding CA ancestries) as a third source only when needed (table S2.12 and fig. S7). By using C_Iberia_CA and Yamnaya_Samara in a two-way model, we could successfully model all Iberian BA groups, but again to the exclusion of the southeastern Almoloya_Argar_Early and Late, and SE_CabezoRedondo_BA groups. Notably, Bastida_Argar also failed for the distal model 2 (table S2.10). However, these three groups returned values ≥0.05 in the proximal local CA substrate model when Iran_N was added as a third source (fig. S7 and table S2.12). In turn, when we exchanged C_Iberia_CA with a local SE_Iberian_CA source, most of the models were rejected, which confirms the genetic substructure in Iberia_CA groups. Using SE_Iberia_CA as substrate only improves the model fit for southeastern Iberian groups, suggesting that the first individuals who brought steppe-related ancestry to Iberia admixed locally with the different CA groups (table S2.12). Notably, SE_Iberia_CA is a sufficiently supported local proxy for Balearic CA_Stp, BA, and LBA groups. We interpret this signal with caution, as Mallorca_CA_Stp and Menorca_LBA are only represented by a single individual each; however, this signal is well supported in the newly typed Aritges_LBA group from Menorca, which consists of six individuals (fig. S7 and table S2.12).

We then tried to further specify the proximal sources without affecting the power of resolution to distinguish between them. We exchanged Yamnaya_Samara by C_Iberia_CA_Stp as the most proximal local source with steppe-related ancestry () and combined these with the different Iberian CA sources (table S2.13 and text S8) in two-way models, only adding Iran_N as a third source when needed. Using increasingly more possible proximal sources for steppe-related ancestry, we could replicate the results of the distal models. However, we can show that Almoloya_Argar_Early and Late, Bastida_Argar, SE_CabezoRedondo_BA, and Aritges_LBA still require Iran_N as a third source in all models with C_Iberia_CA and C_Iberia_CA_Stp. In contrast, Iran_N is not required when SE_Iberia_CA is used as a local substrate, although models for Almoloya_Argar_Early and Late are still rejected at P ≤ 0.05. This finding is best explained by either an unresolved distinctive ancestry of La Almoloya individuals or, alternatively, the ability to detect such subtle signals only in groups consisting of a larger number of individuals, as is the case for Almoloya_Argar_Early (n = 36) and Almoloya_Argar_Late (n = 22), which provide the statistical power to reject simpler models. When SE_Iberia_CA is used instead of C_Iberia_CA to model La Bastida, the model is no longer rejected, which suggests a direct contribution of local CA groups to those of the southeastern Iberian BA. Together, the results hint at an additional, albeit subtle, Iran_N-enriched ancestry that is present in SE_Iberia_CA groups and thus potentially predating the El Argar BA.

We also tested whether Almoloya_Argar_Early and Late can be modeled when Iran_N is exchanged by geographically more proximal Iran_N-rich sources involving central and eastern Mediterranean populations (table S2.14 and text S8) and by moving Iran_N to the outgroups. However, this did not result in any supported models (text S8). While we were not able to find a proximal Mediterranean source, we note that adding a central Mediterranean population to the outgroups (Sicily_EBA, Greece_EBA, or Greece_MBA) decreases the model support (P values) for Almoloya_Argar_Early and Almoloya_Argar_Late, indirectly attesting to the importance of central Mediterranean BA. However, we observed no changes in P values by adding other Mediterranean groups to the outgroups (text S8 and table S2.15).

The complete turnover on the Y chromosome to R1b-Z195 (a lineage derived from P312), observed in all males at La Almoloya (29 male individuals tested) and La Bastida (except for one E1b lineage dated around 2134 to 1947 cal BCE; 7 males tested), is another independent source of evidence of gene flow predominantly during the time of the CA-BA transition (table S1.1 and text S3). Notably, Y lineage R1b-Z195, the most common Y lineage in BA Iberia, ultimately derives from a common ancestor R1b-P312 in central Europe but already differs from other derived Bell Beaker R1b variants reported from central Europe and the British Isles. R1b-Z195 has been found in Sicily_CA_Stp (previously assigned to the EBA and thus considered outliers) and Sicily_EBA (). However, the genealogical and geographic link to other R1b variants still remains unclear. The subtle presence of Iran_N-like ancestry in El Argar and Y haplogroup R1b-Z195 in Sicily opens the possibility of gene flow not only from Iberia to Sicily as previously suggested () but also in the opposite direction, implying reciprocal contact with the western and central Mediterranean during the BA, although direct contacts are hardly noticeable in the archaeological record.

Together, these results suggest a dual genetic contribution to the formation of the BA genetic profile of southeastern Iberians in addition to a local CA genetic substrate. The major additional ancestry source resembles central European Bell Beaker groups, which first contributed ancestry to northern Iberia, followed by a southward spread in the form of C_Iberia_CA_Stp. A second minor ancestry component is an Iran_N-rich/central Mediterranean source, which is restricted to individuals from BA El Argar contexts. The timing of the last contributing component remains unclear and points to either a Neolithic legacy that persisted throughout the local CA or a subtle trace of connections to insular central Mediterranean BA groups.

A late Argar genetic outlier makes links to North Africa and the central Mediterranean

The PCA also revealed a genetic outlier among the newly typed individuals (Fig. 3A). The male individual ZAP002 from the late Argaric site of Lorca-Zapatería falls outside the Iberian_BA cluster on the PCA plot, with a clear shift toward the central Mediterranean. We formally tested for a specific attraction to groups from the eastern Mediterranean, the Near East, and Africa with an f4-admixture test of the form f4(ZAP002, Almoloya_Argar_Early; eastern_Med/Africa, Chimp) (Fig. 5A, table S2.16, and text S9). The resulting f4-values were consistent with zero, showing that ZAP002 and Almoloya_Argar_Early (the Argar group with the largest number of individuals) are symmetrically related to eastern Mediterranean and African groups. However, we noticed a near-significant positive f4-value with a z score = 2.4 for Morocco Iberomaurusian. The attraction to Mbuti in an analogous f4-test suggests either African ancestry or high levels of Basal Eurasian ancestry shared with Morocco Iberomaurusian and Natufians (table S2.16). We also confirmed a negative deviation from zero in f4(test, EHG; CHG, Mbuti), indicative of either African or Basal Eurasian ancestry. Here, we assume that the split of EHG and CHG groups [Caucasus HG; a basal form of HG ancestry south of the Caucasus ()] is sufficiently deep to result in negative values for any test population. Following this rationale, the finding of negative values suggests a high amount of shared deep ancestry with Mbuti, which also explains the different f4-value for ZAP002 in Fig. 3B when compared to other El Argar individuals. We observe negative values for ZAP002 (z score = −2.14), as well as Morocco Iberomaurusian and Natufian individuals, who provide additional support for this claim (Fig. 5B, table S2.17, and text S9).

Fig. 5.

Summary of genetic analyses that explore the genetic outlier status of individual ZAP002.

(A) f4-statistics to detect Basal Eurasian ancestry and/or African ancestry. Cladality test with ZAP002 and individuals from the site La Almoloya (ALM) and a set of test populations including Africans and Near Eastern (NE) groups, and chimpanzee in the outgroup. The results show that only Moroccan, Iberomaurusian, and Mbuti show a shift toward positive f4-values (table S2.16). (B) f4-statistics to detect either high Basal Eurasian or North African ancestry in which negative values imply excess affinity to Mbuti (table S2.17). (C) Supported and nonsupported qpAdm admixture models for outlier individual ZAP002. On the left, C_Iberia_CA is located in the outgroups and the model finds support with Sicily EBA as a source, pointing to a nonlocal origin of individual ZAP002 (table S2.18).

Using qpAdm, we further explored candidate sources of distinct ancestries such as Iberian CA and Germany_Bell_Beaker and a candidate list of central Mediterranean sources from Sardinia and Sicily. Using the same set of outgroups as for the main southeastern Iberian BA groups, we did not find a statistically well-supported model for ZAP002. However, moving Morocco_Iberomaurusian from the outgroups to the sources resulted in a model fit (≥0.05) when a central Mediterranean population was included as a source. We found that ZAP002 does not require a local Iberia_CA source and can be modeled as a three-way mixture of Germany_Bell_Beaker, Morocco_Iberomaurusian, and either Sicily EBA, Sardinia EBA, or Sardinia Nuragic BA (Fig. 5C and table S2.18). These qpAdm results are congruent with the location of ZAP002 in PC space, where he clusters with Sicily EBA individuals (Fig. 3A). Rotating the Iberian CA source to the outgroups in the same model still upholds the model fit (Fig. 5C, table S2.18, and text S9), suggesting an entirely nonlocal origin of individual ZAP002 despite being granted the same burial treatment as other individuals from El Argar (pithos burial and Argaric pottery) (Fig. 5C).

Insights into phenotypic variation, demography, and social correlates of CA and EBA El Argar societies

We explored the phenotypic variants included in the 1240k SNP panel among all well-preserved Iberian CA and BA individuals (>400,000 SNPs covered). In addition to 14 variants associated with the ability to digest lactose and alcohol, adaptation to metabolize fatty acids, predisposition to celiac disease, and resistance to infectious diseases, we also investigated 41 SNPs encoding for skin/hair and eye pigmentation (figs. S8 and S9, tables S2.19 and table S2.20, and text S10). We found no significant increase or decrease in allele frequency of any derived allele between the CA and BA (fig. S9 and table S2.20).

To gain insight into the effective population size during the CA and BA, we applied hapROH to identify runs of homozygosity (, ). We found no notable differences in short ROH (<2 cM) patterns between the two chronocultural periods, indicating that all of the individuals under study were drawn from populations with comparably and sufficiently large effective population sizes. We also found no evidence for a decline of CA population size before the emerging BA. In both groups, we did not detect signs of inbreeding, represented by long ROH segments (20 to 300 cM) (Fig. 6A and table S2.21).

Fig. 6.

Biological kinship and female-male dynamics at the site La Almoloya.

(A) Sum of all ROH segments measured with hapROH for CA and BA above 400,000 SNPs in 1240k panel (table S2.21). (B) qpAdm z scores between autosomes and the X chromosome showing no signal for male bias related to steppe ancestry in the model with distal sources (right) or proximal sources (left) (table S2.22). (C) f3-outgroup statistics of the form f3(female, female; Mbuti), f3(female, male; Mbuti), and f3(male, male; Mbuti), highlighting a closer relationship among males than among females or closer relationships whenever they involve at least one male individual (first- and second-degree related pairs are excluded from the calculation) (table S2.23). (D) PWMR values of individuals compared to the PWMR average of the site La Almoloya. Dashed lines indicate a lower PWMR average and, thus, higher pairwise relatedness in adult males than adult females (table S2.24).

To assess a potential sex bias in the steppe-related ancestry contribution as postulated previously (), we applied distal and proximal qpAdm models (text S8) to the X chromosome and autosomes (Fig. 6B and table S2.22). As males contribute, on average, only one-third of the X chromosomes to the next generation (), a lower proportion of any ancestral component on the X chromosome than on the autosomes would thus be indicative of male bias concerning the respective component (, ). In turn, if the ancestral component is statistically higher on the X chromosome, this indicates a female bias. On the basis of this rationale, we do not observe significant male bias in steppe-related ancestry using either distal or proximal sources (Fig. 6B and table S2.22). The fact that the male bias is not detectable could be indicative of an already balanced ancestral component in both sexes, as is reflected in the work of Mittnik et al. (), where the male bias in the steppe component is only detected in Corded Ware, but no longer in Bell Beaker or BA populations. We also took advantage of the recent extensive excavations and sampled all individuals with suitable morphological preservation (n = 86) from the site La Almoloya, resulting in high-quality genome-wide data for 67 individuals that allowed insight into an El Argar community. We estimated the genetic relatedness using a combination of the three methods: pairwise mismatch rate (PWMR) (), READ (Relationship Estimation from Ancient DNA) (), and lcMLkin () for pairs of individuals with more than 1000 shared SNPs (see Materials and Methods, fig. S2, and text S5). All reconstructed pedigrees will be published in an accompanying paper, fully integrated into the anthropological and archaeological context. Thus, we report here the main observations of biological relatedness at the metalevel. We found no first-degree relationships between adult women (0 of 30 adult women analyzed). All first-degree relationships among adults involved at least one adult male (three first-degree relationships among adult men and women). We also found no second-degree relationships between adult women. The few second-degree relationships involving adults were all between males (4 of 18 adult males analyzed).

We then formally tested whether the males of La Almoloya had more close relatives at the site than the females. We calculated f3-output statistics after removing first- and second-degree pairs of the form f3(female, female; Mbuti) and f3(male, male; Mbuti) and found that despite a similar average value in female:female and male:male comparisons (Fig. 6C and table S2.23), males have a more skewed f3-distribution than females, which suggests a more close genetic relationships (e.g., at third to fourth degree) among males with other males (P = 1.9 × 10−18) and males with other females (P = 2.2 × 10−4) than for females with other females (see Materials and Methods). We specifically tested for signals that suggest patrilocality using the PWMR extracted from READ (P0_Normalized) () and built a matrix to calculate the biological relatedness of each adult individual to the entire group (Fig. 6D and table S2.24). We observe a tendency that males are, on average, more closely related to other individuals from the site than females, despite there also being some adult females with relatives at the site, which renders the formal test for patrilocality nonsignificant (see Materials and Methods). At the intrasite level, this result could also suggest a founder effect from the male line with a higher reproductive rate.

DISCUSSION

Our focused study on the southeastern part of Iberia highlights the potential of genetic analysis on a regional level in the light of broader-scale population dynamics in Europe. Overall, we were able to show that the structure of CA populations shows persistent genetic stability and continuity since the Neolithic. This is coupled with a differential HG ancestry in the South in addition to potential early contact between southeastern Iberian CA and Mediterranean populations, which is also revealed by numerous threads of archaeological evidence (–). We provide nuanced evidence for the arrival of steppe-related ancestry in southeastern Iberia by 2200 cal BCE, when widespread social and cultural changes occurred across much of southern Iberia, including the end of the ditched enclosures system of population aggregation and the emergence of BA groups. In the Southeast, all sampled BA groups can be shown to carry steppe-related ancestry, while the more extensively sampled El Argar groups also carry excess Iran_N-like ancestry, which has also been observed in other CA and BA Mediterranean groups. The analysis of one outlier from El Argar reveals a central Mediterranean migrant individual with additional North African ancestry. Overall, we propose that El Argar has likely formed from a mixture of new groups arriving from north-central Iberia, which already carried central European steppe-related ancestry (and the predominant Y-chromosome lineage) and local southeastern Iberian CA groups that differed from other regions in Iberian in that they carried excess Iran_N-like ancestry similar to eastern and/or central Mediterranean groups. Moreover, the large sample sizes of early and late El Argar groups show that the Iran_N-like ancestry contribution in the local CA stratum was not enough to explain this type of ancestry satisfyingly, which, in turn, argues for a continued connection to and genetic influence from the Mediterranean BA at least until the end of the El Argar period.

Last, we were able to shed light onto the population structure of Early Bronze El Argar societies at the intrasite level. In La Almoloya, a closer genetic relationship among males is a strong indicator of patrilineality.

MATERIALS AND METHODS

Experimental design

Archaeological samples

We processed 244 individuals from the Iberian Peninsula and the Balearic Islands dated to a time span from the CA (3000 cal BCE) to the LBA (1200 cal BCE). A detailed description of the archaeological context, sites, and individuals is reported in text S1 and summarized in table S1.1.

Laboratory work

All the laboratory work was performed in the dedicated aDNA facilities of the Department of Archaeogenetics of the Max Planck Institute for the Science of Human History in Jena, Germany.

Sampling

We preferentially sampled petrous bones and teeth (essentially molars). Details of skeletal elements samples are listed in table S1.1. Before sampling, all samples were irradiated with ultraviolet light on either side for 30 min each, after which the surface was cleaned with low-concentration bleach solution (3%). For the petrous bones, we either drilled from the outside, removing the outer layer, or cut the petrous pyramid longitudinally to drill the dense part directly from either side. Molars were cut below the enamel-dentine junction, and tooth powder was drilled from the pulp chamber. We collected between 50 and 100 mg of bone or tooth powder per sample for DNA extractions.

Extraction

DNA was extracted following the modified version of () described in www.protocols.io/view/ancient-dna-extraction-from-skeletal-material-baksicwe.

Library preparation

We generated double-indexed double-stranded () UDG (“ds UDG-half”) () libraries using 25 μl of DNA extract and following established protocols (). Libraries produced per individual are reported in table S1.1.

Shotgun screening and in-solution enrichment (1240k, MT-capture, and Y-capture)

Libraries were sequenced in-house on an Illumina HiSeq 4000 platform to an average depth of 5 million reads and after demultiplexing processed through EAGER (). After an initial quality filter based on the presence of aDNA damage and endogenous DNA higher than 0.1%, libraries were subsequently enriched using in-solution capture probes synthetized by Agilent Technologies for ~1240k SNPs along the nuclear genome and independently for the complete mitogenome () following () and mappable regions of the Y chromosome (), both in-house protocols. The captured libraries were sequenced for 20 to 40 million reads on average using either a single end (1 × 75-bp reads) or paired end configuration (2 × 50-bp reads). Results are reported in tables S1.2 and S1.3.

aDNA data processing

Read processing and DNA damage

After demultiplexing based on a specific pair of indexes, raw sequencing data were processed with EAGER (1.92.59) (). First, the adaptor sequences were clipped with Adapter removal (v2.3.1) () and then mapped against the Human Reference Genome hs37d5 with BWA (v0.7.12) () aln and samse commands (-l 16500, -n0.01, -q30). We removed duplicates (reads with same start and end and same orientation) with DeDup (v0.12.2) (). Last, we determined the deamination rate pattern in our UDG-half libraries using mapDamage (v.2.0.6) (). All our libraries showed expected deamination patterns for UDG-half treatment, and according to this, we trimmed 2 bp on each terminal side before genotyping. A summary of quality statistics is given in table S1.2.

Genetic sex determination

We determined genetic sex at BAM file level by calculating the coverage on the X and Y chromosomes and the autosomes using a bed file of the regions captured by the 1240k SNP array. We normalized the X and Y reads by the autosomal coverage resulting in the X ratio and likewise for the Y ratio (table S1.1 and fig. S1). For uncontaminated libraries, we expect an X ratio of 1 and a Y ratio of 0 for females and an X and Y ratio of 0.5 for males (). Individuals that fall in an intermediate position indicate the presence of DNA contamination from the opposite sex.

Contamination estimation

As first instance, we inspected the scatterplot of X/autosomes ratio versus Y/autosomes ratio from the sex determination to look for outliers that would indicate contamination of the opposite sex. We then estimated exogenous human contamination at nuclear and mitochondrial level. We used method 2 in ANGSD, which estimates heterozygosity at polymorphic sites on the X chromosome in males (), applying a contamination threshold of 3% in individuals with at least 100 X-SNP positions covered twice. In cases where the number of SNPs was lower and for all other libraries from female individuals, we quantified the heterozygosity on the individual mitochondrial reads using ContamMix 1.0.10 () by comparing the consensus mitogenome sequence reconstructed from the MT-capture data with a comparative mitochondrial dataset of 311 global mitogenomes. Only one individual had a rate of contamination of 7.8%, while all others were below 1% (table S1.1). We visually inspected the reads with Geneious (v9.1.8) (), looking specifically for background contamination at all variable sites but found no consistent pattern of contamination that would affect the haplotype calls.

Genotyping

We used our trimmed BAM files to call genotypes with the tool pilleupCaller (https://github.com/stschiff/sequenceTools), which randomly chooses one allele at every SNP position generating pseudo-haploid genotypes. We generated genotype files using the respective SNP lists of the 1240k panel () and the Affymetrix HO panel with the intersecting SNPs (~600,000 SNPs) (, ). Both of them were restricted to only autosomal SNPs for population genetics downstream analysis. Numbers of SNPs covered at least once in both HO and 1240k are given in table S1.1.

Mitochondrial haplogroup assignment

Following targeted enrichment, we could reconstruct complete mitochondrial genomes of 144 individuals, with an average coverage of 1.3 to 964.6 X (tables S1.3 and S1.4). To process mitochondrial DNA data, we extracted reads from mito-capture data and complemented with reads from 1240k capture and shotgun sequencing when necessary using SAMtools v1.3.1 () and mapping to the mitochondrial reference, exclusively. Then, we remapped reads to the revised Cambridge Reference Sequence (rCRS) using Geneious (v9.1.8) (). To determine haplogroups from the consensus sequences reliably, we applied a threshold of at least 3000 reads per individual (table S1.4). Upon visual inspection of the assembly, we called and exported the consensus sequences in FASTA format to be used in HaploGrep2 v2.1.1 [available via https://haplogrep.uibk.ac.at/ for an automated mitochondrial haplogroup assignment based on PhyloTree (mtDNA tree build 17, available via www.phylotree.org/)] (). Using Geneious (v9.1.8), we also visually double-checked whether specific private SNPs were covered and shared in pairs of individuals with an attested degree of relatedness (table S1.4).

Y-haplogroup assignment

By using an in-house protocol for Y capture (), we were able to assign Y-chromosome haplogroups in male individuals. We genotyped the Y-chromosome reads using a Y-SNP list from the ISOGG (International Society of Genetic Genealogy) dataset included in the 1240k and Y-capture probes and by using an in-house script (A.B.R.) (table S1.1). This procedure allowed us to check in a semi-automatic way which positions were covered to assign an ancestral or derived haplogroup or to make corrections to calls in cases where, for instance, a more derived haplogroup was called because of residual ancient damage (C to T or G to A mismatches) in terminal read positions at diagnostic SNPs.

Kinship estimation

We calculated the PWMR (, ) in all pairs of individuals from our pseudo-haploid dataset to double-check for potential duplicate individuals and to determine first- and second-degree relatives. We split our individuals into two groups according to their genetic ancestry profiles (CA and BA individuals), assuming that the median PWMR of the two groups provides a background for unrelated individuals (). We considered one-half of the unrelated value as baseline for identical samples/twins and the midpoint between the two for first-degree relatives (). After running PWMR, we found no duplicates and removed one of each pair of first-degree relatives for downstream population genetic analyses (in all cases, the individuals with lower SNP coverage) (table S1.5 and text S5).

For the kinship study of La Almoloya, we compared and combined kinship estimates from several methods. We used READ () to determine first- and second-degree relatedness among individuals based on the proportion of nonmatching alleles (P0) in nonoverlapping windows of 1 mega base pairs (Mbps) and calculated standard errors. The normalized value of P0 is expected to be between 0 and 0.625 for identical twins, between 0.625 and 0.8125 for first-degree related individuals, between 0.8125 and 0.90625 for second-degree related individuals, and greater than 0.90625 for all other degrees of relatedness (which are reported as “unrelated”) (table S1.6 and text S5).

We also used the method LcMLkin (), which uses genotype likelihoods to estimate the three k-coefficients (k0, k1, or k2) defined by (), which define the probability that two individuals have zero, one, or two alleles identical by descendent (IBD) at a random site in the genome. As the method has to deal with missing data, it uses fractions of sites (haplotypes) that are identical by state (IBS). The length of these analyzed fractions can be modified with the thinning parameter --thin. We performed LcMLkin, exploring different thinning parameters (100,000, 50,000, 10,000, and 5000) to observe the behavior of first-degree estimates and the potential to distinguish between possible parent-offspring or sibling relationships (table S1.7 and text S5). To inform this approach contextually, we combined the genetic data with the archaeological information following this order of relevance: (i) We defined the anthropological age at death to determine sibling relatedness among subadults (before reaching reproductive age) or to inform on the directionality of parent-offspring relatives, as infants cannot be parents. (ii) When two unrelated adult male and female both are first degree related to (an)other individual(s), we identified these as parent-offspring relationship and used these trios to identify potential cases of half siblings. (iii) We used the stratigraphic information to confirm the directionality (Who died first?) to establish the type of first-degree relatedness, e.g., assuming that mothers cannot have passed away long before neonates. (iv) We used uniparentally inherited haplogroups to corroborate the established pedigrees.

Phenotypic and functional variants

We explored a list of 56 SNPs associated with known phenotypic traits. Forty-one of these are linked to physical appearance (skin, eye, or hair color), and we used the HIrisPlex method to calculate the probabilities of phenotype prediction (fig. S8 and table S2.19) (). The other 15 SNPs are associated with the capacity to digest lactose in adulthood (rs4988235 in LCT), alcohol dehydrogenase (rs1229984 in ADH1B), adaptation to metabolize fatty acids (rs174546 in FADS1), predisposition to celiac disease (rs272872 in SLC22A4 and rs653178 in ATXN2/SH2B3), or resistance against infectious diseases [rs4833103 in the gene cluster TLR1-TLR6-TLR10 and rs2269424 in the major histocompatibility complex (MHC)]. We calculated the genotype likelihood to be homozygous or heterozygous for the ancestral (no effect) or derived (effect) allele based on the number of reads from BAM files after removing duplicates. We restricted the analysis to individuals with more than 400,000 SNPs covered in the 1240k panel, and results are summarized in fig. S9 and table S2.20.

Population genetic analysis

Dataset

We merged our final dataset with previously published datasets of ancient and modern individuals reported by the Reich Lab (https://reich.hms.harvard.edu/datasets; please see table S2.1 for a detailed list of references) as well as genotypes from recently published studies (table S2.1) (–).

Grouping and labeling individuals

A detailed description of new labels used in population genetics analysis is described in text S6 and table S2.1. We only included new individuals with more than 40,000 of the 1240k SNPs covered in the downstream population genetics analysis and lowered this threshold to >30,000 SNPs only for individuals forming part of a consistent cluster of individuals. In addition, we excluded groups that have consistently fewer than 60,000 SNPs in qpAdm modeling (table S2.1).

Principal components analysis

We computed PCA using smartpca software from the EIGENSOFT package (v6.0.1) () with the lsqproject and SHRINKMODE option YES and an extended list of modern populations from Eurasia and North Africa () and the Caucasus (). The newly typed ancient individuals were projected onto PC1 and PC2.

F-statistics

F-statistics were computed with ADMIXTOOLS (https://github.com/DReichLab). For F4-statistics, we used the qpDstat and with the activated f4-mode. F3-statistics were calculated using qp3Pop. We used the 1240k panel in all our f-statistics to increase the number of SNPs covered by the ancient individuals. Standard errors were calculated with the default block jacknife.

Admixture modeling

We used ADMIXTURE () to define the main genetic cluster profiles of CA and BA Iberian individuals. We pruned our data for linkage disequilibrium in PLINK () with parameters --indep-pairwise 200 25 0.4 and –maf 0.01, which retained 243,848 autosomal SNPs for the HO dataset. ADMIXTURE was run with default fivefold cross-validation (--cv = 5), varying the number of ancestral populations between K = 2 and K = 18 with different random seeds. Each K was replicated five times.

We used the qpWave and qpAdm programs from the ADMIXTOOLS v3.0 package (https://github.com/DReichLab), with the “allsnps: YES” option to minimally reduce the number of SNPs used and subsequently increase the power to reject models, to model the ancestry in our new reported individuals. With qpWave, we estimated the minimum number of independent contributing ancestral sources of gene flow needed to explain a target population. With qpAdm, we quantified the proportion of genetic ancestry contributed by each source. The ancestry proportions in the target population are inferred on the basis of how the target population is differentially related to a set of reference/outgroups via the source populations.

For all the models applied here, we have used a set of 12 outgroups (Mbuti.DG, Ethiopia_4500BP.SG, Ust_Ishim.DG, Russia_MA1_HG.SG, Italy_Villabruna, Belgium_GoyetQ116_1, Han.DG, Onge.DG, Papuan.DG, AHG, CHG, and Morocco_Iberomaurusian). When we modeled ancestries with proximal sources, we added additional populations or rotated them to the outgroup set. The corresponding rationale is explained in each section of the Supplementary Materials.

Admixture dating

We used DATES v.753 to estimate the date of admixture among different previously fitted ancestral sources (https://github.com/priyamoorjani/DATES). The method works on modeling genotypes of the admixed test population as a linear combination of two population sources calculating the ancestry covariance coefficient. The results are given in generations (28 years, one generation; table S2.9) ().

Sex bias

To test for potential sex bias during periods of admixture, we ran the distal (model 4; described in text S8) and proximal qpAdm models (models 2 and 4b; described in text S8) for the autosomes and in the X chromosome separately using the parameters allsnps: YES and chrom23: YES. If there was a sex bias related to any of the modeled sources of ancestry, we expect to observe increasing (female bias) or decreasing (male bias) amounts of this ancestry on the X chromosome compared to the autosomes, as the paternal contribution on X chromosome is about one-third when compared to the autosomes (the father only contributes one X chromosome to the female offspring). If the ancestry proportion is significantly higher in the autosomes, it implies that the contribution was male driven. Significance was calculated with z scores following previous studies (, ) (table S2.22).

Runs of homozygosity

We used hapROH to calculate the portion of the genome under runs of homozygosity (), which can be applied to low-coverage genotype data. We included all CA and BA individuals from Iberia, which had more than 400,000 SNPs covered in the 1240k panel (table S2.21).

Testing for skewness in f3-statistics

To test for a significant change in the skewness of the distributions of f3-statistics of the form f3(female, female; Mbuti), f3(male, female; Mbuti), and f3(male, male; Mbuti), we measured the sample skewness for each distribution, denoted b, b, and b, respectively. We then looked at the values for the pairwise skewness differences of the form δ = (b − b), where i, j ∈ (mm, mf, ff). We estimated the standard deviation of the δ, denoted by taking 1000 bootstrap samples. We then tested the null hypothesis that the δ = 0 by calculating a z score of the form and then calculated a P value using a standard normal distribution. Our P values were ad hoc adjusted via the Bonferroni method to account for multiple hypothesis testing (table S2.23).

Testing for patrilocality via mean PWMR

To test the overall mean relatedness for each individual, we considered the PWMR for each pair of individuals, denoted p (for individuals i and j). We then calculated the mean PWMR for each individual such that . For the being attributed to groups A and B (male and female), we tested for a significant difference in sample location using a Wilcoxon rank sum test as implemented in the stats library for the R statistical software package (table S2.24).