Genomic transformation and social organization during the Copper Age–Bronze Age transition in southern Iberia-2
Mediterranean and central European ancestries shaped the genetic profile of southeastern BA groups in Iberia
To explore the genetic turnover and the contribution of the local groups to the newly formed BA genetic profile in Iberia, we systematically tested a series of qpAdm models. We started by using the distal ancestry sources Anatolia_N, WHG, GoyetQ2, Yamnaya_Samara, and Iran_N to model the genetic ancestry components of Iberian BA groups (table S2.10 and fig. S6). We found that the local traces of GoyetQ2, a characteristic but variable component of southern Iberia CA individuals, were no longer detectable, suggesting a dissolution of geographic substructure in BA Iberia with respect to HG ancestry. We explain this by the spread of steppe-related ancestry from North to South (7) that also contributed northern and central Iberian ancestry to the South, diluting the subtle GoyetQ2 signal to a level beyond the limits of detection (text S8). By using the same qpAdm model, we also observed that Almoloya_Argar_Early, Almoloya_Argar_Late, SE_CabezoRedondo_BA, and Bastida_Argar cannot be modeled with Yamnaya_Samara as a single source but find better support with a combination of Iran_N and Yamnaya_Samara, however, without reaching P values ≥0.05 in Almoloya_Argar_Early and Late and SE_CabezoRedondo_BA (table S2.10 and fig. S6). These El Argar groups (Almoloya and Bastida) are also slightly shifted to the right on the PC1 axis, in the direction of Mediterranean BA groups with excess Iran_N-like ancestry, such as “Minoans,” who only carry Iran_N-like ancestry but not steppe-related ancestry, or “Mycenaeans,” who carry a mix of both (71), and that has also been shown for some BA individuals from Sicily_MBA (51) and for Sardinians here (Fig. 3A).
To explore the reasons for the model rejection observed in southeastern BA groups, we tested several qpAdm models (text S8). We used two- and three-way competitive qpAdm models to test whether a third source with Iran_N-like ancestry was needed. In the two-way model (table S2.11), Germany_Bell_Beaker (the largest number of individuals representing BB) was used as a fixed proximal ancestry source together with a local CA source (either N_Iberia_CA, C_Iberia_CA, or SE_Iberia_CA). For target groups in which these two sources were rejected, we iteratively tested central and eastern Mediterranean populations as a third source, as some of them are known to carry Iran_N-like ancestry (Fig. 4A and table S2.11). Using C_Iberia_CA as a local source of ancestry and Germany_Bell_Beaker as a proxy for steppe ancestry, only models with Almoloya_Argar_Early and Late as a target were rejected, suggesting that a third component was required for these Iberian BA groups. We found that adding any Mediterranean population as a third source improves the model fit, albeit without reaching P values >0.05. However, we obtained P values >0.05 when using Iran_Chalcolithic as a third source, which suggests that a higher proportion of Iran_N-like ancestry is needed (text S8 and table S2.11). Notably, when we exchanged C_Iberia_CA with the local SE_Iberia_CA, the model did not find support for a larger number of Iberian BA groups (N_Iberia_BA, SE_Iberia_BA, Almoloya_Argar_Early and Late, and C_Iberia_CA_Stp), and adding a third source in this constellation did not improve the model fit either (Fig. 4A and table S2.11).
FIG. 4. Modeling genetic ancestry in Iberian BA groups.(A) Modeling Iberian BA individuals with steppe-related ancestry as a two- and three-way qpAdm admixture model using proximal sources C_Iberia_CA/SE_Iberia_CA, Germany_Bell_Beaker, and Iran_C (table S2.11). (B) Modeling Iberian BA individuals with steppe-related ancestry as a two- and three-way qpAdm admixture model using proximal sources C_Iberia_CA_Stp, C_Iberia_CA/SE_Iberia_CA, and Iran_N (table S2.1). P values for each group are given inside each column. Faded colors indicate rejected models when applying a P value cutoff of ≤0.05.
To rule out that the rejection of models involving Germany_Bell_Beaker + C_Iberia_CA in La Almoloya (Early and Late) is due to specific north/central European LN admixture in Germany_Bell_Beaker that was not present in the source of steppe-related ancestry in Iberia, we fixed the source of direct steppe-related ancestry to Yamnaya_Samara, as we expect the same distal steppe source contributing to all Iberian groups, and iterated through the proximal local CA sources, which in these models contribute 100% of the Chalcolithic (i.e., Neolithic farmer) ancestry. We added distal Iran_N (instead of Iran_Chalcolithic, which could contribute other confounding CA ancestries) as a third source only when needed (table S2.12 and fig. S7). By using C_Iberia_CA and Yamnaya_Samara in a two-way model, we could successfully model all Iberian BA groups, but again to the exclusion of the southeastern Almoloya_Argar_Early and Late, and SE_CabezoRedondo_BA groups. Notably, Bastida_Argar also failed for the distal model 2 (table S2.10). However, these three groups returned values ≥0.05 in the proximal local CA substrate model when Iran_N was added as a third source (fig. S7 and table S2.12). In turn, when we exchanged C_Iberia_CA with a local SE_Iberian_CA source, most of the models were rejected, which confirms the genetic substructure in Iberia_CA groups. Using SE_Iberia_CA as substrate only improves the model fit for southeastern Iberian groups, suggesting that the first individuals who brought steppe-related ancestry to Iberia admixed locally with the different CA groups (table S2.12). Notably, SE_Iberia_CA is a sufficiently supported local proxy for Balearic CA_Stp, BA, and LBA groups. We interpret this signal with caution, as Mallorca_CA_Stp and Menorca_LBA are only represented by a single individual each; however, this signal is well supported in the newly typed Aritges_LBA group from Menorca, which consists of six individuals (fig. S7 and table S2.12).
We then tried to further specify the proximal sources without affecting the power of resolution to distinguish between them. We exchanged Yamnaya_Samara by C_Iberia_CA_Stp as the most proximal local source with steppe-related ancestry (6) and combined these with the different Iberian CA sources (table S2.13 and text S8) in two-way models, only adding Iran_N as a third source when needed. Using increasingly more possible proximal sources for steppe-related ancestry, we could replicate the results of the distal models. However, we can show that Almoloya_Argar_Early and Late, Bastida_Argar, SE_CabezoRedondo_BA, and Aritges_LBA still require Iran_N as a third source in all models with C_Iberia_CA and C_Iberia_CA_Stp. In contrast, Iran_N is not required when SE_Iberia_CA is used as a local substrate, although models for Almoloya_Argar_Early and Late are still rejected at P ≤ 0.05. This finding is best explained by either an unresolved distinctive ancestry of La Almoloya individuals or, alternatively, the ability to detect such subtle signals only in groups consisting of a larger number of individuals, as is the case for Almoloya_Argar_Early (n = 36) and Almoloya_Argar_Late (n = 22), which provide the statistical power to reject simpler models. When SE_Iberia_CA is used instead of C_Iberia_CA to model La Bastida, the model is no longer rejected, which suggests a direct contribution of local CA groups to those of the southeastern Iberian BA. Together, the results hint at an additional, albeit subtle, Iran_N-enriched ancestry that is present in SE_Iberia_CA groups and thus potentially predating the El Argar BA.
We also tested whether Almoloya_Argar_Early and Late can be modeled when Iran_N is exchanged by geographically more proximal Iran_N-rich sources involving central and eastern Mediterranean populations (table S2.14 and text S8) and by moving Iran_N to the outgroups. However, this did not result in any supported models (text S8). While we were not able to find a proximal Mediterranean source, we note that adding a central Mediterranean population to the outgroups (Sicily_EBA, Greece_EBA, or Greece_MBA) decreases the model support (P values) for Almoloya_Argar_Early and Almoloya_Argar_Late, indirectly attesting to the importance of central Mediterranean BA. However, we observed no changes in P values by adding other Mediterranean groups to the outgroups (text S8 and table S2.15).
The complete turnover on the Y chromosome to R1b-Z195 (a lineage derived from P312), observed in all males at La Almoloya (29 male individuals tested) and La Bastida (except for one E1b lineage dated around 2134 to 1947 cal BCE; 7 males tested), is another independent source of evidence of gene flow predominantly during the time of the CA-BA transition (table S1.1 and text S3). Notably, Y lineage R1b-Z195, the most common Y lineage in BA Iberia, ultimately derives from a common ancestor R1b-P312 in central Europe but already differs from other derived Bell Beaker R1b variants reported from central Europe and the British Isles. R1b-Z195 has been found in Sicily_CA_Stp (previously assigned to the EBA and thus considered outliers) and Sicily_EBA (51). However, the genealogical and geographic link to other R1b variants still remains unclear. The subtle presence of Iran_N-like ancestry in El Argar and Y haplogroup R1b-Z195 in Sicily opens the possibility of gene flow not only from Iberia to Sicily as previously suggested (51) but also in the opposite direction, implying reciprocal contact with the western and central Mediterranean during the BA, although direct contacts are hardly noticeable in the archaeological record.
Together, these results suggest a dual genetic contribution to the formation of the BA genetic profile of southeastern Iberians in addition to a local CA genetic substrate. The major additional ancestry source resembles central European Bell Beaker groups, which first contributed ancestry to northern Iberia, followed by a southward spread in the form of C_Iberia_CA_Stp. A second minor ancestry component is an Iran_N-rich/central Mediterranean source, which is restricted to individuals from BA El Argar contexts. The timing of the last contributing component remains unclear and points to either a Neolithic legacy that persisted throughout the local CA or a subtle trace of connections to insular central Mediterranean BA groups.
A late Argar genetic outlier makes links to North Africa and the central Mediterranean
The PCA also revealed a genetic outlier among the newly typed individuals (Fig. 3A). The male individual ZAP002 from the late Argaric site of Lorca-Zapatería falls outside the Iberian_BA cluster on the PCA plot, with a clear shift toward the central Mediterranean. We formally tested for a specific attraction to groups from the eastern Mediterranean, the Near East, and Africa with an f4-admixture test of the form f4(ZAP002, Almoloya_Argar_Early; eastern_Med/Africa, Chimp) (Fig. 5A, table S2.16, and text S9). The resulting f4-values were consistent with zero, showing that ZAP002 and Almoloya_Argar_Early (the Argar group with the largest number of individuals) are symmetrically related to eastern Mediterranean and African groups. However, we noticed a near-significant positive f4-value with a z score = 2.4 for Morocco Iberomaurusian. The attraction to Mbuti in an analogous f4-test suggests either African ancestry or high levels of Basal Eurasian ancestry shared with Morocco Iberomaurusian and Natufians (table S2.16). We also confirmed a negative deviation from zero in f4(test, EHG; CHG, Mbuti), indicative of either African or Basal Eurasian ancestry. Here, we assume that the split of EHG and CHG groups [Caucasus HG; a basal form of HG ancestry south of the Caucasus (72)] is sufficiently deep to result in negative values for any test population. Following this rationale, the finding of negative values suggests a high amount of shared deep ancestry with Mbuti, which also explains the different f4-value for ZAP002 in Fig. 3B when compared to other El Argar individuals. We observe negative values for ZAP002 (z score = −2.14), as well as Morocco Iberomaurusian and Natufian individuals, who provide additional support for this claim (Fig. 5B, table S2.17, and text S9).
FIG. 5. Summary of genetic analyses that explore the genetic outlier status of individual ZAP002.(A) f4-statistics to detect Basal Eurasian ancestry and/or African ancestry. Cladality test with ZAP002 and individuals from the site La Almoloya (ALM) and a set of test populations including Africans and Near Eastern (NE) groups, and chimpanzee in the outgroup. The results show that only Moroccan, Iberomaurusian, and Mbuti show a shift toward positive f4-values (table S2.16). (B) f4-statistics to detect either high Basal Eurasian or North African ancestry in which negative values imply excess affinity to Mbuti (table S2.17). (C) Supported and nonsupported qpAdm admixture models for outlier individual ZAP002. On the left, C_Iberia_CA is located in the outgroups and the model finds support with Sicily EBA as a source, pointing to a nonlocal origin of individual ZAP002 (table S2.18).
Using qpAdm, we further explored candidate sources of distinct ancestries such as Iberian CA and Germany_Bell_Beaker and a candidate list of central Mediterranean sources from Sardinia and Sicily. Using the same set of outgroups as for the main southeastern Iberian BA groups, we did not find a statistically well-supported model for ZAP002. However, moving Morocco_Iberomaurusian from the outgroups to the sources resulted in a model fit (≥0.05) when a central Mediterranean population was included as a source. We found that ZAP002 does not require a local Iberia_CA source and can be modeled as a three-way mixture of Germany_Bell_Beaker, Morocco_Iberomaurusian, and either Sicily EBA, Sardinia EBA, or Sardinia Nuragic BA (Fig. 5C and table S2.18). These qpAdm results are congruent with the location of ZAP002 in PC space, where he clusters with Sicily EBA individuals (Fig. 3A). Rotating the Iberian CA source to the outgroups in the same model still upholds the model fit (Fig. 5C, table S2.18, and text S9), suggesting an entirely nonlocal origin of individual ZAP002 despite being granted the same burial treatment as other individuals from El Argar (pithos burial and Argaric pottery) (Fig. 5C).
Insights into phenotypic variation, demography, and social correlates of CA and EBA El Argar societies
We explored the phenotypic variants included in the 1240k SNP panel among all well-preserved Iberian CA and BA individuals (>400,000 SNPs covered). In addition to 14 variants associated with the ability to digest lactose and alcohol, adaptation to metabolize fatty acids, predisposition to celiac disease, and resistance to infectious diseases, we also investigated 41 SNPs encoding for skin/hair and eye pigmentation (figs. S8 and S9, tables S2.19 and table S2.20, and text S10). We found no significant increase or decrease in allele frequency of any derived allele between the CA and BA (fig. S9 and table S2.20).
To gain insight into the effective population size during the CA and BA, we applied hapROH to identify runs of homozygosity (73, 74). We found no notable differences in short ROH (<2 cM) patterns between the two chronocultural periods, indicating that all of the individuals under study were drawn from populations with comparably and sufficiently large effective population sizes. We also found no evidence for a decline of CA population size before the emerging BA. In both groups, we did not detect signs of inbreeding, represented by long ROH segments (20 to 300 cM) (Fig. 6A and table S2.21).
FIG. 6. Biological kinship and female-male dynamics at the site La Almoloya.(A) Sum of all ROH segments measured with hapROH for CA and BA above 400,000 SNPs in 1240k panel (table S2.21). (B) qpAdm z scores between autosomes and the X chromosome showing no signal for male bias related to steppe ancestry in the model with distal sources (right) or proximal sources (left) (table S2.22). (C) f3-outgroup statistics of the form f3(female, female; Mbuti), f3(female, male; Mbuti), and f3(male, male; Mbuti), highlighting a closer relationship among males than among females or closer relationships whenever they involve at least one male individual (first- and second-degree related pairs are excluded from the calculation) (table S2.23). (D) PWMR values of individuals compared to the PWMR average of the site La Almoloya. Dashed lines indicate a lower PWMR average and, thus, higher pairwise relatedness in adult males than adult females (table S2.24).
To assess a potential sex bias in the steppe-related ancestry contribution as postulated previously (7), we applied distal and proximal qpAdm models (text S8) to the X chromosome and autosomes (Fig. 6B and table S2.22). As males contribute, on average, only one-third of the X chromosomes to the next generation (75), a lower proportion of any ancestral component on the X chromosome than on the autosomes would thus be indicative of male bias concerning the respective component (76, 77). In turn, if the ancestral component is statistically higher on the X chromosome, this indicates a female bias. On the basis of this rationale, we do not observe significant male bias in steppe-related ancestry using either distal or proximal sources (Fig. 6B and table S2.22). The fact that the male bias is not detectable could be indicative of an already balanced ancestral component in both sexes, as is reflected in the work of Mittnik et al. (77), where the male bias in the steppe component is only detected in Corded Ware, but no longer in Bell Beaker or BA populations. We also took advantage of the recent extensive excavations and sampled all individuals with suitable morphological preservation (n = 86) from the site La Almoloya, resulting in high-quality genome-wide data for 67 individuals that allowed insight into an El Argar community. We estimated the genetic relatedness using a combination of the three methods: pairwise mismatch rate (PWMR) (48), READ (Relationship Estimation from Ancient DNA) (49), and lcMLkin (47) for pairs of individuals with more than 1000 shared SNPs (see Materials and Methods, fig. S2, and text S5). All reconstructed pedigrees will be published in an accompanying paper, fully integrated into the anthropological and archaeological context. Thus, we report here the main observations of biological relatedness at the metalevel. We found no first-degree relationships between adult women (0 of 30 adult women analyzed). All first-degree relationships among adults involved at least one adult male (three first-degree relationships among adult men and women). We also found no second-degree relationships between adult women. The few second-degree relationships involving adults were all between males (4 of 18 adult males analyzed).
We then formally tested whether the males of La Almoloya had more close relatives at the site than the females. We calculated f3-output statistics after removing first- and second-degree pairs of the form f3(female, female; Mbuti) and f3(male, male; Mbuti) and found that despite a similar average value in female:female and male:male comparisons (Fig. 6C and table S2.23), males have a more skewed f3-distribution than females, which suggests a more close genetic relationships (e.g., at third to fourth degree) among males with other males (P = 1.9 × 10−18) and males with other females (P = 2.2 × 10−4) than for females with other females (see Materials and Methods). We specifically tested for signals that suggest patrilocality using the PWMR extracted from READ (P0_Normalized) (49) and built a matrix to calculate the biological relatedness of each adult individual to the entire group (Fig. 6D and table S2.24). We observe a tendency that males are, on average, more closely related to other individuals from the site than females, despite there also being some adult females with relatives at the site, which renders the formal test for patrilocality nonsignificant (see Materials and Methods). At the intrasite level, this result could also suggest a founder effect from the male line with a higher reproductive rate.
DISCUSSION
Our focused study on the southeastern part of Iberia highlights the potential of genetic analysis on a regional level in the light of broader-scale population dynamics in Europe. Overall, we were able to show that the structure of CA populations shows persistent genetic stability and continuity since the Neolithic. This is coupled with a differential HG ancestry in the South in addition to potential early contact between southeastern Iberian CA and Mediterranean populations, which is also revealed by numerous threads of archaeological evidence (13–17). We provide nuanced evidence for the arrival of steppe-related ancestry in southeastern Iberia by 2200 cal BCE, when widespread social and cultural changes occurred across much of southern Iberia, including the end of the ditched enclosures system of population aggregation and the emergence of BA groups. In the Southeast, all sampled BA groups can be shown to carry steppe-related ancestry, while the more extensively sampled El Argar groups also carry excess Iran_N-like ancestry, which has also been observed in other CA and BA Mediterranean groups. The analysis of one outlier from El Argar reveals a central Mediterranean migrant individual with additional North African ancestry. Overall, we propose that El Argar has likely formed from a mixture of new groups arriving from north-central Iberia, which already carried central European steppe-related ancestry (and the predominant Y-chromosome lineage) and local southeastern Iberian CA groups that differed from other regions in Iberian in that they carried excess Iran_N-like ancestry similar to eastern and/or central Mediterranean groups. Moreover, the large sample sizes of early and late El Argar groups show that the Iran_N-like ancestry contribution in the local CA stratum was not enough to explain this type of ancestry satisfyingly, which, in turn, argues for a continued connection to and genetic influence from the Mediterranean BA at least until the end of the El Argar period.
Last, we were able to shed light onto the population structure of Early Bronze El Argar societies at the intrasite level. In La Almoloya, a closer genetic relationship among males is a strong indicator of patrilineality.
MATERIALS AND METHODS
Experimental design
Archaeological samples
We processed 244 individuals from the Iberian Peninsula and the Balearic Islands dated to a time span from the CA (3000 cal BCE) to the LBA (1200 cal BCE). A detailed description of the archaeological context, sites, and individuals is reported in text S1 and summarized in table S1.1.
Laboratory work
All the laboratory work was performed in the dedicated aDNA facilities of the Department of Archaeogenetics of the Max Planck Institute for the Science of Human History in Jena, Germany.
Sampling
We preferentially sampled petrous bones and teeth (essentially molars). Details of skeletal elements samples are listed in table S1.1. Before sampling, all samples were irradiated with ultraviolet light on either side for 30 min each, after which the surface was cleaned with low-concentration bleach solution (3%). For the petrous bones, we either drilled from the outside, removing the outer layer, or cut the petrous pyramid longitudinally to drill the dense part directly from either side. Molars were cut below the enamel-dentine junction, and tooth powder was drilled from the pulp chamber. We collected between 50 and 100 mg of bone or tooth powder per sample for DNA extractions.
Extraction
DNA was extracted following the modified version of (78) described in www.protocols.io/view/ancient-dna-extraction-from-skeletal-material-baksicwe.
https://www.science.org/doi/10.1126/sciadv.abi703[/size]8