SCIENCE ADVANCES | RESEARCH ARTICLE ANTHROPOLOGY Copyright © 2021 The Authors, some Genomic transformation and social organization during rights reserved; exclusive licensee the Copper Age–Bronze Age transition in American Association for the Advancement southern Iberia of Science. No claim to original U.S. Government Works. Distributed Vanessa Villalba-Mouco1,2,3*, Camila Oliart4, Cristina Rihuete-Herrada4, Ainash Childebayeva1,3, under a Creative Adam B. Rohrlach1,3,5, María Inés Fregeiro4, Eva Celdrán Beltrán4, Carlos Velasco-Felipe4, Commons Attribution Franziska Aron1, Marie Himmel1, Caecilia Freund1, Kurt W. Alt6,7, Domingo C. Salazar-García8,9,10, NonCommercial Gabriel García Atiénzar11, Ma. Paz de Miguel Ibáñez11, Mauro S. Hernández Pérez11, License 4.0 (CC BY-NC). Virginia Barciela11, Alejandro Romero11,12, Juana Ponce13, Andrés Martínez13, Joaquín Lomba14, Jorge Soler15, Ana Pujante Martínez16, Azucena Avilés Fernández17, María Haber-Uriarte14, Consuelo Roca de Togores Muñoz15, Iñigo Olalde2,18, Carles Lalueza-Fox2, David Reich18,19,20,21, Johannes Krause1,3, Leonardo García Sanjuán22, Vicente Lull4, Rafael Micó4, Roberto Risch1,4*, Wolfgang Haak1,3,23* The emerging Bronze Age (BA) of southeastern Iberia saw marked social changes. Late Copper Age (CA) settle- ments 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 known as the 4.2k event (1–3). Recently, the possibility of substantial During the last centuries of the third millennium BCE, the societies population movements, causing social instability during the third of Europe, the Near East, and Egypt experienced large-scale social millennium BCE, has been proposed as a further explanation for and political upheavals. Settlement abandonment, depopulation, the changes observed at the end of the Copper Age (CA) in central the disappearance of communication networks, and major political and western Europe (4–7). Signs of social and economic turnover are disruptions at the end of the Akkadian empire and the Old Kingdom particularly marked in the southern half of the Iberian Peninsula (8), in Egypt have often been interpreted in the light of a climatic crisis, where the CA is associated with exceptional demographic growth, a diversity of monumental settlements and funerary structures, wide- 1 Department of Archaeogenetics, Max Planck Institute for the Science of Human spread copper metallurgy, and a sophisticated, large-scale manufac- History, Kahlaische Strasse 10, 07745 Jena, Germany. 2Institute of Evolutionary Biology, ture and exchange of symbolic goods, among others [e.g., (9, 10)]. CSIC–Universitat Pompeu Fabra, Barcelona, Spain. 3Department of Archaeoge- Moreover, this period is characterized by a diversity of settlement netics, Max Planck Institute for Evolutionary Anthropology, Deutscher Platz 6, 04103 Leipzig, Germany. 4Department of Prehistory, Universitat Autònoma de types, including fortified sites, ditched enclosures, and so-called Barcelona, Barcelona, Spain. 5ARC Centre of Excellence for Mathematical and Statis- megasites, some of which exceeded 100 ha in size (e.g., Valencina de la tical Frontiers, School of Mathematical Sciences, The University of Adelaide, Adelaide Concepción and Marroquíes Bajos) and all of which were formed at SA-5005, Australia. 6Center of Natural and Cultural Human History, Danube Private University, Steiner Landstr. 124, A-3500 Krems, Austria. 7Department of Biomedical around 3300 to 2800 BCE, therefore predating the Bell Beaker horizon. Engineering, University of Basel, Gewerbestrasse 14-16, CH-4123 Allschwil, Switzerland. This period is also associated with a major increase in inter- 8 Grupo de investigación en Prehistoria IT-1223-19 (UPV-EHU)/IKERBASQUE—Basque connectedness and mobility. On the basis of available radiogenic Foundation for Science, Vitoria, Spain. 9Departament de Prehistòria, Arqueologia i Història Antiga, Universitat de València, València, Spain. 10Department of Geological (Sr) isotope studies, the percentage of southern Iberian individuals Sciences, University of Cape Town, Cape Town, South Africa. 11Institute for Research who were buried in locations other than where they grew up ranges in Archaeology and Historical Heritage (INAPH), Universidad de Alicante, 03690 between 8 and 74% (11, 12). Ivory from Africa and the Near East Alicante, Spain. 12Departamento de Biotecnología, Universidad de Alicante, 03690 (13–15), amber from Sicily (16), and ostrich eggshells from Africa Alicante, Spain. 13Museo Arqueológico Municipal de Lorca, Murcia, Spain. 14Depart- ment of Prehistory, Universidad de Murcia, Murcia, Spain. 15MARQ, Alicante, Spain. (17) are indicative of transregional connections. However, evidence 16 Arqueología Estudios, Murcia, Spain. 17Arqueología y Diseño Web S.L. (Grupo of a strong political centralization and economic inequality remains Entorno), Floridablanca 14, 1.°D, 30800 Lorca, Murcia, Spain. 18Department of Ge- elusive or inconclusive (18–21). netics, Harvard Medical School, Boston, MA 02115, USA. 19Department of Human Evolutionary Biology, Harvard University, Cambridge, MA 02138, USA. 20Broad In- Archaeogenetics has suggested that the remarkable development stitute of Harvard and MIT, Cambridge, MA 02142, USA. 21Howard Hughes Medical during the (south) Iberian CA was coupled with a strong population Institute, Harvard Medical School, Boston, MA 02115, USA. 22Department of Prehistory continuity attested since the Neolithic [e.g., (4, 6, 7, 22–27)]. How- and Archaeology, University of Seville, 41004 Sevilla, Spain. 23School of Biological ever, the Late CA anthropological and archaeological records from Sciences, The University of Adelaide, Adelaide SA-5005, Australia. *Corresponding author. Email: villalba@shh.mpg.de (V.V.-M.); robert.risch@uab.cat the north and central Iberia show the first individuals carrying (R.R.); haak@shh.mpg.de (W.H.) “steppe-related ancestry” by ~2400 calibrated (cal) BCE, which are Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 1 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE often but not exclusively linked to Bell Beaker–associated artifacts in the intramural burials, monumental architecture (31, 36), and the (6, 7). In parallel, African ancestry was also reported in one individ- spatial relation between both (31, 36, 37). ual, which suggests discrete movement/mobility of people (7). Lull and colleagues (30, 38) have argued that social conflict and The beginning of the Bronze Age [Early Bronze Age (EBA)] in environmental degradation resulting from deforestation and exten- Iberia (2200 to 1550 cal BCE) marks a clear population turnover, sive dry farming could have led to the abandonment or destruction suggested by both the omnipresence of steppe-related ancestry in all of all El Argar settlements around 1550 BCE. However, the rise of a individuals directly postdating 2200 BCE. An even more notable shift similar economic organization, architecture, and funerary record in can be observed in the frequency of Y-chromosome haplogroups in inner Alicante, as highlighted by Cabezo Redondo, suggests that males, which are almost exclusively of the R1b-P312 type that was during the peak of El Argar at least some groups managed to estab- completely absent in Iberia before 2400 BCE (6, 7, 25, 26). lish themselves in a territory that was influenced by the Valencian The turn from Late CA to the EBA at the end of the third millennium BA, a neighboring “cultural group” located in southeastern Iberia. BCE saw the demise of fortified settlements such as Los Millares and In this study, we aimed to understand the importance of popula- ditched megasites such as Valencina and Perdigões in southern Iberia, tion dynamics in the collapse of the highly dynamic Iberian CA while in the southeast it is concomitant with the rise of new hilltop world, the rise and development of El Argar, as well as its relations occupations of smaller sizes (<0.5 ha) (8, 28, 29). Substantial and between neighboring BA groups in western Europe and the Balearic densely built hilltop settlements, distinguished by a specific intra- Late BA (LBA). We also explore the genetic makeup of the BA groups mural burial rite and characteristic ceramic and metal types, appear in Iberia and the Balearic Islands in relation to the other western around 2200 BCE in the fertile tertiary basins framed by mountain ranges, and central Mediterranean islands, such as Sardinia and Sicily. In running parallel to the coast of southeastern Iberia (8). This area of circa total, we characterize the genomic profiles of 96 BA individuals (c.) 3500 km2 is considered to be the core of the El Argar “culture,” associated with El Argar and contemporaneous societies, as well as which is one of the most outstanding examples of an early complex 34 CA and 6 LBA individuals. society in European prehistory with evidence for social stratification (8, 30, 31). The origins of El Argar are still unclear, as there are no hybrid contexts where El Argar elements appear in Late CA settlements RESULTS or vice versa. Although the early El Argar material record shares Here, we report genome-wide data in the form of 1.24 million some traits with the Bell Beaker complex (30), such as V-perforated ancestry-informative single-nucleotide polymorphisms (SNPs) for buttons, Palmela-type points, or perforated stone plaques so-called 136 southern Iberian individuals, covering a time span of 2000 years “archers’ wristguards,” the characteristic Bell Beaker pottery is absent. from the Late Neolithic (LN)/CA (3300 cal BCE) to the LBA Upon the discovery of the monumental fortification in the 5-ha hilltop (1200/1000 cal BCE) (Fig. 1, table S1.1, and text S1). Five of these settlement of La Bastida, dated to around 2200 BCE, a possible eastern individuals were also shotgun-sequenced to 1.5× to 5.2× coverage Mediterranean contribution was reconsidered (32). The intramural to reconstruct the complete genomes to add to the growing record burials in large storage vessels (pithoi), the circulation of silver rings of ancient western Eurasian genomes for future studies. Our new and bracelets, and the characteristic footed Argar cup have also dataset includes various groups and types of sites such as sepulchral been interpreted as signs of Aegean or Near Eastern contacts (33), caves from the LN/CA (Cova d’en Pardo, n = 7; Cueva de las Lechuzas, although all these features emerged during later phases of El Argar. n = 10), simple pit burials of pre–Bell Beaker CA megasites [Valencina The genealogy of different characteristic material traits of early (PP4-Montelirio), n = 11], a collective hypogeum with Late CA El Argar, such as apsidal buildings, intramural burial, metal casting burials (Camino del Molino, n = 6), and EBA pit graves (Molinos de technology, and the halberds as distinguished weapons, is reminiscent Papel, n = 3) to provide a diverse comparative dataset for the Ar- of several social developments in southeastern, central, and western garic societies, which are at the center of our study. For the latter, Europe, with possible links of still uncertain origin (31). we include the almost complete and archaeologically well-defined Between 2000 and 1800/1750 BCE, El Argar expanded across a sites of La Almoloya (n = 67) and La Bastida (n = 10), alongside wider region of southeastern Iberia and entered the central Spanish Cerro del Morrón (n = 3) and Lorca town [Madres Mercedarias Meseta. Characteristic El Argar items, such as halberds, are also Church (n = 1), Los Tintes (n = 1), and Zapatería (n = 1) sectors]. present beyond this territory. The leading figures within El Argar We also analyzed individuals of the so-called Valencian BA, including society seem to have been a class of warriors armed with halberds Cabezo Redondo (n = 1), Peñón de la Zorra (n = 1), Puntal de los and daggers. These weapons, occasionally associated with golden Carniceros (n = 4), and La Horna (n = 3), individuals of the BA from arm rings, also became the insignia of political domination in the Catalonia (Miquel Vives, n = 1), as well as LBA individuals from the male elite burials of central Europe around 2000 to 1800 BCE, a time Balearic Islands of Menorca (Es Forat de ses Aritges, n = 6) (Fig. 1, when this specialized weapon had already gone out of use in the rest table S1.1, and text S1). of Europe (34). At the same time, a growing part of the population, We initially screened 244 individuals for DNA preservation by particularly children, was buried in cists, artificial caves, pottery shotgun sequencing ~5 million Illumina single end reads per partial vessels, or pits inside the distinctive settlements of El Argar. uracil-DNA glycosylase (UDG)–treated DNA library (39), followed During its final phase (1800/1750 to 1550 BCE), El Argar reached by assessment of % endogenous human DNA (>0.1%), average read an outstanding level of economic and social development. A vast length between 40 and 75 base pairs (bp), and characteristic ancient amount of grinding tools, large workshops, and storage facilities DNA (aDNA) damage patterns >3% (tables S1.1 and S1.2). Libraries found in the larger hilltop settlements (1 to 6 ha) imply that certain fulfilling these criteria were retained for hybridization capture of groups managed to control the flow of resources and the workforce of 1.24 million informative sites (1240k SNP panel) in the human ge- wider regions (35). The establishment of a dominant and hereditary nome (22). Libraries were independently captured for mitochondrial class, as well as increasing social asymmetries, become recognizable genomes (tables S1.3 and S1.4 and text S2) [(40) modified following Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 2 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE 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). Direct- ly 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. (41)] and mappable regions of the Y chromosome (table S1.1 and one XXX (trisomy X) individual (table S1.1, text S4, and fig. S1). We text S3) (42). Following standard aDNA processing pipelines (43), obtained endogenous human aDNA contents ranging from 1.3 to we quantified contamination rates at the autosomal (in males) (44) 29.5% on the targeted 1240k SNPs, which correspond to 2004 to and mitochondrial levels (both genetic sexes) (45), showing low 831,086 SNP positions covered (tables S1.1 and S1.2). For popula- contamination estimates for the individuals in this study (<3% for tion genetic analyses, we only retained individuals with more than nuclear and <1% for mitochondrial DNA; see Materials and Methods). 40,000 SNPs on the 1240k panel, excluding 18 low-coverage indi- By determining the genetic sex (46), we observed two cases of sex viduals in downstream analyses. We performed kinship analysis chromosome aneuploidies, one XXY (Klinefelter) individual and (47–49) for all newly typed individuals (fig. S2; tables S1.5, S1.6, and Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 3 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE S1.7; and text S5), the result of which will be featured in a separate Iran_N or Jordan_PPNB as a fourth source (Anatolia, WHG, GoyetQ2, study. From each pair of first-degree relatives, we then excluded the and Iran_N/Jordan_PPNB) further improves the model fit for individual with lower SNP coverage. In summary, for in-depth SE_Iberia_CA by an order of magnitude, i.e., from P = 0.008 to population genetic analysis, we merged genome-wide data from P = 0.014 (with Iran_N) or P = 0.027 (with Jordan_PPNB) (table S2.4, 122 newly typed individuals with a comparative set of previously text S7, and fig. S3). In turn, adding other populations as a fourth published ancient and modern-day individuals, which are reported source does not improve the model fit (text S7). The unknown source by the Reich Lab (https://reich.hms.harvard.edu/datasets; please see contributing to SE_Iberia_CA groups is likely to have carried an the Supplementary Materials for a detailed list of references), as well excess of Levantine and/or Iran_N-like ancestries compared to the as genotypes from recently published studies (table S2.1) (50–53). distal source Anatolia_Neolithic, as these together have been found admixed in Anatolian and Levantine CA groups (~6000 to 5000 BCE) Genetic substructure in the Iberian CA (53). This finding suggests a subtle contribution that was spread early We first explored the genetic affinities of the newly typed CA indi- along the Mediterranean or, alternatively, different sources of early viduals by performing a principal components analysis (PCA) with farmer ancestry during the Neolithic transition with varying pro- relevant ancient individuals projected onto PCs calculated from a portions of Levantine and/or Iran_N-like components when com- set of modern-day West Eurasian populations genotyped with the pared to Anatolia_N used here. Removing Anatolian HG (AHG) from Human Origins (HO) SNP panel (Fig. 2A) (54). The new CA indi- the outgroups also improves the model fit (P = 0.045; table S2.4), viduals from southern Iberia fall onto a position that partially overlaps indicating that the Neolithic ancestry is not well represented by using with previous Middle Neolithic (MN), Middle/Late Neolithic (MLN), Anatolia_N as a distal proxy and might come from another farmer and CA (non-steppe) groups from Iberia but are slightly positively group more similar to AHG than Anatolia_N (text S7). More indi- shifted in their coordinates for PC1 toward previously published viduals from the Neolithic and CA across the Mediterranean will be Early Neolithic (EN) groups from Iberia and later groups such as needed to track this contribution more confidently. Sardinia Chalcolithic, suggesting an equally small hunter-gatherer The Magdalenian-related ancestry was not detectable in other (HG) ancestry contribution in the CA individuals of southern Iberia. contemporaneous Mediterranean populations (e.g., Sardinia_Neolithic, These results run counter to previous suggestions of steppe-related Sardinia_Chalcolithic, Sardinia_EBA, Sicily_EBA, and Italy_CA), ancestry diffusion into southern Iberia in the early CA based on which renders gene flow in the direction from southern Iberia to the supposed formal similarities in one element of the Valencina mate- central Mediterranean less plausible. However, applying the same rial culture (55). To formally test the differential contribution of the qpAdm model (text S7, fig. S3, and table S2.5), we detect a previously main HG source in European Neolithic farmers, called Western HG unreported amount of Iran_N-like ancestry in central Mediterranean (WHG), which represents the predominant form of ancestry in pre- groups from Sardinia, ranging from 2.8 ± 1.2% in Sardinia_Chalcolithic farming western Europe (4, 56), we calculated f4-statistics of the form to 5.8 ± 1% in Sardinia_Nuragic_BA. Adding Iran_N as a source to f4(SE/SW_Iberia_CA, N/NE/C_Iberia_CA; WHG, Mbuti) and ob- model Sicily_EBA improves the model fit but without reaching tained a significantly negative value for N_Iberia_CA and C_Iberia_CA P values ≥0.05, making gene flow from the Italian Peninsula to (Fig. 2B and table S2.2), indicative of higher WHG ancestry in these Sardinia more likely (Italy_CA also shows Iran_N-like ancestry) than regions, a signal that is also present in the preceding MLN period, from Sicily (Sicily_EBA) (text S7, fig. S3, and table S2.5). Notably, when as shown by f4-statistics of the form f4(SE/SW_Iberia_MLN, modeling Sicily_EBA, we obtain P values ≥0.05 by removing AHG N_Iberia_MLN; WHG, Mbuti) (table S2.2). from the outgroups in a three-source model (Anatolia_Neolithic, WHG, Leveraging insights from previous studies about the differential and Yamnaya_Samara) or by removing Morocco_Iberomaurusian HG legacy in western Europe (7, 57–59), we explored the HG ancestry in a four-source model (Anatolia_Neolithic, WHG, Yamnaya_Samara, of geographically diverse CA Iberians in different ways. We first and Iran_N), which points to genetic substructure in the Mediterranean calculated f4-statistics of the form f4(GoyetQ2, WHG; test, Mbuti), before the arrival of steppe-related ancestry. However, this finding where test represents all MLN and CA groups (Fig. 2C and table S2.3). does not rule out limited genetic input from other Mediterranean Despite the resulting f4-values being significantly negative for all test populations to southeastern CA Iberians (text S7). groups (indicating WHG admixture), we observed a geographically related gradient with southern groups (MLN and CA) showing lower Genetic turnover in the southern Iberian BA and negative values, indicating a differential relationship with GoyetQ2, the rise of El Argar who acts as a proxy for Magdalenian-related ancestry (Fig. 2C). We When compared to the preceding CA groups, individuals associ- confirmed this subtle signal using qpAdm outgroup-based ancestry ated with El Argar BA and SE_Iberia_BA form a dense cloud in PC modeling with Anatolia_Neolithic, WHG (Loschbour and KO1), space that is shifted toward populations with steppe-related ances- EHG [Eastern European HGs; the predominant form of ancestry try from central Europe and positioned between individuals from in pre-agropastoralist eastern Europe; (4)], and GoyetQ2 as distal the Iberian CA and those of the north Iberian BA (Fig. 3A). A sim- sources (table S2.4). Here, we found that the CA groups from southern ilar genetic shift is noticeable in the ADMIXTURE (60) results (text S8 Iberia differ with respect to the overall quantity of HG ancestry, but and fig. S4), where Iberian BA individuals show an additional com- the latter is also qualitatively different. Estimates of Magdalenian-­ ponent when compared to preceding CA individuals, and this is fur- associated HG ancestry in southern Iberian CA individuals vary ther supported by an f4-test of the form f4 (Anatolia_Neolithic, test; from 6.1 ± 1.3% to 7.2 ± 1.3%, which reflects the geographic struc- Yamnaya_Samara, Mbuti), in which test iterates through newly re- ture of local HG groups (7, 57) and suggests a certain amount of ported CA and BA individuals (Fig. 3B and table S2.6). Here, the genetic continuity since the Neolithic (Fig. 2D and table S2.4). In the shift to negative f4-values in EBA individuals indicates shared case of SE_Iberia_CA, adding GoyetQ2 as a third source improves drift with Yamnaya_Samara individuals, which is absent in earlier the model fit slightly from P = 2.53 × 10−6 to P = 0.008. Adding CA individuals. Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 4 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE A New ancients SE Iberia CA (Cueva de las Lechuzas) SE Iberia CA (Cueva d'En Pardo) SW Iberia CA (Valencina) 0.006 SE Iberia CA (Camino del Molino) 0.005 0.004 0.003 0.002 0.001 PC2 0 0.001 0.002 Published ancients 0.003 CHG WHG N Iberia CA EHG NW Iberia CA Stp Iran N NE Iberia CA 0.004 Yamnaya Samara C Iberia CA Iberia EN C Iberia CA Stp Iberia MN SW Iberia CA 0.005 N Iberia MLN SE Iberia CA NE Iberia MLN Sardinia Chalcolithic SE Iberia MLN Italy CA SW Iberia MLN Italy Remedello Late CA 0.006 0.010 0.005 0 0.005 PC1 B C D GoyetQ2 NE Iberia CA SW Iberia CA AnatoliaN WHG GoyetQ-2 C Iberia CA SE Iberia CA 0.004 C Iberia CA (0.115) N Iberia CA SE Iberia MLN NE Iberia CA (0.107) 0.002 0 SW Iberia MLN NE Iberia MLN (0.948) f4(SW Iberia CA, test; WHG, Mbuti) NE Iberia CA 0.005 N Iberia CA (0.108) C Iberia CA N Iberia MLN (0.422) NE Iberia CA NE Iberia MLN SE Iberia CA (0.009) N Iberia CA 0.006 SW Iberia CA (0.834) C Iberia CA SE Iberia MLN (0.567) N Iberia CA N Iberia MLN SW Iberia MLN (0.120) 0.002 0 0.005 0 WHG f4(SE Iberia CA, test; WHG, Mbuti) f4(GoyetQ2, WHG; test, Mbuti) 0 1 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). This observation suggests a substantial amount of steppe-related Significantly positive f4-values confirmed the presence of steppe-­ ancestry in El Argar BA individuals, which we tested formally and di- related ancestry in all BA individuals. We then tested for differ- rectly with f4-statistics of the form f4(Argar_Iberia_BA/SE_Iberia_BA, ences in affinity to steppe-related ancestry by contrasting northern SE_Iberia_CA; Yamnaya_Samara, Mbuti) (fig. S5A and table S2.7). versus southern BA individuals using f4(N/NE/C_Iberia_BA, Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 5 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE A Published ancients 0.004 Iberia CA (published) Iberia CA Stp (published) N Iberia BA 0.003 NE Iberia BA C Iberia BA SW Iberia BA SE Iberia BA Valencian 0.002 France Bell Beaker Germany Bell Beaker Italy Bell Beaker 0.001 Hungary Bell Beaker Sardinia EBA Greece Minoan Greece Mycenaean 0 Sardinia Nuragic BA Sicily BA Sicily EBA 0.001 Sicily LBA PC2 Sicily MBA Iran Seh Gabi Chalcolithic Alalakh MLBA 0.002 Ebla EMBA ZAP002 Anatolia Chalcolithic Anatolia EBA 0.003 Caucasus Kura Araxes Russia Yamnaya Samara Greece EBA Greece MBA 0.004 New ancients SE Iberia CA (C. Lechuzas and C. Pardo) La Bastida (SE Iberia BA Argar) 0.005 SE Iberia CA (Camino del Molino) Madre Mercedarias (SE Iberia BA Argar) Miquel Vives (NE Iberia BA) Zapateria (SE Iberia BA Argar) Forat des ses Aritges (Menorca LBA) 0.006 Molinos de Papel (SE Iberia EBA) Cabezo Redondo (SE Iberia BA Valencian) Los Tintes (SE Iberia BA Argar) La Almoloya (SE Iberia BA Argar Early) Puntal Carniceros (SE Iberia BA Valencian) 0.007 La Almoloya (SE Iberia BA Argar Late) La Horna (SE Iberia BA Valencian) 0.004 0.002 0 0.002 0.004 PC1 B f4 (Anatolia Neolithic, test; Yamnaya Samara, Mbuti) 0.002 H2 E1b Male I2a R1b Female 0.001 G2a 0 0.001 0.002 0.003 0.004 CDM001 CDM002 CDM003 CDM004 CDM005 CDM006 CDP001 CDP002 CDP003 CDP006 CDP008 CDP009 CDP011 CLL001 CLL002 CLL003 CLL004 CLL005 CLL006 CLL008 CLL009 CLL010 CLL011 MON013 MON014 MON015 MON016 MON017 MON019 MON020 MON021 MON029 MON033 MON036 ALM001 ALM002 ALM003 ALM004 ALM006 ALM007 ALM008 ALM014 ALM015 ALM016 ALM017 ALM018 ALM019 ALM020 ALM021 ALM024 ALM025 ALM026 ALM027 ALM028 ALM029 ALM030 ALM031 ALM032 ALM034 ALM035 ALM036 ALM038 ALM039 ALM040 ALM041 ALM042 ALM043 ALM044 ALM046 ALM047 ALM048 ALM049 ALM050 ALM051 ALM052 ALM053 ALM055 ALM056 ALM057 ALM058 ALM060 ALM062 ALM063 ALM064 ALM067 ALM068 ALM069 ALM070 ALM071 ALM073 ALM075 ALM076 ALM077 ALM078 ALM079 ALM080 ALM081 ALM084 ALM086 ALM087 ALM088 BAS002 BAS003 BAS017 BAS018 BAS022 BAS023 BAS024 BAS025 BAS026 BAS027 CBR004 CMO001 CMO002 CMO003 EFA006 EFA007 EFA008 EFA009 EFA010 EFA011 LHO001 LHO002 LHO003 LOT001 MDP001 MDP002 MDP003 MIV001 MMI004 PLZ001 PUC001 PUC002 PUC003 PUC004 ZAP002 Copper Age Bronze Age Fig. 3. Evidence for steppe-related ancestry in EBA individuals from southeastern Iberia. (A) West Eurasian PCA calculated with modern populations (gray dots) (23) 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). Argar_Iberia_BA/SE_Iberia_BA; Yamnaya_Samara, Mbuti) (fig. S5B f4-statistics (Fig. 3A and fig. S5), which suggests that the contribu- and table S2.8). The resulting f4-values confirmed a smaller amount tion is homogenized across the population. of steppe-related ancestry in individuals from the Argaric sites La The newly analyzed LBA individuals from the Balearic Islands Almoloya and La Bastida compared to the rest of Iberia_BA groups, (Aritges_LBA) showed less steppe-related ancestry than the published especially when compared to those from northern Iberia (fig. S5B and Late CA individual from Mallorca (previously named Mallorca_EBA), table S2.8), despite the complete turnover to lineage R1b-P312 (except confirming a decrease of steppe-related ancestry over time in for one subadult male in La Bastida) visible in the Y-chromosome agreement with (51). The amount of steppe-related ancestry in record (Fig. 3B, table S2.6, and text S8). However, at the intrasite Mallorca_CA_Stp is similar to the first CA groups in Iberia level, we observe no significant differences with respect to the (C_Iberia_CA_Stp and NW_Iberia_CA_Stp) (fig. S5A) and could amount of steppe-related ancestry between the early and late have arrived at the same time, as there is no clear evidence of human phase of La Almoloya and La Bastida based on PCA and formal occupation in the archipelago before their arrival (61, 62). Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 6 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE Olalde and colleagues (7) had proposed for the Iberian Peninsula ~2100 BCE, coincides with the foundation of monumental settle- that the first contribution of steppe-related ancestry was diluted ments, such as Las Motillas fortifications. during the BA by admixture with descendants of local CA groups, This shows that while the genetic contribution of steppe-related but increased again during the LBA–Iron Age during a second pulse. ancestry to Iberia was a long-term process starting around 2400 cal We do not detect this second pulse in the Balearic Islands, in agree- BCE in the northern and central regions (7), from where it spread ment with the findings of (51), albeit based on one LBA individual southward over ~300 years. At a local scale, this change might have from Menorca (fig. S5, A and B). occurred faster. A similar situation might have existed in central The situation observed in southern Iberia offers a specific example Portugal, where we still find individuals with no steppe-related of the demographic changes occurring during the transition from ancestry in collective CA burials (Galería da Cisterna and Cova da CA to EBA societies around 2200 BCE. The absence of steppe-related Moura) around 2300 to 2200 cal BCE. However, after 2100 cal BCE, ancestry in our newly typed individuals from the Early CA at all individuals from all sites carry steppe-related ancestry, in line Valencina (c. 2900 to 2800 BCE) and the Late CA collective burial with R1b-P312 becoming the predominant Y-chromosomal lineage cave Camino del Molino, which is contemporaneous with the Beaker present not only in El Argar but also in the rest of BA Iberia. Horizon (albeit without Beaker pottery) and the first appearance of steppe-related ancestry in a double tomb with Bell Beaker–associated Mediterranean and central European ancestries shaped artifacts at Molinos de Papel, a site only 400 m away from the latter, the genetic profile of southeastern BA groups in Iberia provides the temporal framework for the arrival of steppe-related To explore the genetic turnover and the contribution of the local ancestry in southern Iberia. Taking the youngest date from Camino groups to the newly formed BA genetic profile in Iberia, we system- del Molino [Beta-261524, 3830 ± 40 years before present (BP)] and atically tested a series of qpAdm models. We started by using the the oldest from Molinos de Papel (MAMS-11826, 3780 ± 30 BP) as distal ancestry sources Anatolia_N, WHG, GoyetQ2, Yamnaya_Samara, reference points, we can show that both C14 dates overlap at the and Iran_N to model the genetic ancestry components of Iberian 95% level and date the arrival of steppe-related ancestry in southeastern BA groups (table S2.10 and fig. S6). We found that the local traces Iberia to ~2200 cal BCE. At a cross-regional scale, none of the south of GoyetQ2, a characteristic but variable component of southern Iberian CA individuals analyzed so far show evidence of any steppe-­ Iberia CA individuals, were no longer detectable, suggesting a dis- related ancestry, whereas the oldest individuals from the El Argar sites solution of geographic substructure in BA Iberia with respect to HG of La Bastida (BAS024 and BAS025) and La Almoloya (ALM019), ancestry. We explain this by the spread of steppe-related ancestry directly dated to the 21st century BCE based on the radiocarbon evi- from North to South (7) that also contributed northern and central dence, show clear evidence of such ancestry. Consequently, the Iberian ancestry to the South, diluting the subtle GoyetQ2 signal to population mixture must have occurred in the preceding 150 years a level beyond the limits of detection (text S8). By using the same at the latest. Estimating the time of admixture of the putative sources qpAdm model, we also observed that Almoloya_Argar_Early, using the program DATES (Distribution of Ancestry Tracts of Evo- Almoloya_Argar_Late, SE_CabezoRedondo_BA, and Bastida_Argar lutionary Signals; https://github.com/priyamoorjani/DATES) pro- cannot be modeled with Yamnaya_Samara as a single source but find duced compatible admixture dates in all new BA groups, ranging better support with a combination of Iran_N and Yamnaya_Samara, from ~2550 BCE ± 189 years to ~1642 BCE ± 133 years (table S2.9 however, without reaching P values ≥0.05 in Almoloya_Argar_Early and text S8). This could have occurred after the collapse of the CA and Late and SE_CabezoRedondo_BA (table S2.10 and fig. S6). These social system in some aggregation megasites such as Valencina at El Argar groups (Almoloya and Bastida) are also slightly shifted to around 2350 BCE (63) and more or less contemporary in others, the right on the PC1 axis, in the direction of Mediterranean BA such as Los Millares at around 2300/2200 BCE (64). The demise of groups with excess Iran_N-like ancestry, such as “Minoans,” who CA cultures ~2350 to 2200 cal BCE has been linked to the so-called only carry Iran_N-like ancestry but not steppe-related ancestry, or 4.2–thousand year cal BP climatic event, as there is a temporal overlap “Mycenaeans,” who carry a mix of both (71), and that has also been between this event and the collapse of large megasites (63, 65–68), shown for some BA individuals from Sicily_MBA (51) and for although direct, conclusive, paleoclimatic evidence for Iberia remains Sardinians here (Fig. 3A). sparse (8, 69). The social and economic changes that lead into the To explore the reasons for the model rejection observed in BA, as well as the arrival of the new genetic ancestry, appear to have southeastern BA groups, we tested several qpAdm models (text S8). spread from North to South Iberia (7, 70) and thus might be linked. We used two- and three-way competitive qpAdm models to test However, it remains unclear whether this expansion was opportu- whether a third source with Iran_N-like ancestry was needed. In the nistic (i.e., followed the consequences of a potential climatic deteri- two-way model (table S2.11), Germany_Bell_Beaker (the largest oration) or causal to the actual shifts observed in archaeology. number of individuals representing BB) was used as a fixed proximal Notably, the arrival of new genetic ancestry is not paralleled by the ancestry source together with a local CA source (either N_Iberia_CA, same social changes in all regions of Iberia. In the coastal regions C_Iberia_CA, or SE_Iberia_CA). For target groups in which these of southeastern Iberia, the shift to single or double burial practices two sources were rejected, we iteratively tested central and eastern inside settlements occurred ~2200 BCE, coinciding with the be- Mediterranean populations as a third source, as some of them are ginning of El Argar. However, in inner Iberia, the associated grave known to carry Iran_N-like ancestry (Fig. 4A and table S2.11). Using goods and settlement pottery are more related to late Bell Beaker CA C_Iberia_CA as a local source of ancestry and Germany_Bell_Beaker as (leaf-shaped tanged point, V-perforated button, and incised Beaker a proxy for steppe ancestry, only models with Almoloya_Argar_Early pottery) than to the material culture of the early El Argar “core” ter- and Late as a target were rejected, suggesting that a third component ritory (8). In La Mancha, the southern part of the central Spanish was required for these Iberian BA groups. We found that adding plateau, the first individual with a confirmed steppe-ancestry com- any Mediterranean population as a third source improves the ponent (7), a male from tomb 4 of Castillejo del Bonete dated to model fit, albeit without reaching P values >0.05. However, we obtained Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 7 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE A Iran CA C Iberia CA Germany BB Iran CA SE Iberia CA Germany BB 1.0 1.0 0.8 0.8 0.6 0.6 Not feasible Not feasible Not feasible Not feasible 0.4 0.4 0.2 0.2 0.319 0.179 0.720 0.128 0.053 0.089 0.380 0.108 0.657 0.755 0.540 0.303 0.872 0.742 0.870 0.539 0.169 0.284 0.616 0.230 0.374 0.002 0.742 0.605 0.187 0.445 0.197 0.148 0.092 0.140 0.548 0.182 0.365 0.854 0 0 1.0 1.0 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0.319 0.179 0.720 0.008 3.6E-06 0.089 0.380 0.026 0.657 0.755 0.540 0.303 0.872 0.742 0.870 0.539 0.165 0.284 0.616 0.020 0.230 0.374 0.004 0.005 0.742 0.605 0.187 0.445 0.001 0.197 0.040 0.148 0.092 0.140 0.548 0.182 0.365 0.854 0 0 C NW Mo Alm Alm Ba SE SE Ce SE SW N I NE C C Ma Fo Me Ar C NW Mo Alm Alm Ba SE SE Ce SE SW N I NE C C Ma Fo Me Ar Ibe li s r b Ib Ib ll r n it Ibe li s r b Ib Ib ll r n it ria Iber nos oloy oloy tida Iberia Cab ro V Iberia Iber eria Iberia eria eria orca men orca ges L ria Iber nos oloy oloy tida Iberia Cab ro V Iberia Iber eria Iberia eria eria orca men orca ges L CA ia C Pap a A aA A ez irg ia B ter CA ia C Pap a A aA A ez irg ia ter rga rgar BA A oRe en L BA V BA A BA BA BA C L Co A St a MB BA BA rga rgar BA A oRe en L BA V BA A B BA BA BA C L Co A St a MB BA BA Stp A S el E rga d Stp A S el E rga d tp BA r r rga on a a gI p A tp BA r r rga on a a gI p A Ea L rly ate r do te Ar lencia Ea L rly ate r do te Ar lencia BA gar n BA gar n B Iran N C Iberia CA C Iberia CA Stp Iran N SE Iberia CA C Iberia CA Stp 1.0 1.0 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0.412 0.908 0.119 0.092 0.900 0.351 0.133 0.810 0.847 0.346 0.283 0.092 0.798 0.522 0.589 0.076 0.174 0.396 0.412 0.725 0.014 0.002 0.634 0.663 0.139 0.783 0.188 0.221 0.369 0.092 0.442 0.192 0.588 0.076 0.174 0.106 0 0 1.0 1.0 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0.412 0.908 5.1E-04 2.5E-08 0.046 0.351 0.036 0.810 0.847 0.346 0.283 0.092 0.798 0.522 0.588 0.076 0.174 0.045 0.412 0.725 0.020 2.4E-04 0.634 0.663 0.139 0.783 0.189 0.221 0.369 0.092 0.442 0.192 0.588 0.076 0.174 0.106 0 0 NW Mo Alm Alm Ba SE S Ce S S N N CI CI Ma Fo Me A NW Mo Alm Alm Ba SE S Ce S S N N CI CI Ma Fo Me A Ibe lino o o st Ib E C rro E Ib W Ib Ibe E Ib be b ll ria eria orca men r no ritge Ibe lino o o st Ib E C rro E Ib W Ib Ibe E Ib be b ll ria eria orca men r no ritge ria sP loya loya ida A eria abe V eria eria ria e BA ria B ter rca L s LB ria sP loya loya ida A eria abe V eria eria ria e BA ria B ter rca L s LB CA ape Arg Arg rga BA zoR irgen B B BA BA C Co A St a MB BA A CA ape Arg Arg rga BA zoR irgen B B BA BA C Co A St a MB BA A Stp l EB ar E ar L r Arg edo La A Va A A g p Stp l EB ar E ar L r Arg edo La A Va A A g p A arly ate ar nd t e le I A A arly ate ar nd t e le I A o B Arg ncia o B Arg ncia A ar n A ar n 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. P values >0.05 when using Iran_Chalcolithic as a third source, which contributing to all Iberian groups, and iterated through the proxi- suggests that a higher proportion of Iran_N-like ancestry is needed mal local CA sources, which in these models contribute 100% of (text S8 and table S2.11). Notably, when we exchanged C_Iberia_CA the Chalcolithic (i.e., Neolithic farmer) ancestry. We added distal with the local SE_Iberia_CA, the model did not find support for a Iran_N (instead of Iran_Chalcolithic, which could contribute other larger number of Iberian BA groups (N_Iberia_BA, SE_Iberia_BA, confounding CA ancestries) as a third source only when needed Almoloya_Argar_Early and Late, and C_Iberia_CA_Stp), and adding (table S2.12 and fig. S7). By using C_Iberia_CA and Yamnaya_Samara a third source in this constellation did not improve the model fit in a two-way model, we could successfully model all Iberian BA either (Fig. 4A and table S2.11). groups, but again to the exclusion of the southeastern Almoloya_ To rule out that the rejection of models involving Germany_Bell_ Argar_Early and Late, and SE_CabezoRedondo_BA groups. Notably, Beaker + C_Iberia_CA in La Almoloya (Early and Late) is due to Bastida_Argar also failed for the distal model 2 (table S2.10). How- specific north/central European LN admixture in Germany_Bell_Beaker ever, these three groups returned values ≥0.05 in the proximal local that was not present in the source of steppe-related ancestry in CA substrate model when Iran_N was added as a third source (fig. Iberia, we fixed the source of direct steppe-related ancestry to S7 and table S2.12). In turn, when we exchanged C_Iberia_CA with Yamnaya_Samara, as we expect the same distal steppe source a local SE_Iberian_CA source, most of the models were rejected, Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 8 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE which confirms the genetic substructure in Iberia_CA groups. Using R1b-Z195 has been found in Sicily_CA_Stp (previously assigned to SE_Iberia_CA as substrate only improves the model fit for south- the EBA and thus considered outliers) and Sicily_EBA (51). How- eastern Iberian groups, suggesting that the first individuals who ever, the genealogical and geographic link to other R1b variants still brought steppe-related ancestry to Iberia admixed locally with the remains unclear. The subtle presence of Iran_N-like ancestry in different CA groups (table S2.12). Notably, SE_Iberia_CA is a suffi- El Argar and Y haplogroup R1b-Z195 in Sicily opens the possibility of ciently supported local proxy for Balearic CA_Stp, BA, and LBA gene flow not only from Iberia to Sicily as previously suggested (51) groups. We interpret this signal with caution, as Mallorca_CA_Stp but also in the opposite direction, implying reciprocal contact with and Menorca_LBA are only represented by a single individual the western and central Mediterranean during the BA, although each; however, this signal is well supported in the newly typed direct contacts are hardly noticeable in the archaeological record. Aritges_LBA group from Menorca, which consists of six individuals Together, these results suggest a dual genetic contribution to the (fig. S7 and table S2.12). formation of the BA genetic profile of southeastern Iberians in ad- We then tried to further specify the proximal sources without dition to a local CA genetic substrate. The major additional ancestry affecting the power of resolution to distinguish between them. We source resembles central European Bell Beaker groups, which first exchanged Yamnaya_Samara by C_Iberia_CA_Stp as the most proxi- contributed ancestry to northern Iberia, followed by a southward mal local source with steppe-related ancestry (6) and combined these spread in the form of C_Iberia_CA_Stp. A second minor ancestry with the different Iberian CA sources (table S2.13 and text S8) in component is an Iran_N-rich/central Mediterranean source, which two-way models, only adding Iran_N as a third source when needed. is restricted to individuals from BA El Argar contexts. The timing of Using increasingly more possible proximal sources for steppe-related the last contributing component remains unclear and points to either ancestry, we could replicate the results of the distal models. However, a Neolithic legacy that persisted throughout the local CA or a subtle we can show that Almoloya_Argar_Early and Late, Bastida_Argar, trace of connections to insular central Mediterranean BA groups. 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. A late Argar genetic outlier makes links to North Africa In contrast, Iran_N is not required when SE_Iberia_CA is used as a and the central Mediterranean local substrate, although models for Almoloya_Argar_Early and Late The PCA also revealed a genetic outlier among the newly typed are still rejected at P ≤ 0.05. This finding is best explained by either individuals (Fig. 3A). The male individual ZAP002 from the late an unresolved distinctive ancestry of La Almoloya individuals or, Argaric site of Lorca-Zapatería falls outside the Iberian_BA cluster alternatively, the ability to detect such subtle signals only in groups on the PCA plot, with a clear shift toward the central Mediterranean. consisting of a larger number of individuals, as is the case for We formally tested for a specific attraction to groups from the eastern Almoloya_Argar_Early (n = 36) and Almoloya_Argar_Late (n = 22), Mediterranean, the Near East, and Africa with an f4-admixture test which provide the statistical power to reject simpler models. When of the form f4(ZAP002, Almoloya_Argar_Early; eastern_Med/Africa, SE_Iberia_CA is used instead of C_Iberia_CA to model La Bastida, Chimp) (Fig. 5A, table S2.16, and text S9). The resulting f4-values were the model is no longer rejected, which suggests a direct contribution consistent with zero, showing that ZAP002 and Almoloya_Argar_Early of local CA groups to those of the southeastern Iberian BA. Together, (the Argar group with the largest number of individuals) are sym- the results hint at an additional, albeit subtle, Iran_N-enriched metrically related to eastern Mediterranean and African groups. ancestry that is present in SE_Iberia_CA groups and thus potentially However, we noticed a near-significant positive f4-value with a predating the El Argar BA. z score = 2.4 for Morocco Iberomaurusian. The attraction to Mbuti We also tested whether Almoloya_Argar_Early and Late can be in an analogous f4-test suggests either African ancestry or high levels modeled when Iran_N is exchanged by geographically more proximal of Basal Eurasian ancestry shared with Morocco Iberomaurusian Iran_N-rich sources involving central and eastern Mediterranean and Natufians (table S2.16). We also confirmed a negative devia- populations (table S2.14 and text S8) and by moving Iran_N to the tion from zero in f4(test, EHG; CHG, Mbuti), indicative of either outgroups. However, this did not result in any supported models African or Basal Eurasian ancestry. Here, we assume that the split of (text S8). While we were not able to find a proximal Mediterranean EHG and CHG groups [Caucasus HG; a basal form of HG ancestry source, we note that adding a central Mediterranean population to south of the Caucasus (72)] is sufficiently deep to result in negative the outgroups (Sicily_EBA, Greece_EBA, or Greece_MBA) de- values for any test population. Following this rationale, the finding creases the model support (P values) for Almoloya_Argar_Early of negative values suggests a high amount of shared deep ancestry and Almoloya_Argar_Late, indirectly attesting to the importance of with Mbuti, which also explains the different f4-value for ZAP002 in central Mediterranean BA. However, we observed no changes in Fig. 3B when compared to other El Argar individuals. We observe P values by adding other Mediterranean groups to the outgroups negative values for ZAP002 (z score = −2.14), as well as Morocco (text S8 and table S2.15). Iberomaurusian and Natufian individuals, who provide additional The complete turnover on the Y chromosome to R1b-Z195 (a support for this claim (Fig. 5B, table S2.17, and text S9). lineage derived from P312), observed in all males at La Almoloya Using qpAdm, we further explored candidate sources of distinct (29 male individuals tested) and La Bastida (except for one E1b ancestries such as Iberian CA and Germany_Bell_Beaker and a lineage dated around 2134 to 1947 cal BCE; 7 males tested), is candidate list of central Mediterranean sources from Sardinia and another independent source of evidence of gene flow predominantly Sicily. Using the same set of outgroups as for the main southeastern during the time of the CA-BA transition (table S1.1 and text S3). Iberian BA groups, we did not find a statistically well-supported Notably, Y lineage R1b-Z195, the most common Y lineage in BA model for ZAP002. However, moving Morocco_Iberomaurusian Iberia, ultimately derives from a common ancestor R1b-P312 in from the outgroups to the sources resulted in a model fit (≥0.05) central Europe but already differs from other derived Bell Beaker when a central Mediterranean population was included as a source. R1b variants reported from central Europe and the British Isles. We found that ZAP002 does not require a local Iberia_CA source Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 9 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE A f4 (ZAP002, Almoloya Argar Early; test, Chimp) B f4 (test, EHG; CHG, Mbuti) TellKurdu EC Sicily EBA Sicily EBA Zapateria Argar Morocco LN Russia Yamnaya Samara Morocco EN Morocco LN Morocco Iberomaurusian Morocco Iberomaurusian England Bell Beaker Morocco EN Germany Bell Beaker Mbuti.DG Germany Corded Ware Jordan PPNB Mallorca CA Stp Formentera MBA Menorca LBA Iran Seh Gabi Chalcolithic C Iberia CA Stp Iran Ganj Dareh Neolithic NW Iberia CA Stp Aritges LBA Ikiztepe LC Bastida Argar Greece Mycenaean Almoloya Argar Late Greece Minoan Almoloya Argar Early SE Iberia BA Argar Ebla EMBA MolinosPapel EBA CamlibelTarlasi LC SW Iberia BA Bastida Argar SE Iberia BA Valencian Cerro Virgen Late Argar Arslantepe LC SE CabezoRedondo BA Arslantepe EBA C Iberia BA CogI Almoloya Argar Late C Iberia BA N Iberia BA Alalakh MLBA NE Iberia BA 0.004 0.002 0 0.002 0.02 0.01 0 3 2 1 0 3 2 1 0 1 Low African attraction High African attraction High Basal Eurasian Low Basal Eurasian Supported model Nonsupported model C Germany Bell Beaker Morocco Iberomaurusian Germany Bell Beaker Morocco Iberomaurusian Sicily EBA C Iberia CA ZAP002 (P = 0.118) ZAP002 (P = 0.002) 0 0.5 1.0 0 0.5 1.0 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). and can be modeled as a three-way mixture of Germany_Bell_Beaker, We found no notable differences in short ROH (<2 cM) patterns Morocco_Iberomaurusian, and either Sicily EBA, Sardinia EBA, or between the two chronocultural periods, indicating that all of the indi- Sardinia Nuragic BA (Fig. 5C and table S2.18). These qpAdm results viduals under study were drawn from populations with comparably and are congruent with the location of ZAP002 in PC space, where he sufficiently large effective population sizes. We also found no evidence clusters with Sicily EBA individuals (Fig. 3A). Rotating the Iberian for a decline of CA population size before the emerging BA. In both CA source to the outgroups in the same model still upholds the model groups, we did not detect signs of inbreeding, represented by long fit (Fig. 5C, table S2.18, and text S9), suggesting an entirely nonlocal ROH segments (20 to 300 cM) (Fig. 6A and table S2.21). origin of individual ZAP002 despite being granted the same burial To assess a potential sex bias in the steppe-related ancestry con- treatment as other individuals from El Argar (pithos burial and tribution as postulated previously (7), we applied distal and proximal Argaric pottery) (Fig. 5C). qpAdm models (text S8) to the X chromosome and autosomes (Fig. 6B and table S2.22). As males contribute, on average, only Insights into phenotypic variation, demography, and social one-third of the X chromosomes to the next generation (75), a lower correlates of CA and EBA El Argar societies proportion of any ancestral component on the X chromosome than We explored the phenotypic variants included in the 1240k SNP panel on the autosomes would thus be indicative of male bias concerning among all well-preserved Iberian CA and BA individuals (>400,000 the respective component (76, 77). In turn, if the ancestral compo- SNPs covered). In addition to 14 variants associated with the ability nent is statistically higher on the X chromosome, this indicates a to digest lactose and alcohol, adaptation to metabolize fatty acids, pre- female bias. On the basis of this rationale, we do not observe signif- disposition to celiac disease, and resistance to infectious diseases, we icant male bias in steppe-related ancestry using either distal or also investigated 41 SNPs encoding for skin/hair and eye pigmenta- proximal sources (Fig. 6B and table S2.22). The fact that the male tion (figs. S8 and S9, tables S2.19 and table S2.20, and text S10). We bias is not detectable could be indicative of an already balanced found no significant increase or decrease in allele frequency of any ancestral component in both sexes, as is reflected in the work of derived allele between the CA and BA (fig. S9 and table S2.20). Mittnik et al. (77), where the male bias in the steppe component To gain insight into the effective population size during the CA and is only detected in Corded Ware, but no longer in Bell Beaker or BA, we applied hapROH to identify runs of homozygosity (73, 74). BA populations. We also took advantage of the recent extensive Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 10 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE A B Russia Yamnaya Samara WHG C Iberia CA Stp C Iberia CA Iran N Anatolia Neolithic Iran N Female bias Male bias Female bias Male bias SE Iberia BA Valencian SE Iberia BA Valencian Almoloya Argar Late Almoloya Argar Late Almoloya Argar Early Almoloya Argar Early N Iberia BA N Iberia BA 4 2 0 2 4 4 2 0 2 4 Sex bias z score Sex bias z score C D 0.252 Mean pairwise mismatch rate (READ) Category 0.250 Female adult Male adult 0.248 ALM052 ALM048 ALM049 ALM034 ALM058 ALM001 ALM043 ALM004 ALM087 ALM015 ALM055 ALM025 ALM032 ALM018 ALM086 ALM028 ALM064 ALM016 ALM071 ALM020 ALM027 ALM051 ALM077 ALM056 ALM053 ALM021 ALM040 ALM017 ALM084 ALM050 ALM024 ALM002 ALM006 ALM038 ALM035 ALM019 ALM008 ALM076 ALM036 ALM003 ALM042 ALM039 ALM007 ALM026 ALM029 ALM044 ALM088 ALM073 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). excavations and sampled all individuals with suitable morphological 30 adult women analyzed). All first-degree relationships among adults preservation (n = 86) from the site La Almoloya, resulting in high-­ involved at least one adult male (three first-degree relationships among quality genome-wide data for 67 individuals that allowed insight into adult men and women). We also found no second-­degree relationships an El Argar community. We estimated the genetic relatedness using a between adult women. The few second-degree relationships involv- combination of the three methods: pairwise mismatch rate (PWMR) ing adults were all between males (4 of 18 adult males analyzed). (48), READ (Relationship Estimation from Ancient DNA) (49), and We then formally tested whether the males of La Almoloya had lcMLkin (47) for pairs of individuals with more than 1000 shared SNPs more close relatives at the site than the females. We calculated (see Materials and Methods, fig. S2, and text S5). All reconstructed pedi- f3-output statistics after removing first- and second-degree pairs of grees will be published in an accompanying paper, fully integrated the form f3(female, female; Mbuti) and f3(male, male; Mbuti) and into the anthropological and archaeological context. Thus, we report found that despite a similar average value in female:female and here the main observations of biological relatedness at the metalevel. male:male comparisons (Fig. 6C and table S2.23), males have a more We found no first-degree relationships between adult women (0 of skewed f3-distribution than females, which suggests a more close Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 11 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE genetic relationships (e.g., at third to fourth degree) among males Laboratory work with other males (P = 1.9 × 10−18) and males with other females All the laboratory work was performed in the dedicated aDNA (P = 2.2 × 10−4) than for females with other females (see Materials and facilities of the Department of Archaeogenetics of the Max Planck Methods). We specifically tested for signals that suggest patrilocality Institute for the Science of Human History in Jena, Germany. using the PWMR extracted from READ (P0_Normalized) (49) and Sampling built a matrix to calculate the biological relatedness of each adult We preferentially sampled petrous bones and teeth (essentially individual to the entire group (Fig. 6D and table S2.24). We observe molars). Details of skeletal elements samples are listed in table S1.1. a tendency that males are, on average, more closely related to other Before sampling, all samples were irradiated with ultraviolet light individuals from the site than females, despite there also being some on either side for 30 min each, after which the surface was cleaned adult females with relatives at the site, which renders the formal test with low-concentration bleach solution (3%). For the petrous bones, for patrilocality nonsignificant (see Materials and Methods). At the we either drilled from the outside, removing the outer layer, or cut intrasite level, this result could also suggest a founder effect from the petrous pyramid longitudinally to drill the dense part directly the male line with a higher reproductive rate. from either side. Molars were cut below the enamel-dentine junc- tion, and tooth powder was drilled from the pulp chamber. We col- lected between 50 and 100 mg of bone or tooth powder per sample DISCUSSION for DNA extractions. Our focused study on the southeastern part of Iberia highlights the Extraction potential of genetic analysis on a regional level in the light of broader-­ DNA was extracted following the modified version of (78) described scale population dynamics in Europe. Overall, we were able to show in www.protocols.io/view/ancient-dna-extraction-from-skeletal-­ that the structure of CA populations shows persistent genetic stability material-baksicwe. and continuity since the Neolithic. This is coupled with a differen- Library preparation tial HG ancestry in the South in addition to potential early contact We generated double-indexed double-stranded (79) UDG (“ds between southeastern Iberian CA and Mediterranean populations, UDG-half”) (39) libraries using 25 l of DNA extract and following which is also revealed by numerous threads of archaeological evidence established protocols (80). Libraries produced per individual are (13–17). We provide nuanced evidence for the arrival of steppe-­ reported in table S1.1. related ancestry in southeastern Iberia by 2200 cal BCE, when Shotgun screening and in-solution enrichment (1240k, MT-capture, widespread social and cultural changes occurred across much of and Y-capture) southern Iberia, including the end of the ditched enclosures system Libraries were sequenced in-house on an Illumina HiSeq 4000 platform of population aggregation and the emergence of BA groups. In the to an average depth of 5 million reads and after demultiplexing pro- Southeast, all sampled BA groups can be shown to carry steppe-­ cessed through EAGER (43). After an initial quality filter based on related ancestry, while the more extensively sampled El Argar groups the presence of aDNA damage and endogenous DNA higher than also carry excess Iran_N-like ancestry, which has also been observed 0.1%, libraries were subsequently enriched using in-solution capture in other CA and BA Mediterranean groups. The analysis of one outlier probes synthetized by Agilent Technologies for ~1240k SNPs along from El Argar reveals a central Mediterranean migrant individual the nuclear genome and independently for the complete mitogenome with additional North African ancestry. Overall, we propose that (40) following (41) and mappable regions of the Y chromosome (42), El Argar has likely formed from a mixture of new groups arriving both in-house protocols. The captured libraries were sequenced for from north-central Iberia, which already carried central European 20 to 40 million reads on average using either a single end (1 × steppe-related ancestry (and the predominant Y-chromosome 75-bp reads) or paired end configuration (2 × 50-bp reads). Results lineage) and local southeastern Iberian CA groups that differed are reported in tables S1.2 and S1.3. from other regions in Iberian in that they carried excess Iran_N-like ancestry similar to eastern and/or central Mediterranean groups. aDNA data processing Moreover, the large sample sizes of early and late El Argar groups Read processing and DNA damage show that the Iran_N-like ancestry contribution in the local CA stratum After demultiplexing based on a specific pair of indexes, raw se- was not enough to explain this type of ancestry satisfyingly, which, in quencing data were processed with EAGER (1.92.59) (43). First, the turn, argues for a continued connection to and genetic influence from adaptor sequences were clipped with Adapter removal (v2.3.1) (81) the Mediterranean BA at least until the end of the El Argar period. and then mapped against the Human Reference Genome hs37d5 with Last, we were able to shed light onto the population structure of Early BWA (v0.7.12) (82) aln and samse commands (-l 16500, -n0.01, -q30). Bronze El Argar societies at the intrasite level. In La Almoloya, a closer We removed duplicates (reads with same start and end and same genetic relationship among males is a strong indicator of patrilineality. orientation) with DeDup (v0.12.2) (43). Last, we determined the de- amination rate pattern in our UDG-half libraries using mapDamage (v.2.0.6) (83). All our libraries showed expected deamination patterns MATERIALS AND METHODS for UDG-half treatment, and according to this, we trimmed 2 bp on Experimental design each terminal side before genotyping. A summary of quality statis- Archaeological samples tics is given in table S1.2. We processed 244 individuals from the Iberian Peninsula and the Genetic sex determination Balearic Islands dated to a time span from the CA (3000 cal BCE) to We determined genetic sex at BAM file level by calculating the cover- the LBA (1200 cal BCE). A detailed description of the archaeological age on the X and Y chromosomes and the autosomes using a bed context, sites, and individuals is reported in text S1 and summarized file of the regions captured by the 1240k SNP array. We normalized in table S1.1. the X and Y reads by the autosomal coverage resulting in the X ratio Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 12 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE and likewise for the Y ratio (table S1.1 and fig. S1). For uncontami- (table S1.1). This procedure allowed us to check in a semi-automatic nated libraries, we expect an X ratio of 1 and a Y ratio of 0 for fe- way which positions were covered to assign an ancestral or derived males and an X and Y ratio of 0.5 for males (46). Individuals that haplogroup or to make corrections to calls in cases where, for in- fall in an intermediate position indicate the presence of DNA stance, a more derived haplogroup was called because of residual contamination from the opposite sex. ancient damage (C to T or G to A mismatches) in terminal read Contamination estimation positions at diagnostic SNPs. As first instance, we inspected the scatterplot of X/autosomes ratio Kinship estimation versus Y/autosomes ratio from the sex determination to look for We calculated the PWMR (48, 87) in all pairs of individuals from outliers that would indicate contamination of the opposite sex. We our pseudo-haploid dataset to double-check for potential duplicate then estimated exogenous human contamination at nuclear and individuals and to determine first- and second-degree relatives. We mitochondrial level. We used method 2 in ANGSD, which estimates split our individuals into two groups according to their genetic heterozygosity at polymorphic sites on the X chromosome in males ancestry profiles (CA and BA individuals), assuming that the median (44), applying a contamination threshold of 3% in individuals with PWMR of the two groups provides a background for unrelated at least 100 X-SNP positions covered twice. In cases where the number individuals (77). We considered one-half of the unrelated value as of SNPs was lower and for all other libraries from female individuals, baseline for identical samples/twins and the midpoint between the we quantified the heterozygosity on the individual mitochondrial two for first-degree relatives (77). After running PWMR, we found reads using ContamMix 1.0.10 (45) by comparing the consensus no duplicates and removed one of each pair of first-degree relatives mitogenome sequence reconstructed from the MT-capture data with for downstream population genetic analyses (in all cases, the indi- a comparative mitochondrial dataset of 311 global mitogenomes. viduals with lower SNP coverage) (table S1.5 and text S5). Only one individual had a rate of contamination of 7.8%, while all For the kinship study of La Almoloya, we compared and com- others were below 1% (table S1.1). We visually inspected the reads bined kinship estimates from several methods. We used READ (49) with Geneious (v9.1.8) (84), looking specifically for background to determine first- and second-degree relatedness among individuals contamination at all variable sites but found no consistent pattern based on the proportion of nonmatching alleles (P0) in nonoverlapping of contamination that would affect the haplotype calls. windows of 1 mega base pairs (Mbps) and calculated standard errors. Genotyping The normalized value of P0 is expected to be between 0 and 0.625 We used our trimmed BAM files to call genotypes with the tool for identical twins, between 0.625 and 0.8125 for first-degree related pilleupCaller (https://github.com/stschiff/sequenceTools), which individuals, between 0.8125 and 0.90625 for second-degree related randomly chooses one allele at every SNP position generating individuals, and greater than 0.90625 for all other degrees of relatedness pseudo-haploid genotypes. We generated genotype files using the (which are reported as “unrelated”) (table S1.6 and text S5). respective SNP lists of the 1240k panel (22) and the Affymetrix HO We also used the method LcMLkin (47), which uses genotype panel with the intersecting SNPs (~600,000 SNPs) (23, 54). Both of likelihoods to estimate the three k-coefficients (k0, k1, or k2) de- them were restricted to only autosomal SNPs for population genetics fined by (88), which define the probability that two individuals have downstream analysis. Numbers of SNPs covered at least once in zero, one, or two alleles identical by descendent (IBD) at a random both HO and 1240k are given in table S1.1. site in the genome. As the method has to deal with missing data, it Mitochondrial haplogroup assignment uses fractions of sites (haplotypes) that are identical by state (IBS). Following targeted enrichment, we could reconstruct complete The length of these analyzed fractions can be modified with the mitochondrial genomes of 144 individuals, with an average coverage thinning parameter --thin. We performed LcMLkin, exploring dif- of 1.3 to 964.6 X (tables S1.3 and S1.4). To process mitochondrial ferent thinning parameters (100,000, 50,000, 10,000, and 5000) to DNA data, we extracted reads from mito-capture data and comple- observe the behavior of first-degree estimates and the potential to mented with reads from 1240k capture and shotgun sequencing when distinguish between possible parent-offspring or sibling relationships necessary using SAMtools v1.3.1 (85) and mapping to the mito- (table S1.7 and text S5). To inform this approach contextually, we chondrial reference, exclusively. Then, we remapped reads to the combined the genetic data with the archaeological information fol- revised Cambridge Reference Sequence (rCRS) using Geneious lowing this order of relevance: (i) We defined the anthropological (v9.1.8) (84). To determine haplogroups from the consensus sequences age at death to determine sibling relatedness among subadults (be- reliably, we applied a threshold of at least 3000 reads per individual fore reaching reproductive age) or to inform on the directionality of (table S1.4). Upon visual inspection of the assembly, we called and parent-offspring relatives, as infants cannot be parents. (ii) When exported the consensus sequences in FASTA format to be used in two unrelated adult male and female both are first degree related to HaploGrep2 v2.1.1 [available via https://haplogrep.uibk.ac.at/ for an (an)other individual(s), we identified these as parent-offspring rela- automated mitochondrial haplogroup assignment based on PhyloTree tionship and used these trios to identify potential cases of half (mtDNA tree build 17, available via www.phylotree.org/)] (86). Us- siblings. (iii) We used the stratigraphic information to confirm the ing Geneious (v9.1.8), we also visually double-checked whether spe- directionality (Who died first?) to establish the type of first-degree cific private SNPs were covered and shared in pairs of individuals relatedness, e.g., assuming that mothers cannot have passed away long with an attested degree of relatedness (table S1.4). before neonates. (iv) We used uniparentally inherited haplogroups Y-haplogroup assignment to corroborate the established pedigrees. By using an in-house protocol for Y capture (42), we were able to Phenotypic and functional variants assign Y-chromosome haplogroups in male individuals. We geno- We explored a list of 56 SNPs associated with known phenotypic typed the Y-chromosome reads using a Y-SNP list from the ISOGG traits. Forty-one of these are linked to physical appearance (skin, eye, (International Society of Genetic Genealogy) dataset included in the or hair color), and we used the HIrisPlex method to calculate the 1240k and Y-capture probes and by using an in-house script (A.B.R.) probabilities of phenotype prediction (fig. S8 and table S2.19) (89). Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 13 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE The other 15 SNPs are associated with the capacity to digest lactose in are inferred on the basis of how the target population is differentially adulthood (rs4988235 in LCT), alcohol dehydrogenase (rs1229984 in related to a set of reference/outgroups via the source populations. ADH1B), adaptation to metabolize fatty acids (rs174546 in FADS1), For all the models applied here, we have used a set of 12 out- predisposition to celiac disease (rs272872 in SLC22A4 and rs653178 groups (Mbuti.DG, Ethiopia_4500BP.SG, Ust_Ishim.DG, Russia_ in ATXN2/SH2B3), or resistance against infectious diseases [rs4833103 MA1_HG.SG, Italy_Villabruna, Belgium_GoyetQ116_1, Han.DG, in the gene cluster TLR1-TLR6-TLR10 and rs2269424 in the major Onge.DG, Papuan.DG, AHG, CHG, and Morocco_Iberomaurusian). histocompatibility complex (MHC)]. We calculated the genotype When we modeled ancestries with proximal sources, we added likelihood to be homozygous or heterozygous for the ancestral (no additional populations or rotated them to the outgroup set. The effect) or derived (effect) allele based on the number of reads from corresponding rationale is explained in each section of the Supple- BAM files after removing duplicates. We restricted the analysis to mentary Materials. individuals with more than 400,000 SNPs covered in the 1240k panel, Admixture dating and results are summarized in fig. S9 and table S2.20. We used DATES v.753 to estimate the date of admixture among different previously fitted ancestral sources (https://github.com/ Population genetic analysis priyamoorjani/DATES). The method works on modeling genotypes Dataset of the admixed test population as a linear combination of two We merged our final dataset with previously published datasets of population sources calculating the ancestry covariance coefficient. ancient and modern individuals reported by the Reich Lab (https:// The results are given in generations (28 years, one generation; reich.hms.harvard.edu/datasets; please see table S2.1 for a detailed table S2.9) (93). list of references) as well as genotypes from recently published studies Sex bias (table S2.1) (50–53). To test for potential sex bias during periods of admixture, we ran the Grouping and labeling individuals distal (model 4; described in text S8) and proximal qpAdm models A detailed description of new labels used in population genetics (models 2 and 4b; described in text S8) for the autosomes and in the analysis is described in text S6 and table S2.1. We only included new X chromosome separately using the parameters allsnps: YES and individuals with more than 40,000 of the 1240k SNPs covered in the chrom23: YES. If there was a sex bias related to any of the modeled downstream population genetics analysis and lowered this threshold sources of ancestry, we expect to observe increasing (female bias) or to >30,000 SNPs only for individuals forming part of a consistent decreasing (male bias) amounts of this ancestry on the X chromo- cluster of individuals. In addition, we excluded groups that have con- some compared to the autosomes, as the paternal contribution on sistently fewer than 60,000 SNPs in qpAdm modeling (table S2.1). X chromosome is about one-third when compared to the autosomes Principal components analysis (the father only contributes one X chromosome to the female We computed PCA using smartpca software from the EIGENSOFT offspring). If the ancestry proportion is significantly higher in the package (v6.0.1) (90) with the lsqproject and SHRINKMODE option autosomes, it implies that the contribution was male driven. Signifi- YES and an extended list of modern populations from Eurasia and cance was calculated with z scores following previous studies (76, 77) North Africa (23) and the Caucasus (91). The newly typed ancient (table S2.22). individuals were projected onto PC1 and PC2. Runs of homozygosity F-statistics We used hapROH to calculate the portion of the genome under runs F-statistics were computed with ADMIXTOOLS (https://github. of homozygosity (73), which can be applied to low-coverage geno- com/DReichLab). For F4-statistics, we used the qpDstat and with the type data. We included all CA and BA individuals from Iberia, which activated f4-mode. F3-statistics were calculated using qp3Pop. We had more than 400,000 SNPs covered in the 1240k panel (table S2.21). used the 1240k panel in all our f-statistics to increase the number of Testing for skewness in f3-statistics SNPs covered by the ancient individuals. Standard errors were cal- To test for a significant change in the skewness of the distributions culated with the default block jacknife. of f3-statistics of the form f3(female, female; Mbuti), f3(male, female; Admixture modeling Mbuti), and f3(male, male; Mbuti), we measured the sample skew- We used ADMIXTURE (60) to define the main genetic cluster pro- ness for each distribution, denoted bff, bmf, and bmm, respectively. files of CA and BA Iberian individuals. We pruned our data for We then looked at the values for the pairwise skewness differences linkage disequilibrium in PLINK (92) with parameters --indep-pairwise of the form i, j = (bi − bj), where i, j ∈ (mm, mf, ff). We estimated the 200 25 0.4 and –maf 0.01, which retained 243,848 autosomal SNPs standard deviation of the i, j, denoted ​​  ˆ ​​​ i,j​​​by taking 1000 bootstrap for the HO dataset. ADMIXTURE was run with default fivefold samples. We then tested the null hypothesis that the i, j = 0 by cal- cross-validation (--cv = 5), varying the number of ancestral popula- culating a z score of the form ​​Z​  i,j​​ = ​​  i,j​​  / ​​ˆ ​​  i,j​​​and then calculated a P tions between K = 2 and K = 18 with different random seeds. Each value using a standard normal distribution. Our P values were ad K was replicated five times. hoc adjusted via the Bonferroni method to account for multiple hy- We used the qpWave and qpAdm programs from the ADMIX- pothesis testing (table S2.23). TOOLS v3.0 package (https://github.com/DReichLab), with the Testing for patrilocality via mean PWMR “allsnps: YES” option to minimally reduce the number of SNPs used To test the overall mean relatedness for each individual, we considered and subsequently increase the power to reject models, to model the the PWMR for each pair of individuals, denoted pij (for individuals ancestry in our new reported individuals. With qpWave, we esti- i and j). We _ then calculated the _ mean PWMR for each individual mated the minimum number of independent contributing ancestral such that p​​​  ​​  i​​ = ​∑ i≠j​ ​​ ​p​ ij​​​. For the p​​​  ​​  i​​​being attributed to groups A and B sources of gene flow needed to explain a target population. With (male and female), we tested for a significant difference in sample qpAdm, we quantified the proportion of genetic ancestry contributed location using a Wilcoxon rank sum test as implemented in the stats by each source. The ancestry proportions in the target population library for the R statistical software package (table S2.24). Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 14 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE SUPPLEMENTARY MATERIALS C. Blasco, C. Liesau, P. Ríos, J. Daura, M. J. de Pedro Michó, A. A. Diez-Castillo, Supplementary material for this article is available at https://science.org/doi/10.1126/ R. Flores Fernández, J. Francès Farré, R. Garrido-Pena, V. S. Gonçalves, E. Guerra-Doce, sciadv.abi7038 A. M. Herrero-Corral, J. Juan-Cabanilles, D. López-Reyes, S. B. McClure, M. Merino Pérez, View/request a protocol for this paper from Bio-protocol. A. Oliver Foix, M. Sanz Borràs, A. C. Sousa, J. M. Vidal Encinas, D. J. Kennett, M. B. Richards, K. Werner Alt, W. Haak, R. Pinhasi, C. Lalueza-Fox, D. Reich, Science 363, 1230–1234 (2019). REFERENCES AND NOTES 8. V. Lull, R. Micó, C. Rihuete Herrada, R. Risch, Transition and Conflict at the End of the 3rd 1. H. Weiss, M.-A. Courty, W. Wetterstrom, F. Guichard, L. Senior, R. Meadow, A. Curnow, The Millennium BC in South Iberia, H. Meller, H. W. Arz, R. Jung, R. Risch, Eds. (Landesmuseum genesis and collapse of third millennium north Mesopotamian civilization. Science 261, für Vorgeschichte, 2015). 995–1004 (1993). 9. R. Chapman, Archaeologies of Complexity (Routledge, 2003). 2. C. Kuzucuoǧlu, C. Marro, Sociétés Humaines Et Changement Climatique À la Fin Du 10. P. B. Ramírez, J. A. Soler Díaz, Ídolos: Miradas Milenarias (Museo Arqueológico de Alicante, Troisième Millénaire: Une Crise A-t-elle Eu Lieu en Haute Mésopotamie?: Actes Du Colloque de 2020). Lyon, 5–8 Décembre 2005 (Institut français d’études anatolienne Georges-Dumézil, 2007). 11. M. D.-Z. Bonilla, J. Beck, H. Bocherens, P. Díaz-del-Río, Isotopic evidence for mobility at 3. H. Meller, H. W. Arz, R. Jung, R. Risch, 2200 BC—A Climatic Breakdown as a Cause for the large-scale human aggregations in Copper Age Iberia: The mega-site of Marroquíes. Collapse of the Old World? (Tagungen des Landesmuseums für Vorgeschichte Halle, Antiquity 92, 991–1007 (2018). 12. Landesamt für Denkmalpflege und Archäologie Sachsen-Anhalt, Landesmuseum für 12. A. C. Valera, I. Žalaitė, A. F. Maurer, V. Grimes, A. M. Silva, S. Ribeiro, J. F. Santos, Vorgeschichte, 2015). C. Barrocas Dias, Addressing human mobility in Iberian Neolithic and Chalcolithic ditched 4. W. Haak, I. Lazaridis, N. Patterson, N. Rohland, S. Mallick, B. Llamas, G. Brandt, enclosures: The case of Perdigões (South Portugal). J. Archaeol. Sci. Rep. 30, 102264 S. Nordenfelt, E. Harney, K. Stewardson, Q. Fu, A. Mittnik, E. Bánffy, C. Economou, (2020). M. Francken, S. Friederich, R. G. Pena, F. Hallgren, V. Khartanovich, A. Khokhlov, M. Kunst, 13. A. Banerjee, J. A. López Padilla, T. X. Schumacher, Marfil y elefantes en la Península Ibérica P. Kuznetsov, H. Meller, O. Mochalov, V. Moiseyev, N. Nicklisch, S. L. Pichler, R. Risch, y el Mediterráneo occidental, Iberia Archaeologica 16.1. (Zabern, 2012). M. A. Rojo Guerra, C. Roth, A. Szécsényi-Nagy, J. Wahl, M. Meyer, J. Krause, D. Brown, 14. T. X. Schuhmacher, Frühbronzezeitliche Kontakte im westlichen und zentralen D. Anthony, A. Cooper, K. W. Alt, D. Reich, Massive migration from the steppe Mittelmeerraum und die Rolle der Iberischen Halbinsel. Madrider Mitteilungen 45, was a source for Indo-European languages in Europe. Nature 522, 207–211 (2015). 147–180 (2004). 5. M. E. Allentoft, M. Sikora, K.-G. Sjögren, S. Rasmussen, M. Rasmussen, J. Stenderup, 15. T. X. Schuhmacher, Elfenbeinobjekte des Chalkolithikums und der Frühen Bronzezeit auf P. B. Damgaard, H. Schroeder, T. Ahlström, L. Vinner, A.-S. Malaspinas, A. Margaryan, der Iberischen Halbinsel. Studien zu Herkunft, Austausch, Verarbeitung und sozialer T. Higham, D. Chivall, N. Lynnerup, L. Harvig, J. Baron, P. D. Casa, P. Dąbrowski, P. R. Duffy, Bedeutung von Elfenbein. Iberia Archaeologica 16, 610 (2012). A. V. Ebel, A. Epimakhov, K. Frei, M. Furmanek, T. Gralak, A. Gromov, S. Gronkiewicz, 16. M. Murillo-Barroso, M. Martinón-Torres, Amber sources and trade in the prehistory G. Grupe, T. Hajdu, R. Jarysz, V. Khartanovich, A. Khokhlov, V. Kiss, J. Kolář, A. Kriiska, of the Iberian peninsula. Eur. J. Archaeol. 15, 187–216 (2012). I. Lasak, C. Longhi, G. McGlynn, A. Merkevicius, I. Merkyte, M. Metspalu, R. Mkrtchyan, 17. R. J. Harrison, A. Gilman, Trade in the second and third millennia BC between the V. Moiseyev, L. Paja, G. Pálfi, D. Pokutta, Ł. Pospieszny, T. D. Price, L. Saag, M. Sablin, Maghreb and Iberia, in Ancient Europe and the Mediterranean: studies in honour of Hugh N. Shishlina, V. Smrčka, V. I. Soenov, V. Szeverényi, G. Tóth, S. V. Trifanova, L. Varul, Hencken, V. Markotic, Ed. (Aris & Phillips, 1977), pp. 90–104. M. Vicze, L. Yepiskoposyan, V. Zhitenev, L. Orlando, T. Sicheritz-Pontén, S. Brunak, 18. V. Hurtado Pérez, L. García Sanjuán, Los inicios de la Jerarquización Social en el Suroeste R. Nielsen, K. Kristiansen, E. Willerslev, Population genomics of Bronze Age Eurasia. Nature de la Península Ibérica (c. 2500–1700 ane): Problemas conceptuales y empíricos (1997); 522, 167–172 (2015). http://roderic.uv.es/bitstream/handle/10550/29498/2134.pdf?sequence=1. 6. I. Olalde, S. Brace, M. E. Allentoft, I. Armit, K. Kristiansen, T. Booth, N. Rohland, S. Mallick, A. Szécsényi-Nagy, A. Mittnik, E. Altena, M. Lipson, I. Lazaridis, T. K. Harper, N. Patterson, 19. R. Chapman, Producing inequalities: Regional sequences in later prehistoric Southern N. Broomandkhoshbacht, Y. Diekmann, Z. Faltyskova, D. Fernandes, M. Ferry, E. Harney, Spain. J. World Prehist. 21, 195–260 (2008). P. de Knijff, M. Michel, J. Oppenheimer, K. Stewardson, A. Barclay, K. W. Alt, C. Liesau, 20. P. Díaz del Río, Labor in the making of Iberian Copper Age lineages (2011); https://digital. P. Ríos, C. Blasco, J. V. Miguel, R. M. García, A. A. Fernández, E. Bánffy, M. Bernabò-Brea, csic.es/handle/10261/35427. D. Billoin, C. Bonsall, L. Bonsall, T. Allen, L. Büster, S. Carver, L. C. Navarro, O. E. Craig, 21. R. Risch, in Überschuss ohne Staat-Politische Formen in der Vorgeschichte: 10. G. T. Cook, B. Cunliffe, A. Denaire, K. E. Dinwiddy, N. Dodwell, M. Ernée, C. Evans, Mitteldeutscher Archäologentag vom 19. bis 21. Oktober 2017 in Halle (Saale) M. Kuchařík, J. F. Farré, C. Fowler, M. Gazenbeek, R. G. Pena, M. Haber-Uriarte, E. Haduch, (Landesmuseum für Vorgeschichte, 2018), pp. 45–65. G. Hey, N. Jowett, T. Knowles, K. Massy, S. Pfrengle, P. Lefranc, O. Lemercier, A. Lefebvre, 22. I. Mathieson, I. Lazaridis, N. Rohland, S. Mallick, N. Patterson, S. A. Roodenberg, E. Harney, C. H. Martínez, V. G. Olmo, A. B. Ramírez, J. L. Maurandi, T. Majó, J. I. McKinley, K. Stewardson, D. Fernandes, M. Novak, K. Sirak, C. Gamba, E. R. Jones, B. Llamas, K. McSweeney, B. G. Mende, A. Modi, G. Kulcsár, V. Kiss, A. Czene, R. Patay, A. Endrődi, S. Dryomov, J. Pickrell, J. L. Arsuaga, J. M. B. de Castro, E. Carbonell, F. Gerritsen, K. Köhler, T. Hajdu, T. Szeniczey, J. Dani, Z. Bernert, M. Hoole, O. Cheronet, D. Keating, A. Khokhlov, P. Kuznetsov, M. Lozano, H. Meller, O. Mochalov, V. Moiseyev, P. Velemínský, M. Dobeš, F. Candilio, F. Brown, R. F. Fernández, A.-M. Herrero-Corral, M. A. R. Guerra, J. Roodenberg, J. M. Vergès, J. Krause, A. Cooper, K. W. Alt, D. Brown, S. Tusa, E. Carnieri, L. Lentini, A. Valenti, A. Zanini, C. Waddington, G. Delibes, D. Anthony, C. Lalueza-Fox, W. Haak, R. Pinhasi, D. Reich, Genome-wide patterns E. Guerra-Doce, B. Neil, M. Brittain, M. Luke, R. Mortimer, J. Desideri, M. Besse, G. Brücken, of selection in 230 ancient Eurasians. Nature 528, 499–503 (2015). M. Furmanek, A. Hałuszko, M. Mackiewicz, A. Rapiński, S. Leach, I. Soriano, K. T. Lillios, 23. I. Lazaridis, D. Nadel, G. Rollefson, D. C. Merrett, N. Rohland, S. Mallick, D. Fernandes, J. L. Cardoso, M. P. Pearson, P. Włodarczak, T. D. Price, P. Prieto, P.-J. Rey, R. Risch, M. Novak, B. Gamarra, K. Sirak, S. Connell, K. Stewardson, E. Harney, Q. Fu, M. A. Rojo Guerra, A. Schmitt, J. Serralongue, A. M. Silva, V. Smrčka, L. Vergnaud, J. Zilhão, G. Gonzalez-Fortes, E. R. Jones, S. A. Roodenberg, G. Lengyel, F. Bocquentin, B. Gasparian, D. Caramelli, T. Higham, M. G. Thomas, D. J. Kennett, H. Fokkens, V. Heyd, A. Sheridan, J. M. Monge, M. Gregg, V. Eshed, A.-S. Mizrahi, C. Meiklejohn, F. Gerritsen, L. Bejenaru, K.-G. Sjögren, P. W. Stockhammer, J. Krause, R. Pinhasi, W. Haak, I. Barnes, C. Lalueza-Fox, M. Blüher, A. Campbell, G. Cavalleri, D. Comas, P. Froguel, E. Gilbert, S. M. Kerr, P. Kovacs, D. Reich, The Beaker phenomenon and the genomic transformation of northwest Europe. J. Krause, D. McGettigan, M. Merrigan, D. A. Merriwether, S. O’Reilly, M. B. Richards, Nature 555, 190–196 (2018). O. Semino, M. Shamoon-Pour, G. Stefanescu, M. Stumvoll, A. Tönjes, A. Torroni, 7. I. Olalde, S. Mallick, N. Patterson, N. Rohland, V. Villalba-Mouco, M. Silva, K. Dulias, J. F. Wilson, L. Yengo, N. A. Hovhannisyan, N. Patterson, R. Pinhasi, D. Reich, Genomic C. J. Edwards, F. Gandini, M. Pala, P. Soares, M. Ferrando-Bernal, N. Adamski, insights into the origin of farming in the ancient Near East. Nature 536, 419–424 (2016). N. Broomandkhoshbacht, O. Cheronet, B. J. Culleton, D. Fernandes, A. M. Lawson, M. Mah, 24. M. Lipson, A. Szécsényi-Nagy, S. Mallick, A. Pósa, B. Stégmár, V. Keerl, N. Rohland, J. Oppenheimer, K. Stewardson, Z. Zhang, J. M. Jiménez Arenas, I. J. Toro Moyano, K. Stewardson, M. Ferry, M. Michel, J. Oppenheimer, N. Broomandkhoshbacht, E. Harney, D. C. Salazar-García, P. Castanyer, M. Santos, J. Tremoleda, M. Lozano, P. García Borja, S. Nordenfelt, B. Llamas, B. Gusztáv Mende, K. Köhler, K. Oross, M. Bondár, T. Marton, J. Fernández-Eraso, J. A. Mujika-Alustiza, C. Barroso, F. J. Bermúdez, E. Viguera Mínguez, A. Osztás, J. Jakucs, T. Paluch, F. Horváth, P. Csengeri, J. Koós, K. Sebők, A. Anders, J. Burch, N. Coromina, D. Vivó, A. Cebrià, J. M. Fullola, O. García-Puchol, J. I. Morales, P. Raczky, J. Regenye, J. P. Barna, S. Fábián, G. Serlegi, Z. Toldi, E. Gyöngyvér Nagy, J. Dani, F. X. Oms, T. Majó, J. M. Vergès, A. Díaz-Carvajal, I. Ollich-Castanyer, F. J. López-Cachero, E. Molnár, G. Pálfi, L. Márk, B. Melegh, Z. Bánfai, L. Domboróczki, J. Fernández-Eraso, A. M. Silva, C. Alonso-Fernández, G. Delibes de Castro, J. Jiménez Echevarría, J. Antonio Mujika-Alustiza, C. Alonso Fernández, J. Jiménez Echevarría, R. Bollongino, A. Moreno-Márquez, G. Pascual Berlanga, P. Ramos-García, J. Ramos-Muñoz, E. Vijande Vila, J. Orschiedt, K. Schierhold, H. Meller, A. Cooper, J. Burger, E. Bánffy, K. W. Alt, G. Aguilella Arzo, Á. Esparza Arroyo, K. T. Lillios, J. Mack, J. Velasco-Vázquez, A. Waterman, C. Lalueza-Fox, W. Haak, D. Reich, Parallel palaeogenomic transects reveal complex L. Benítez de Lugo Enrich, M. Benito Sánchez, B. Agustí, F. Codina, G. de Prado, genetic history of early European farmers. Nature 551, 368–372 (2017). A. Estalrrich, Á. Fernández Flores, C. Finlayson, G. Finlayson, S. Finlayson, F. Giles-Guzmán, 25. R. Martiniano, L. M. Cassidy, R. Ó’Maoldúin, R. McLaughlin, N. M. Silva, L. Manco, A. Rosas, V. Barciela González, G. García Atiénzar, M. S. Hernández Pérez, A. Llanos, D. Fidalgo, T. Pereira, M. J. Coelho, M. Serra, J. Burger, R. Parreira, E. Moran, A. C. Valera, Y. Carrión Marco, I. Collado Beneyto, D. López-Serrano, M. Sanz Tormo, A. C. Valera, E. Porfirio, R. Boaventura, A. M. Silva, D. G. Bradley, The population genomics Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 15 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE of archaeological transition in west Iberia: Investigation of ancient substructure using 46. A. Mittnik, C.-C. Wang, J. Svoboda, J. Krause, A molecular approach to the sexing imputation and haplotype-based methods. PLOS Genet. 13, e1006852 (2017). of the triple burial at the upper paleolithic site of Dolní Věstonice. PLOS ONE 11, e0163019 26. C. Valdiosera, T. Günther, J. C. Vera-Rodríguez, I. Ureña, E. Iriarte, R. Rodríguez-Varela, (2016). L. G. Simões, R. M. Martínez-Sánchez, E. M. Svensson, H. Malmström, L. Rodríguez, 47. M. Lipatov, K. Sanjeev, R. Patro, K. R. Veeramah, Maximum likelihood estimation of J.-M. Bermúdez de Castro, E. Carbonell, A. Alday, J. A. Hernández Vera, A. Götherström, biological relatedness from low coverage sequencing data. bioRxiv 023374 [Preprint]. J.-M. Carretero, J. L. Arsuaga, C. I. Smith, M. Jakobsson, Four millennia of Iberian 29 July 2015. https://doi.org/10.1101/023374. biomolecular prehistory illustrate the impact of prehistoric migrations at the far end 48. D. J. Kennett, S. Plog, R. J. George, B. J. Culleton, A. S. Watson, P. Skoglund, N. Rohland, of Eurasia. Proc. Natl. Acad. Sci. U.S.A. 115, 3428–3433 (2018). S. Mallick, K. Stewardson, L. Kistler, S. A. LeBlanc, P. M. Whiteley, D. Reich, G. H. Perry, 27. S. Brace, Y. Diekmann, T. J. Booth, L. van Dorp, Z. Faltyskova, N. Rohland, S. Mallick, Archaeogenomic evidence reveals prehistoric matrilineal dynasty. Nat. Commun. 8, I. Olalde, M. Ferry, M. Michel, J. Oppenheimer, N. Broomandkhoshbacht, K. Stewardson, 14115 (2017). R. Martiniano, S. Walsh, M. Kayser, S. Charlton, G. Hellenthal, I. Armit, R. Schulting, 49. J. M. Monroy Kuhn, M. Jakobsson, T. Günther, Estimating genetic kin relationships O. E. Craig, A. Sheridan, M. Parker Pearson, C. Stringer, D. Reich, M. G. Thomas, I. Barnes, in prehistoric populations. PLOS ONE 13, e0195491 (2018). Ancient genomes indicate population replacement in Early Neolithic Britain. Nat. Ecol. 50. M. L. Antonio, Z. Gao, H. M. Moots, M. Lucci, F. Candilio, S. Sawyer, V. Oberreiter, Evol. 3, 765–771 (2019). D. Calderon, K. Devitofranceschi, R. C. Aikens, S. Aneli, F. Bartoli, A. Bedini, O. Cheronet, 28. F. Molina González, J. A. Cámara Serrano, J. Capel Martínez, T. Nájera Colino, L. Sáez Pérez, D. J. Cotter, D. M. Fernandes, G. Gasperetti, R. Grifoni, A. Guidi, F. La Pastina, E. Loreti, Los Millares y la periodización de la Prehistoria Reciente del Sudeste, in II - III Simposios de D. Manacorda, G. Matullo, S. Morretta, A. Nava, V. Fiocchi Nicolai, F. Nomi, C. Pavolini, Prehistoria Cueva de Nerja (Fundación Cueva de Nerja, Nerja, 2004), pp. 142–158. M. Pentiricci, P. Pergola, M. Piranomonte, R. Schmidt, G. Spinola, A. Sperduti, M. Rubini, 29. A. C. Valera, Social change in the late 3rd millennium BC in Portugal: The twilight of L. Bondioli, A. Coppa, R. Pinhasi, J. K. Pritchard, Ancient Rome: A genetic crossroads enclosures, in 2200 BC—A Climatic Breakdown as a Cause for the Collapse of the Old World?, of Europe and the Mediterranean. Science 366, 708–714 (2019). H. Meller, R. Risch, R. Jung, H. W. Arz, Eds. (Landesmuseum für vorgeschichte, 2015), 51. D. M. Fernandes, A. Mittnik, I. Olalde, I. Lazaridis, O. Cheronet, N. Rohland, S. Mallick, pp. 409–427. R. Bernardos, N. Broomandkhoshbacht, J. Carlsson, B. J. Culleton, M. Ferry, B. Gamarra, 30. V. Lull, La“ cultura” de El Argar: Un modelo para el estudio de las formaciones económico- M. Lari, M. Mah, M. Michel, A. Modi, M. Novak, J. Oppenheimer, K. A. Sirak, K. Stewardson, sociales prehistóricas (Ediciones Akal, 1983). K. Mandl, C. Schattke, K. T. Özdoğan, M. Lucci, G. Gasperetti, F. Candilio, G. Salis, S. Vai, 31. V. Lull, R. Micó, C. Rihuete Herrada, R. Risch, in Sozialarchäologische Perspektiven: E. Camarós, C. Calò, G. Catalano, M. Cueto, V. Forgia, M. Lozano, E. Marini, M. Micheletti, Gesellschaftlicher Wandel 5000–1500 v. Chr. Zwischen Atlantik und Kaukasus, S. Hansen, R. M. Miccichè, M. R. Palombo, D. Ramis, V. Schimmenti, P. Sureda, L. Teira, M. Teschler-Nicola, J. Müller, Eds. (Archäologie in Eurasien, 2011), vol. 24, pp. 381–414. D. J. Kennett, C. Lalueza-Fox, N. Patterson, L. Sineo, A. Coppa, D. Caramelli, R. Pinhasi, 32. V. Lull, R. Micó, C. Rihuete-Herrada, R. Risch, The La Bastida fortification: New light D. Reich, The spread of steppe and Iranian-related ancestry in the islands of the western and new questions on Early Bronze Age societies in the western Mediterranean. Antiquity Mediterranean. Nat. Ecol. Evol. 4, 334–345 (2020). 88, 395–410 (2014). 52. J. H. Marcus, C. Posth, H. Ringbauer, L. Lai, R. Skeates, C. Sidore, J. Beckett, A. Furtwängler, 33. H. Schubart, Mediterrane Beziehungen der El Argar-Kultur. Madrider Mitteilungen 14, A. Olivieri, C. W. K. Chiang, H. Al-Asadi, K. Dey, T. A. Joseph, C.-C. Liu, C. Der Sarkissian, 41–59 (1973). R. Radzevičiūtė, M. Michel, M. G. Gradoli, P. Marongiu, S. Rubino, V. Mazzarello, D. Rovina, 34. V. Lull, R. Micó, C. Rihuete-Herrada, R. Risch, N. Escanilla, in Actas del Congreso de A. La Fragola, R. M. Serra, P. Bandiera, R. Bianucci, E. Pompianu, C. Murgia, M. Guirguis, Cronometrías Para la Historia de la Península Ibérica (CEUR Workshop Proceedings 2024), R. P. Orquin, N. Tuross, P. van Dommelen, W. Haak, D. Reich, D. Schlessinger, F. Cucca, J. A. Barceló, I. Bogdanovic. B. Morell, Eds. (IberCrono, Barcelona, 2017), pp. 143–162. J. Krause, J. Novembre, Genetic history from the Middle Neolithic to present 35. R. Risch, Recursos naturales, medios de producción y explotación social. Un análisis on the Mediterranean island of Sardinia. Nat. Commun. 11, 939 (2020). económico de la industria lítica de Fuente Álamo (Almería), 2250-1400 ANE. (P. von Zabern. 53. E. Skourtanioti, Y. S. Erdal, M. Frangipane, F. Balossi Restelli, K. A. Yener, F. Pinnock, Mainz, 2002). P. Matthiae, R. Özbal, U.-D. Schoop, F. Guliyev, T. Akhundov, B. Lyonnet, E. L. Hammer, 36. V. Lull, R. M. Pérez, C. R. Herrada, R. Risch, Property relations in the Bronze Age S. E. Nugent, M. Burri, G. U. Neumann, S. Penske, T. Ingman, M. Akar, R. Shafiq, G. Palumbi, of South-western Europe: An archaeological analysis of infant burials from El Argar S. Eisenmann, M. D’Andrea, A. B. Rohrlach, C. Warinner, C. Jeong, P. W. Stockhammer, (Almeria, Spain). Proc. Prehist. Soc. 71, 247–268 (2005). W. Haak, J. Krause, Genomic history of Neolithic to Bronze Age Anatolia, Northern Levant, 37. V. Lull, J. Estévez. Propuesta metodológica para el estudio de las necrópolis argáricas, in and Southern Caucasus. Cell 181, 1158–1175.e28 (2020). Homenaje a Luis Siret (1934–1984), Consejería de Cultura de la Junta de Andalucía 54. N. Patterson, P. Moorjani, Y. Luo, S. Mallick, N. Rohland, Y. Zhan, T. Genschoreck, (Dirección General de Bellas Artes, 1986), pp. 441–452. T. Webster, D. Reich, Ancient admixture in human history. Genetics 192, 1065–1093 (2012). 38. V. Lull, R. Micó, C. Rihuete-Herrada, R. Risch, Funerary practices and kinship in an Early 55. V. Heyd, Kossinna’s smile. Antiquity 91, 348–359 (2017). Bronze Age society: A Bayesian approach applied to the radiocarbon dating of Argaric 56. Q. Fu, C. Posth, M. Hajdinjak, M. Petr, S. Mallick, D. Fernandes, A. Furtwängler, W. Haak, double tombs. J. Archaeol. Sci. 40, 4626–4634 (2013). M. Meyer, A. Mittnik, B. Nickel, A. Peltzer, N. Rohland, V. Slon, S. Talamo, I. Lazaridis, 39. N. Rohland, E. Harney, S. Mallick, S. Nordenfelt, D. Reich, Partial uracil–DNA–glycosylase M. Lipson, I. Mathieson, S. Schiffels, P. Skoglund, A. P. Derevianko, N. Drozdov, treatment for screening of ancient DNA. Philos. Trans. R. Soc. Lond. B Biol. Sci. 370, V. Slavinsky, A. Tsybankov, R. G. Cremonesi, F. Mallegni, B. Gély, E. Vacca, M. R. G. Morales, 20130624 (2015). L. G. Straus, C. Neugebauer-Maresch, M. Teschler-Nicola, S. Constantin, O. T. Moldovan, 40. T. Maricic, M. Whitten, S. Pääbo, Multiplexed DNA sequence capture of mitochondrial S. Benazzi, M. Peresani, D. Coppola, M. Lari, S. Ricci, A. Ronchitelli, F. Valentin, C. Thevenet, genomes using PCR products. PLOS ONE 5, e14004 (2010). K. Wehrberger, D. Grigorescu, H. Rougier, I. Crevecoeur, D. Flas, P. Semal, M. A. Mannino, 41. Q. Fu, M. Hajdinjak, O. T. Moldovan, S. Constantin, S. Mallick, P. Skoglund, N. Patterson, C. Cupillard, H. Bocherens, N. J. Conard, K. Harvati, V. Moiseyev, D. G. Drucker, J. Svoboda, N. Rohland, I. Lazaridis, B. Nickel, B. Viola, K. Prüfer, M. Meyer, J. Kelso, D. Reich, S. Pääbo, M. P. Richards, D. Caramelli, R. Pinhasi, J. Kelso, N. Patterson, J. Krause, S. Pääbo, D. Reich, An early modern human from Romania with a recent Neanderthal ancestor. Nature 524, The genetic history of Ice Age Europe. Nature 534, 200–205 (2016). 216–219 (2015). 57. V. Villalba-Mouco, M. S. van de Loosdrecht, C. Posth, R. Mora, J. Martínez-Moreno, 42. A. B. Rohrlach, L. Papac, A. Childebayeva, M. Rivollat, V. Villalba-Mouco, G. U. Neumann, M. Rojo-Guerra, D. C. Salazar-García, J. I. Royo-Guillén, M. Kunst, H. Rougier, I. Crevecoeur, S. Penske, E. Skourtanioti, M. van de Loosdrecht, M. Akar, K. Boyadzhiev, Y. Boyadzhiev, H. Arcusa-Magallón, C. Tejedor-Rodríguez, I. García-Martínez de Lagrán, R. Garrido-Pena, M.-F. Deguilloux, M. Dobeš, Y. S. Erdal, M. Ernée, M. Frangipane, M. Furmanek, K. W. Alt, C. Jeong, S. Schiffels, P. Utrilla, J. Krause, W. Haak, Survival of Late Pleistocene S. Friederich, E. Ghesquière, A. Haluszko, S. Hansen, M. Küßner, M. Mannino, R. Özbal, hunter-gatherer ancestry in the Iberian peninsula. Curr. Biol. 29, 1169–1177.e7 (2019). S. Reinhold, S. Rottier, D. C. Salazar-García, J. Soler Diaz, P. W. Stockhammer, 58. M. Rivollat, C. Jeong, S. Schiffels, İ. Küçükkalıpçı, M.-H. Pemonge, A. B. Rohrlach, K. W. Alt, C. Roca de Togores Muñoz, K. Aslihan Yener, C. Posth, J. Krause, A. Herbig, W. Haak, Using D. Binder, S. Friederich, E. Ghesquière, D. Gronenborn, L. Laporte, P. Lefranc, H. Meller, Y-chromosome capture enrichment to resolve haplogroup H2 shows new evidence H. Réveillas, E. Rosenstock, S. Rottier, C. Scarre, L. Soler, J. Wahl, J. Krause, M.-F. Deguilloux, for a two-Path Neolithic expansion to Western Europe. Sci. Rep. 11, 15005 (2021). W. Haak, Ancient genome-wide DNA from France highlights the complexity 43. A. Peltzer, G. Jäger, A. Herbig, A. Seitz, C. Kniep, J. Krause, K. Nieselt, EAGER: Efficient of interactions between Mesolithic hunter-gatherers and Neolithic farmers. Sci. Adv. 6, ancient genome reconstruction. Genome Biol. 17, 60 (2016). eaaz5344 (2020). 44. T. S. Korneliussen, A. Albrechtsen, R. Nielsen, ANGSD: Analysis of next generation 59. M. S. van de Loosdrecht, M. A. Mannino, S. Talamo, V. Villalba-Mouco, C. Posth, F. Aron, sequencing data. BMC Bioinformatics 15, 356 (2014). G. Brandt, M. Burri, C. Freund, R. Radzeviciute, R. Stahl, A. Wissgott, L. Klausnitzer, S. Nagel, 45. Q. Fu, A. Mittnik, P. L. F. Johnson, K. Bos, M. Lari, R. Bollongino, C. Sun, L. Giemsch, M. Meyer, A. Tagliacozzo, M. Piperno, S. Tusa, C. Collina, V. Schimmenti, R. D. Salvo, R. Schmitz, J. Burger, A. M. Ronchitelli, F. Martini, R. G. Cremonesi, J. Svoboda, P. Bauer, K. Prüfer, J.-J. Hublin, S. Schiffels, C. Jeong, W. Haak, J. Krause, Genomic and dietary D. Caramelli, S. Castellano, D. Reich, S. Pääbo, J. Krause, A revised timescale for human transitions during the Mesolithic and Early Neolithic in Sicily. bioRxiv 2020.03.11.986158 evolution based on ancient mitochondrial genomes. Curr. Biol. 23, 553–559 (2013). [Preprint]. 12 March 2020; https://doi.org/10.1101/2020.03.11.986158. Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 16 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE 60. D. H. Alexander, J. Novembre, K. Lange, Fast model-based estimation of ancestry C. Thevenet, I. Valchev, F. Valentin, S. Vasilyev, F. Veljanovska, S. Venelinova, in unrelated individuals. Genome Res. 19, 1655–1664 (2009). E. Veselovskaya, B. Viola, C. Virag, J. Zaninović, S. Zäuner, P. W. Stockhammer, G. Catalano, 61. J. A. Alcover, The first Mallorcans: Prehistoric colonization in the Western Mediterranean. R. Krauß, D. Caramelli, G. Zariņa, B. Gaydarska, M. Lillie, A. G. Nikitin, I. Potekhina, J. World Prehist. 21, 19–84 (2008). A. Papathanasiou, D. Borić, C. Bonsall, J. Krause, R. Pinhasi, D. Reich, The genomic history 62. R. Micó, Radiocarbon dating and Balearic prehistory: Reviewing the periodization of southeastern Europe. Nature 555, 197–203 (2018). of the prehistoric sequence. Radiocarbon 48, 421–434 (2006). 77. A. Mittnik, K. Massy, C. Knipper, F. Wittenborn, R. Friedrich, S. Pfrengle, M. Burri, 63. L. García Sanjuán, J. M. Vargas Jiménez, L. M. Cáceres Puro, M. E. Costa Caramé, N. Carlichi-Witjes, H. Deeg, A. Furtwängler, M. Harbeck, K. von Heyking, C. Kociumaka, M. Díaz-Guardamino Uribe, M. Díaz-Zorita Bonilla, Á. Fernández Flores, V. Hurtado Pérez, I. Kucukkalipci, S. Lindauer, S. Metz, A. Staskiewicz, A. Thiel, J. Wahl, W. Haak, E. Pernicka, P. M. López Aldana, E. Méndez Izquierdo, A. Pajuelo Pando, J. Rodríguez Vidal, S. Schiffels, P. W. Stockhammer, J. Krause, Kinship-based social inequality in Bronze Age D. Wheatley, C. Bronk Ramsey, A. Delgado-Huertas, E. Dunbar, A. Mora González, Europe. Science 366, 731–734 (2019). A. Bayliss, N. Beavan, D. Hamilton, A. Whittle, Assembling the dead, gathering the living: 78. J. Dabney, M. Knapp, I. Glocke, M.-T. Gansauge, A. Weihmann, B. Nickel, C. Valdiosera, Radiocarbon dating and Bayesian modelling for Copper Age Valencina de la Concepción N. García, S. Pääbo, J.-L. Arsuaga, M. Meyer, Complete mitochondrial genome sequence (Seville, Spain). J. World Prehist. 31, 179–313 (2018). of a Middle Pleistocene cave bear reconstructed from ultrashort DNA fragments. Proc. 64. G. Aranda Jiménez, M. Díaz-Zorita Bonilla, D. Hamilton, L. Milesi, M. Sánchez Romero, The Natl. Acad. Sci. U.S.A. 110, 15758–15763 (2013). radiocarbon chronology and temporality of the megalithic cemetery of Los Millares 79. M. Kircher, S. Sawyer, M. Meyer, Double indexing overcomes inaccuracies in multiplex (Almería, Spain). Archaeol. Anthropol. Sci. 12, 104 (2020). sequencing on the Illumina platform. Nucleic Acids Res. 40, e3 (2012). 65. A. C. Valera, Social change in the late 3rd millennium BC in Portugal: The twilight of enclosures, 80. M. Meyer, M. Kircher, Illumina sequencing library preparation for highly multiplexed in 2200 BC: A Climatic Break-down as a Cause for the Collapse of the Old World, H. Meller, target capture and sequencing. Cold Spring Harb. Protoc. 2010, pdb.prot5448 (2010). H. Arz, R. Jung, R. Risch, Eds. (Landesmuseum für Vorgeschichte, Halle, 2014), 81. M. Schubert, S. Lindgreen, L. Orlando, AdapterRemoval v2: Rapid adapter trimming, pp. 409–428. identification, and read merging. BMC Res. Notes 9, 88 (2016). 66. A. Blanco-González, K. T. Lillios, J. A. López-Sáez, B. L. Drake, Cultural, demographic 82. H. Li, R. Durbin, Fast and accurate short read alignment with Burrows–Wheeler transform. and environmental dynamics of the Copper and Early Bronze Age in Iberia (3300–1500 Bioinformatics 25, 1754–1760 (2009). BC): Towards an interregional multiproxy comparison at the time of the 4.2 ky BP event. 83. A. Ginolhac, H. Jónsson, L. Orlando, M. Schubert, P. L. F. Johnson, mapDamage2.0: Fast J. World Prehist. 31, 1–79 (2018). approximate Bayesian estimates of ancient DNA damage parameters. Bioinformatics 29, 67. M. Hinz, J. Schirrmacher, J. Kneisel, C. Rinne, M. Weinelt, The Chalcolithic–Bronze Age 1682–1684 (2013). transition in southern Iberia under the influence of the 4.2 kyr event? A correlation 84. M. Kearse, R. Moir, A. Wilson, S. Stones-Havas, M. Cheung, S. Sturrock, S. Buxton, of climatological and demographic proxies. J. Neolithic Archaeol. 21, 1–26 (2019). A. Cooper, S. Markowitz, C. Duran, T. Thierer, B. Ashton, P. Meintjes, A. Drummond, 68. J. Schirrmacher, J. Kneisel, D. Knitter, W. Hamer, M. Hinz, R. R. Schneider, M. Weinelt, Geneious Basic: An integrated and extendable desktop software platform Spatial patterns of temperature, precipitation, and settlement dynamics on the Iberian for the organization and analysis of sequence data. Bioinformatics 28, 1647–1649 Peninsula during the Chalcolithic and the Bronze Age. Quat. Sci. Rev. 233, 106220 (2012). (2020). 85. H. Li, B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer, G. Marth, G. Abecasis, 69. K. T. Lillios, A. Blanco-González, B. L. Drake, J. A. López-Sáez, Mid-late Holocene climate, R. Durbin; 1000 Genome Project Data Processing Subgroup, The sequence alignment/ demography, and cultural dynamics in Iberia: A multi-proxy approach. Quat. Sci. Rev. 135, Map format and SAMtools. Bioinformatics 25, 2078–2079 (2009). 138–153 (2016). 86. H. Weissensteiner, D. Pacher, A. Kloss-Brandstätter, L. Forer, G. Specht, H.-J. Bandelt, F. Kronenberg, A. Salas, S. Schönherr, HaploGrep 2: Mitochondrial haplogroup 70. V. Lull, R. Micó, C. Rihuete, R. Risch, in The Matter of Prehistory. Papers in Honour of Antonio classification in the era of high-throughput sequencing. Nucleic Acids Res. 44, W58–W63 Gilman Guillén, P. D. del Río, K. Lillios, I. Sastre, Eds. (Biblioteca Praehistorica Hispana, (2016). 36—Consejo Superior de Investigaciones Científicas, 2020), pp. 191–209. 87. M. van de Loosdrecht, A. Bouzouggar, L. Humphrey, C. Posth, N. Barton, A. Aximu-Petri, 71. I. Lazaridis, A. Mittnik, N. Patterson, S. Mallick, N. Rohland, S. Pfrengle, A. Furtwängler, B. Nickel, S. Nagel, E. H. Talbi, M. A. El Hajraoui, S. Amzazi, J.-J. Hublin, S. Pääbo, S. Schiffels, A. Peltzer, C. Posth, A. Vasilakis, P. J. P. McGeorge, E. Konsolaki-Yannopoulou, G. Korres, M. Meyer, W. Haak, C. Jeong, J. Krause, Pleistocene North African genomes link Near H. Martlew, M. Michalodimitrakis, M. Özsait, N. Özsait, A. Papathanasiou, M. Richards, Eastern and sub-Saharan African human populations. Science 360, 548–552 (2018). S. A. Roodenberg, Y. Tzedakis, R. Arnott, D. M. Fernandes, J. R. Hughey, D. M. Lotakis, 88. C. C. Cockerham, Higher order probability functions of identity of allelles by descent. P. A. Navas, Y. Maniatis, J. A. Stamatoyannopoulos, K. Stewardson, P. Stockhammer, Genetics 69, 235–246 (1971). R. Pinhasi, D. Reich, J. Krause, G. Stamatoyannopoulos, Genetic origins of the Minoans 89. S. Walsh, F. Liu, A. Wollstein, L. Kovatsi, A. Ralf, A. Kosiniak-Kamysz, W. Branicki, M. Kayser, and Mycenaeans. Nature 548, 214–218 (2017). The HIrisPlex system for simultaneous prediction of hair and eye colour from DNA. 72. E. R. Jones, G. Gonzalez-Fortes, S. Connell, V. Siska, A. Eriksson, R. Martiniano, Forensic Sci. Int. Genet. 7, 98–115 (2013). R. L. McLaughlin, M. Gallego Llorente, L. M. Cassidy, C. Gamba, T. Meshveliani, 90. N. Patterson, A. L. Price, D. Reich, Population structure and eigenanalysis. PLOS Genet. 2, O. Bar-Yosef, W. Müller, A. Belfer-Cohen, Z. Matskevich, N. Jakeli, T. F. G. Higham, e190 (2006). M. Currat, D. Lordkipanidze, M. Hofreiter, A. Manica, R. Pinhasi, D. G. Bradley, Upper Palaeolithic genomes reveal deep roots of modern Eurasians. Nat. Commun. 6, 8912 91. C.-C. Wang, S. Reinhold, A. Kalmykov, A. Wissgott, G. Brandt, C. Jeong, O. Cheronet, (2015). M. Ferry, E. Harney, D. Keating, S. Mallick, N. Rohland, K. Stewardson, A. R. Kantorovich, V. E. Maslov, V. G. Petrenko, V. R. Erlikh, B. C. Atabiev, R. G. Magomedov, P. L. Kohl, 73. H. Ringbauer, J. Novembre, M. Steinrücken, Parental relatedness through time revealed K. W. Alt, S. L. Pichler, C. Gerling, H. Meller, B. Vardanyan, L. Yeganyan, A. D. Rezepkin, by runs of homozygosity in ancient DNA. Nat. Commun. 12, 5425 (2021). D. Mariaschk, N. Berezina, J. Gresky, K. Fuchs, C. Knipper, S. Schiffels, E. Balanovska, 74. H. Ringbauer, M. Steinrücken, L. Fehren-Schmitz, D. Reich, Increased rate of close-kin O. Balanovsky, I. Mathieson, T. Higham, Y. B. Berezin, A. Buzhilova, V. Trifonov, R. Pinhasi, unions in the central Andes in the half millennium before European contact. Curr. Biol. 30, A. B. Belinskij, D. Reich, S. Hansen, J. Krause, W. Haak, Ancient human genome-wide data R980–R981 (2020). from a 3000-year interval in the Caucasus corresponds with eco-geographic regions. 75. A. Goldberg, N. A. Rosenberg, Beyond 2/3 and 1/3: The complex signatures of sex-biased Nat. Commun. 10, 590 (2019). admixture on the X chromosome. Genetics 201, 263–279 (2015). 92. S. Purcell, B. Neale, K. Todd-Brown, L. Thomas, M. A. R. Ferreira, D. Bender, J. Maller, 76. I. Mathieson, S. Alpaslan-Roodenberg, C. Posth, A. Szécsényi-Nagy, N. Rohland, S. Mallick, P. Sklar, P. I. W. de Bakker, M. J. Daly, P. C. Sham, PLINK: A tool set for whole-genome I. Olalde, N. Broomandkhoshbacht, F. Candilio, O. Cheronet, D. Fernandes, M. Ferry, association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 B. Gamarra, G. G. Fortes, W. Haak, E. Harney, E. Jones, D. Keating, B. Krause-Kyora, (2007). I. Kucukkalipci, M. Michel, A. Mittnik, K. Nägele, M. Novak, J. Oppenheimer, N. Patterson, 93. P. Moorjani, S. Sankararaman, Q. Fu, M. Przeworski, N. Patterson, D. Reich, A genetic S. Pfrengle, K. Sirak, K. Stewardson, S. Vai, S. Alexandrov, K. W. Alt, R. Andreescu, method for dating ancient genomes provides a direct estimate of human generation D. Antonović, A. Ash, N. Atanassova, K. Bacvarov, M. B. Gusztáv, H. Bocherens, M. Bolus, interval in the last 45,000 years. Proc. Natl. Acad. Sci. U.S.A. 113, 5652–5657 (2016). A. Boroneanţ, Y. Boyadzhiev, A. Budnik, J. Burmaz, S. Chohadzhiev, N. J. Conard, 94. V. Lull, C. Rihuete Herrada, R. Risch, E. Celdrán Beltrán, M. I. Fregeiro Morador, C. Oliart R. Cottiaux, M. Čuka, C. Cupillard, D. G. Drucker, N. Elenski, M. Francken, B. Galabova, Caravatti, C. Velasco Felipe, La Almoloya (Pliego-Murcia): Murcia: Integral, Sociedad Para el G. Ganetsovski, B. Gély, T. Hajdu, V. Handzhyiska, K. Harvati, T. Higham, S. Iliev, I. Janković, Desarrollo Rural (Integral - Sociedad para el Desarrollo Rural, Bullas, Spain, 2015). I. Karavanić, D. J. Kennett, D. Komšo, A. Kozak, D. Labuda, M. Lari, C. Lazar, M. Leppek, 95. V. Lull, C. Rihuete Herrada, R. Risch, B. Bonora, E. Celdrán Beltrán, M. I. Fregeiro Morador, K. Leshtakov, D. L. Vetro, D. Los, I. Lozanov, M. Malina, F. Martini, K. McSweeney, H. Meller, C. Molero, C. Moreno, C. Oliart Caravatti, C. Velasco-Felipe, L. Andúgar, W. Haak, M. Menđušić, P. Mirea, V. Moiseyev, V. Petrova, T. D. Price, A. Simalcsik, L. Sineo, M. Šlaus, V. Villalba-Mouco, R. Micó, Emblems and spaces of power during the Argaric Bronze Age V. Slavchev, P. Stanev, A. Starović, T. Szeniczey, S. Talamo, M. Teschler-Nicola, at La Almoloya, Murcia. Antiquity 95, 329–348 (2021). Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 17 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE 96. V. Lull, R. Micó, C. Rihuete Herrada, R. Risch, La Bastida y Tira del Lienzo (Totana-Murcia): 108. V. M. Narasimhan, N. Patterson, P. Moorjani, N. Rohland, R. Bernardos, S. Mallick, Murcia: Integral, Sociedad Para el Desarrollo Rural (Integral - Sociedad para el Desarrollo I. Lazaridis, N. Nakatsuka, I. Olalde, M. Lipson, A. M. Kim, L. M. Olivieri, A. Coppa, Rural, Bullas, Spain, 2015). M. Vidale, J. Mallory, V. Moiseyev, E. Kitov, J. Monge, N. Adamski, N. Alex, 97. C. Knipper, C. Rihuete-Herrada, J. Voltas, P. Held, V. Lull, R. Micó, R. Risch, K. W. Alt, N. Broomandkhoshbacht, F. Candilio, K. Callan, O. Cheronet, B. J. Culleton, M. Ferry, Reconstructing Bronze Age diets and farming strategies at the early Bronze Age sites D. Fernandes, S. Freilich, B. Gamarra, D. Gaudio, M. Hajdinjak, É. Harney, T. K. Harper, of La Bastida and Gatas (southeast Iberia) using stable isotope analysis. PLOS ONE 15, D. Keating, A. M. Lawson, M. Mah, K. Mandl, M. Michel, M. Novak, J. Oppenheimer, e0229398 (2020). N. Rai, K. Sirak, V. Slon, K. Stewardson, F. Zalzala, Z. Zhang, G. Akhatov, A. N. Bagashev, 98. A. Szécsényi-Nagy, C. Roth, G. Brandt, C. Rihuete-Herrada, C. Tejedor-Rodríguez, P. Held, A. Bagnera, B. Baitanayev, J. Bendezu-Sarmiento, A. A. Bissembaev, G. L. Bonora, Í. García-Martínez-de-Lagrán, H. Arcusa Magallón, S. Zesch, C. Knipper, E. Bánffy, T. T. Chargynov, T. Chikisheva, P. K. Dashkovskiy, A. Derevianko, M. Dobeš, K. Douka, S. Friederich, H. Meller, P. Bueno Ramírez, R. Barroso Bermejo, R. de Balbín Behrmann, N. Dubova, M. N. Duisengali, D. Enshin, A. Epimakhov, A. V. Fribus, D. Fuller, A. M. Herrero-Corral, R. Flores Fernández, C. Alonso Fernández, J. Jiménez Echevarria, A. Goryachev, A. Gromov, S. P. Grushin, B. Hanks, M. Judd, E. Kazizov, A. Khokhlov, L. Rindlisbacher, C. Oliart, M.-I. Fregeiro, I. Soriano, O. Vicente, R. Micó, V. Lull, J. Soler Díaz, A. P. Krygin, E. Kupriyanova, P. Kuznetsov, D. Luiselli, F. Maksudov, A. M. Mamedov, J. A. López Padilla, C. Roca de Togores Muñoz, M. S. Hernández Pérez, F. J. Jover Maestre, T. B. Mamirov, C. Meiklejohn, D. C. Merrett, R. Micheli, O. Mochalov, S. Mustafokulov, J. Lomba Maurandi, A. Avilés Fernández, K. T. Lillios, A. M. Silva, M. Magalhães Ramalho, A. Nayak, D. Pettener, R. Potts, D. Razhev, M. Rykun, S. Sarno, T. M. Savenkova, L. M. Oosterbeek, C. Cunha, A. J. Waterman, J. Roig Buxó, A. Martínez, J. Ponce Martínez, K. Sikhymbaeva, S. M. Slepchenko, O. A. Soltobaev, N. Stepanova, S. Svyatko, M. Hunt Ortiz, J. C. Mejías-García, J. C. Pecero Espín, R. Cruz-Auñón Briones, T. Tomé, K. Tabaldiev, M. Teschler-Nicola, A. A. Tishkin, V. V. Tkachev, S. Vasilyev, P. Velemínský, E. Carmona Ballestero, J. L. Cardoso, A. C. Araújo, C. Liesau von Lettow-Vorbeck, D. Voyakin, A. Yermolayeva, M. Zahir, V. S. Zubkov, A. Zubova, V. S. Shinde, C. Lalueza-Fox, C. Blasco Bosqued, P. Ríos Mendoza, A. Pujante, J. I. Royo-Guillén, M. A. Esquembre Beviá, M. Meyer, D. Anthony, N. Boivin, K. Thangaraj, D. J. Kennett, M. Frachetti, R. Pinhasi, V. M. Dos Santos Goncalves, R. Parreira, E. Morán Hernández, E. Méndez Izquierdo, D. Reich, The formation of human populations in South and Central Asia. Science 365, J. Vega y Miguel, R. Menduiña García, V. Martínez Calvo, O. López Jiménez, J. Krause, eaat7487 (2019). S. L. Pichler, R. Garrido-Pena, M. Kunst, R. Risch, M. A. Rojo-Guerra, W. Haak, K. W. Alt, The 109. H. F. Klinefelter Jr., E. C. Reifenstein Jr., Syndrome characterized by gynecomastia, maternal genetic make-up of the Iberian Peninsula between the Neolithic and the Early aspermatogenesis without A-Leydigism, and increased excretion of follicle-stimulating Bronze Age. Sci. Rep. 7, 15644 (2017). hormone. J. Clin. Endocrinol. 2, 615–627 (1942). 99. J. A. Soler Díaz, C. Roca de Togores, Ritual funerario en la Cova d’En Pardo ca. 3.350–2.850 110. P. Jacobs, A. G. Baikie, T. N. MacGregor, N. Maclean, D. G. Harnden, Evidence CAL ANE: Espacialidad, cronología y territorio cultural, in Cova d’En Pardo: Arqueología en la for the existence of the human “super female.” Lancet 274, 423–425 (1959). Memoria: Excavaciones de M. Tarradell, V. Pascual y E. Llobregat (1961–1965), Catálogo de 111. R. Boada, J. Janusz, C. Hutaff-Lee, N. Tartaglia, The cognitive phenotype in Klinefelter Materiales del Museo de Alcoy y Estudios a Partir de las Campañas del MARQ (1993–2007) syndrome: A review of the literature including genetic and hormonal factors. Dev. Disabil. en la Cavidad de Planes, Alicante. (Museo Arqueológico de Alicante­—MARQ, 2012), Res. Rev. 15, 284–294 (2009). pp. 205–248. 112. S. S. Ebenesersdóttir, M. Sandoval-Velasco, E. D. Gunnarsdóttir, A. Jagadeesan, 100. J. A. S. Díaz, C. Roca de Togores Muñoz, Estudio de los restos humanos encontrados en V. B. Guðmundsdóttir, E. L. Thordardóttir, M. S. Einarsdóttir, K. H. S. Moore, las intervenciones practicadas en los años 1961 y 1965 en la Cova d’En Pardo, Planes, Á. Sigurðsson, D. N. Magnúsdóttir, H. Jónsson, S. Snorradóttir, E. Hovig, P. Møller, Alicante: análisis antropológico y aproximación a su contexto cultural. Saguntum Extra 2 I. Kockum, T. Olsson, L. Alfredsson, T. F. Hansen, T. Werge, G. L. Cavalleri, E. Gilbert, (1999), pp. 369–377. C. Lalueza-Fox, J. W. Walser III, S. Kristjánsdóttir, S. Gopalakrishnan, L. Árnadóttir, 101. C. Roca de Togores Muñoz, J. A. S. Díaz, Restos humanos en la Cova d’en Pardo (Planes). Ó. Þ. Magnússon, M. T. P. Gilbert, K. Stefánsson, A. Helgason, Ancient genomes Problemática y avance de resultados de la investigación antropológica en una cavidad from Iceland reveal the making of a human population. Science 360, 1028–1032 de inhumación múltiple excavada en dos etapas: 1961-1965 y 1993-2007, in Cova d’En (2018). Pardo: Arqueología en la Memoria: Excavaciones de M. Tarradell, V. Pascual y E. Llobregat 113. M. Otter, C. T. Schrander-Stumpel, L. M. G. Curfs, Triple X syndrome: A review (1961–1965), Catálogo de Materiales del Museo de Alcoy y Estudios a Partir de las Campañas of the literature. Eur. J. Hum. Genet. 18, 265–271 (2010). del MARQ (1993–2007) en la Cavidad de Planes, Alicante (Museo Arqueológico de 114. M. Rafique, S. AlObaid, D. Al-Jaroudi, 47, XXX syndrome with infertility, premature ovarian Alicante—MARQ, 2012), pp. 193–204. insufficiency, and streak ovaries. Clin. Case Rep. 7, 1238–1241 (2019). 102. V. Lull, R. Micó, C. Rihuete, R. Risch, La Cova des Càrritx y la Cova des Mussol. Ideología y 115. H. Schroeder, A. Margaryan, M. Szmyt, B. Theulot, P. Włodarczak, S. Rasmussen, Sociedad en la Prehistoria de de Menorca (Consell Insular de Menorca, 1999). S. Gopalakrishnan, A. Szczepanek, T. Konopka, T. Z. T. Jensen, B. Witkowska, S. Wilk, 103. M. De Cet, Long-term Social Development on a Mediterranean Island: Menorca Between M. M. Przybyła, Ł. Pospieszny, K.-G. Sjögren, Z. Belka, J. Olsen, K. Kristiansen, E. Willerslev, 1600 BCE and 1900 CE (In Kommission bei Verlag Dr. Rudolf Habelt GmbH, 2017). K. M. Frei, M. Sikora, N. N. Johannsen, M. E. Allentoft, Unraveling ancestry, kinship, 104. T. Günther, C. Valdiosera, H. Malmström, I. Ureña, R. Rodriguez-Varela, Ó. O. Sverrisdóttir, and violence in a Late Neolithic mass grave. Proc. Natl. Acad. Sci. U.S.A. 116, 10705–10710 E. A. Daskalaki, P. Skoglund, T. Naidoo, E. M. Svensson, J. M. Bermúdez de Castro, (2019). E. Carbonell, M. Dunn, J. Storå, E. Iriarte, J. L. Arsuaga, J.-M. Carretero, A. Götherström, 116. S. Walsh, L. Chaitanya, K. Breslin, C. Muralidharan, A. Bronikowska, E. Pospiech, J. Koller, M. Jakobsson, Ancient genomes link early farmers from Atapuerca in Spain to modern- L. Kovatsi, A. Wollstein, W. Branicki, F. Liu, M. Kayser, Global skin colour prediction day Basques. Proc. Natl. Acad. Sci. U.S.A. 112, 11917–11922 (2015). from DNA. Hum. Genet. 136, 847–863 (2017). 105. P. de Barros Damgaard, N. Marchi, S. Rasmussen, M. Peyrot, G. Renaud, T. Korneliussen, 117. L. Chaitanya, K. Breslin, S. Zuñiga, L. Wirken, E. Pośpiech, M. Kukla-Bartoszek, T. Sijen, J. V. Moreno-Mayar, M. W. Pedersen, A. Goldberg, E. Usmanova, N. Baimukhanov, P. de Knijff, F. Liu, W. Branicki, M. Kayser, S. Walsh, The HIrisPlex-S system for eye, hair V. Loman, L. Hedeager, A. G. Pedersen, K. Nielsen, G. Afanasiev, K. Akmatov, A. Aldashev, and skin colour prediction from DNA: Introduction and forensic developmental A. Alpaslan, G. Baimbetov, V. I. Bazaliiskii, A. Beisenov, B. Boldbaatar, B. Boldgiv, validation. Forensic Sci. Int. Genet. 35, 123–135 (2018). C. Dorzhu, S. Ellingvag, D. Erdenebaatar, R. Dajani, E. Dmitriev, V. Evdokimov, K. M. Frei, A. Gromov, A. Goryachev, H. Hakonarson, T. Hegay, Z. Khachatryan, R. Khaskhanov, E. Kitov, A. Kolbina, T. Kubatbek, A. Kukushkin, I. Kukushkin, N. Lau, A. Margaryan, I. Merkyte, I. V. Mertz, V. K. Mertz, E. Mijiddorj, V. Moiyesev, G. Mukhtarova, Acknowledgments: We thank all members of the Archaeogenetics Department of the Max B. Nurmukhanbetov, Z. Orozbekova, I. Panyushkina, K. Pieta, V. Smrčka, I. Shevnina, Planck Institute for the Science of Human History, especially the Population Genetics and A. Logvin, K.-G. Sjögren, T. Štolcová, A. M. Taravella, K. Tashbaeva, A. Tkachev, PALEORIDER groups, for continued support and discussions. We thank G. Brandt, T. Tulegenov, D. Voyakin, L. Yepiskoposyan, S. Undrakhbold, V. Varfolomeev, A. Weber, A. Wissgott, S. Clayton, and K. Prüfer for sequencing and handling of the raw data; M. A. W. Sayres, N. Kradin, M. E. Allentoft, L. Orlando, R. Nielsen, M. Sikora, E. Heyer, A. Ghalichi for help with ADMIXTURE; and S. O’Reilly for the design of Fig. 1. We are indebted K. Kristiansen, E. Willerslev, 137 ancient human genomes from across the Eurasian to the Museo de Lorca, Museo de Alcoy, Museo de Villena, Museo de Sevilla, and Museo de steppes. Nature 557, 369–374 (2018). Alicante and all archaeologists and students involved in the excavations. Funding: This work 106. A. Mittnik, C.-C. Wang, S. Pfrengle, M. Daubaras, G. Zariņa, F. Hallgren, R. Allmäe, was supported by the Max Planck Society (V.V.-M. and W.H.); European Research Council V. Khartanovich, V. Moiseyev, M. Tõrv, A. Furtwängler, A. Andrades Valtueña, M. Feldman, (ERC) grant 771234—PALEoRIDER (W.H.); Spanish Ministry of Economy, Industry and C. Economou, M. Oinonen, A. Vasks, E. Balanovska, D. Reich, R. Jankauskas, W. Haak, S. Schiffels, Competitiveness project HAR2017-85962-P (C.O., C.R.-H., M.I.F., E.C.B., C.V.-F., V.L., R.M., and J. Krause, The genetic prehistory of the Baltic Sea region. Nat. Commun. 9, 442 (2018). R.R.); AGAUR 2017SGR1044 (C.O., C.R.-H., M.I.F., E.C.B., C.V.-F., V.L., R.M., and R.R.); ICREA 107. R. Rodríguez-Varela, T. Günther, M. Krzewińska, J. Storå, T. H. Gillingwater, M. MacCallum, Academia program (R.R.); John Templeton Foundation grant 61220 (D.R.); and Paul Allen J. L. Arsuaga, K. Dobney, C. Valdiosera, M. Jakobsson, A. Götherström, L. Girdland-Flink, Family Foundation (D.R.). D.R. is an Investigator of the Howard Hughes Medical Institute. Genomic analyses of pre-European conquest human remains from the Canary Islands Author contributions: V.V.-M., R.R., and W.H. conceived the study. C.O., C.R.-H., R.M., M.I.F., reveal close affinity to modern North Africans. Curr. Biol. 28, 1677–1679 (2018). E.C.B., C.V.-F., K.W.A., D.C.S.-G., G.G.A., M.S.H.P., M.P.d.M.I., A.R., V.B., J.P., A.M., J.L., J.S., A.P.M., Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 18 of 19 SCIENCE ADVANCES | RESEARCH ARTICLE J.K., L.G.S., V.L., and R.R. provided archaeological material and advised on the archaeological paper and/or the Supplementary Materials. Genomic data (BAM format) are available at the background and interpretation. F.A., M.H., W.H., C.F., V.V.-M., and C.O. performed the European Nucleotide Archive under accession number PRJEB46907. laboratory work. V.V.-M., W.H., A.B.R., and A.C. performed data analysis, with J.K. providing feedback in the analysis and W.H. providing guidance and methodologies. V.V.-M., R.R., Submitted 26 March 2021 W.H., and A.B.R. wrote the manuscript with input from all coauthors. Competing interests: Accepted 28 September 2021 The authors declare that they have no competing interests. Data and materials Published 17 November 2021 availability: All data needed to evaluate the conclusions in the paper are present in the 10.1126/sciadv.abi7038 Villalba-Mouco et al., Sci. Adv. 7, eabi7038 (2021) 17 November 2021 19 of 19