A 23,000-year-old southern Iberian individual links human groups that lived in Western Europe before and after the Last Glacial Maximum

Human populations underwent range contractions during the Last Glacial Maximum (LGM) which had lasting and dramatic effects on their genetic variation. The genetic ancestry of individuals associated with the post-LGM Magdalenian technocomplex has been interpreted as being derived from groups associated with the pre-LGM Aurignacian. However, both these ancestries differ from that of central European individuals associated with the chronologically intermediate Gravettian. Thus, the genomic transition from pre- to post-LGM remains unclear also in western Europe, where we lack genomic data associated with the intermediate Solutrean, which spans the height of the LGM. Here we present genome-wide data from sites in Andalusia in southern Spain, including from a Solutrean-associated individual from Cueva del Malalmuerzo, directly dated to ~23,000 cal yr bp. The Malalmuerzo individual carried genetic ancestry that directly connects earlier Aurignacian-associated individuals with post-LGM Magdalenian-associated ancestry in western Europe. This scenario differs from Italy, where individuals associated with the transition from pre- and post-LGM carry different genetic ancestries. This suggests different dynamics in the proposed southern refugia of Ice Age Europe and posits Iberia as a potential refugium for western European pre-LGM ancestry. More, individuals from Cueva Ardales, which were thought to be of Palaeolithic origin, date younger than expected and, together with individuals from the Andalusian sites Caserones and Aguilillas, fall within the genetic variation of the Neolithic, Chalcolithic and Bronze Age individuals from southern Iberia.

various chambers, of which the Lower Galleries stand out: Saco Room, Star Room, Calvary Gallery and Labyrinths. Another area of the cave is elevated and known as High Galleries. Inside and near the entrance, there is a heavy/large cone of sediments, which forms part of the stratification and which has been used as a stairway access area since the mid-19th century.
Since its discovery in 1821, the cave system has a long research history with tourist exploitation by Doña Trinidad Grund followed by the integration into prehistoric sciences through the first studies carried out by Henri Breuil and Miguel Such from 1918, who documented the first panels with Paleolithic art 4 . After a long period of abandonment and after the Spanish Civil War, the cave was re-opened in the mid-1980s and various projects were carried out to clean up and document the archaeological record and rock art. It was thus possible to document more than 1,000 painted and engraved motifs, with animals, female human figures, hands, lines, points, and stains. We refer to Cantalejo and colleagues 5,6 for a historiographic synthesis and to numerous studies on the prehistoric art [6][7][8][9][10][11][12][13][14][15][16][17] .
In recent years a General Research Project (GRP) was developed by the authors of this study, with authorization from the Andalusian Government, which allowed the examination and excavation of four areas. Archaeological surveys have obtained a significant chronostratigraphic sequence from the Middle Paleolithic to recent Prehistory (Neolithic and Chalcolithic). Ongoing work by the GRP team includes further archaeological surveys, a variety of chronostratigraphic and interdisciplinary studies on the analysis of archaeological products, and also involves new technologies to study Paleolithic rock art and anthropological remains [18][19][20][21][22][23][24] .
We sampled individuals from occupation layers in zone 2 at Cueva de Ardales, which were directly dated to the Neolithic and Chalcolithic, and thus surprisingly younger than the stratigraphy had indicated. Ongoing excavations revealed that they were intrusive in a Solutrean layer due to disturbance by an animal burrow 25

Necrópolis de las Aguilillas
Las Aguilillas is a necropolis located between the towns of Ardales and Campillos Aguilillas [26][27][28] . The necropolis is artificially carved into an Upper Miocene sandstone outcrop, where seven tombs have been documented. These tombs correspond to oval-shaped chambers, which sometimes have niches and entrance corridors, with dimensions between 3.2m and 2.6m and horizontal domes between 1.9m and 1.6m.
One tomb, specifically structure 6, is a megalithic construction formed by an access corridor 7.5m long, a first rectangular chamber and a small corridor through which a second chamber of smaller proportions is accessed. The two chambers have a flat roof made of large horizontal sandstone slabs. The grave goods are very characteristic, with bowls, pots and ceramic vessels of good quality. The lithic industry includes foliaceous points of bifacial carving and flat retouches with a concave base. There are pressure carving blades and sickle elements. An important set of possible peaks for the elaboration of artificial sandstone caves -more than 200 specimens -were documented in a niche. Among the metallic elements there are palmela points and punches 26,27 .
This type of necropolis with collective burials is very characteristic of the central zone of Andalusia (the interior of Málaga, the plains of Antequera, and the current countryside of Seville, Córdoba and Jaén) with very defined chronologies between the end of the 4 th and the beginning of the 2 nd millennia cal. BCE. However, new direct radiocarbon dates from the individual analyzed from cave 2 have revealed an older date, suggesting the beginning of funerary use during the Neolithic.
It is interesting to note that in front of the necropolis is the settlement of El Castillón, located about 500m to the east, on the other side of the Guadalteba River, which is considered the town and habitat to which the necropolis is associated 29

Necrópolis de los Caserones
El Caserón is a necropolis located in the town of Ardales (Málaga), at the base level of a dam next to the Turón River Turón. The site has single and double burials, made up of a box of slabs, or orthostats, in a quadrangular or rectangular shape with cover slabs, which correspond to the Bronze Age 29,30 . These tombs contain burnished ceramic vessels, with faired shapes, together with metallic elements such as daggers with arsenic copper rivets, swords, quadrangular awls, and silver rings 31 . Culturally, these tombs mark a process of important social hierarchization, with a gradual abandonment of collective burials of Chalcolithic tradition followed by burials that emphasise the individual and smaller families 29,30 .
It is interesting to note that these cist necropoles date after the collective burials in the Guadalteba region, such as the one mentioned in Aguilillas. Los Caserones, located next to the Guadalteba and Turón rivers (Las Grajeras, Morenito, Raja del Boquerón, Caserones, Retamar, La Bolina) 30,32 , may correspond to the westernmost limit of the El Argar Culture in southern Iberia [33][34][35] . However, the 14

Quality controls
A combination of several quality controls were applied at library and/or individual level to evaluate the authenticity of the aDNA sequences.

PMD filtering
We used PMDtools 37 to identify reads with high rates of DNA damage (PDM score > 3) as these are less likely to come from modern DNA contamination. Supplementary Although there was no indication of contamination using other approaches (see

Sex determination
We used sex determination to detect contamination from the opposite sex 38 . Here we calculated the X-ratio and Y-ratio (see methods) in all individual libraries and merged libraries before and after PMD filtering (Supplementary Fig. 1c After PMD filtering, the X-ratio and Y-ratio remains similar for all the individuals ( Supplementary Fig. 1c).
Supplementary Fig. 1 I Summary of quality controls. a,Damage rate per single-stranded library and individual at 5' end. b, Sex determination using X/autosomes-ratio vs Y/autosomes ratio as indicator of genetic contamination from the opposite sex performed at library level. c, Sex determination plot using the X/autosomes-ratio versus Y/autosomes ratio as indicator of genetic contamination from the opposite sex. We show the merged libraries from the same sample and their comparison with the PMD damage-filtered version of each merged library.

Pairwise mismatch rate (PMR) to identify duplicates and contamination among libraries
Due to the fact that MLZ003 and MLZ005 showed the same uniparental markers (MT and Y-chromosome haplogroups, (Supplementary  Fig. 2b). A graphical summary of the baselines and PMR values is shown in Supplementary Fig. 2a-2b, before and after the exclusion of the contaminated libraries from MLZ003, respectively.

X-contamination estimation
We used ANGSD (method 2) to estimate heterozygosity at polymorphic sites on the X chromosome in males 41 . We used the threshold of 100 X-SNPs covered twice to give a confident contamination estimation. As many of our libraries did not reach that threshold, we estimate contamination in the merged bam file which contains all libraries from the same individual. Following this criteria, only the Paleolithic individual MLZ005 reached 115 X-SNPs covered twice and gave an X-contamination estimation of 0.041 ±0.041. MLZ003 only reached 56 X-SNPs covered twice, respectively, and thus x-contamination estimation was not possible. The nuclear contamination rate of the merged MLZ005003 was estimated at 3.7%. In the rest of the male individuals, X-contamination was measurable and below 0.032 ±0.025 (Supplementary Table 1.8).

Mitochondrial DNA-contamination estimation
We merged several mitogenome-captured libraries (Supplementary Table 1

Comparison of standard and PMD-filtered genotype data, and data restricted to transversions (TV).
For the UP individuals, we performed several PCAs to evaluate the deviations in PC-space between standard and PMD-filtered genotype data of MLZ and used these as an indirect quality control, assuming that a highly contaminated individual should plot far from the PMD-filtered version. Additionally, we also restricted the genotyping to transversions only (MLZ005003.TV), and plotted all versions to evaluate whether DNA damage had an effect on the position in PC space. From the PCA plots we consistently observed that MLZ005003 does not deviate from MLZ005003.PMD or .TV (Supplementary Fig. 3). We thus retained individual MLZ005003, but we nonetheless provide separate results for the three versions in all statistics presented in this study.

Relationship of MLZ to other ancient individuals
First, we used f 4 -statistics of the form f 4 Fig. 5b). This could also be the case for some WHG, Scandinavian HGs (SHG) and EHG individuals (e.g. Serbia and Romania_Iron_Gates HG, Motala_HG, Russia_EHG, and Latvia HG groups

Traces of Early Asian ancestry
The 40 ka cal BP Tianyuan individual from East Asia was shown to share more ancestry with present-day East Asians and Native Americans than with present-day or ancient Europeans, which led to the conclusion that the split between ancient Asians and Europeans must have happened at least 40 ka ago 47 . Following these findings, the ancestry found in Tianyuan was defined as Early Asian ancestry.
Intriguingly, among Paleolithic individuals published by Fu et al. 43 , an Aurignacian-associated individual from Goyet cave (in today's Belgium) named Goyet Q116-1 was found to be closely related to ancient Europeans but also shared more alleles with Tianyuan than any other individual in Europe. Yang et al. 48  Magdalenians share most of their ancestry (Fig. 2b) Europe as the observed f 3 -outgroup value for Moita do Sebastião is even higher than those of EHGs, who carry more ANE ancestry and also ancestry from Asian populations 49 , and thus also show elevated positive f 3 -statistics (lighter colors in Fig.   5b). This finding highlights the genetic legacy of the first early Eurasians in southern Europe, and a genetic connection with the Tianyuan individual, which can be traced in European individuals for more than 30,000 years.

Exploration of Basal Eurasian and Near Eastern ancestry
The results of the f 4 -statistics shown in Fig. 4a  Of note, this f 4 -statistic could be negative due to an attraction of Natufians to Mbuti, and driven by the Basal Eurasian ancestry present in Natufians. The f 4 -statistic above is also negative when Zlatý kůň is used as an outgroup (Supplementary  Fig. 6a-6b

Phylogenetic position of MLZ (qpGraph modeling)
We used qpGraph modeling to establish the phylogenetic position of MLZ. We ( Supplementary Fig. 8a). The decision to add a Villabruna-like lineage next was and also added a theoretical node called Basal_WHG.
We continued with the tree represented by Supplementary Fig. 8a by adding MLZ as a result of the admixture between an ancestral node of the Villabruna clade and the sister branch of Goyet Q116-1. This resulted in a worst Z-score of -1.154 ( Supplementary Fig. 8b) We continued adding the individual from El Mirón as a representative for the Magdalenian-associated ancestry in the Iberian Peninsula as well as Kostenki14. El Mirón could be fitted as an admixture of Villabruna and MLZ lineages, without any additional contribution with a Z-score = 1.424 (Supplementary Fig. 8c). Finally, we added Vestonice16, representing the central European Gravettian as a mixture of Kostenki14 clade and Villabruna-like ancestry (Supplementary Fig. 9).   Fig. 10a, Supplementary Table 2.16 Sub-Saharan ancestry, which also contributed to Morocco Iberomaurusian and thus would be a confirmation of an ancestry contribution from Morocco Iberomaurusian.

Testing for North African ancestry in MLZ and
However, significantly negative f 4 -statistics for all test populations rule out a contribution of Sub-Saharan ancestry, and thus Morocco Iberomaurusian, or at least at much lower levels than in Natufians, which were used as baseline population in the test (Supplementary Fig. 10b

Neolithic individuals
In PC space, we noticed a slightly shifted position of Southern_Iberia_EN towards negative values of PC1 and positive values of PC2 (Supplementary Fig. 11a).  62 , and characterized by a subtle proportion of Iran_N-like ancestry 63 .
To formally test these hypotheses we used f 4 -statistics and qpWave/qpAdm following strategies described below, in which we iterated through various sources and rotated these into outgroups.  Regarding qpAdm, several approaches have been used. The classical approach iterates through sources. and retains the models which return a good fit to model the target group (qpAdm p-value > 0.05). Using this approach we look for nested models,

Absence of North African ancestry in Southern Iberia EN
i.e. most parsimonious, fitting models with minimum number of sources needed. The alternative approach is known as rotating strategy 64 . The aim of this approach is to distinguish the true sources needed in the model by moving them to the list of outgroups. If, by doing, so the associated p-value becomes <0.05, the model is rejected (assumptions violated), which in turn indicates that the rotated population is better suited as a direct source. This strategy is helpful when little is known about the populations that are tested.
Here, we address a specific archeological hypothesis that can be translated into differential sources of ancestry, for which we decided to use the source iteration approach and the rotating strategy only when the two sources tested in analogous models cannot be distinguished. We chose a set of basic outgroups, including African and distal HG ancestral sources which we deemed representative of the main ancestral components in our target groups: All results are reported in Supplementary Table 2. 19 and Supplementary Fig. 11d.

B) Rotating strategy to distinguish the best source (Iran_N or MLZ)
We first rotated MLZ from the sources to the outgroups to model When MLZ was added as an outgroup, we observed that none of the previous models were supported, including those for Iberian_EN_North. This suggests that MLZ ancestry is critical in both target groups (northern and southern). In contrast, when Iran_N was added to the outgroups and MLZ used as a source, we still obtained a good fitting value for the three-way model which includes MLZ as a third source, both in Iberian_EN_North and _South. We also observe that the estimated proportion of MLZ ancestry is bigger in Iberia_EN_South ((12.5% ± 3) of the total ancestry profile) than in Iberian_EN_North ((6.5% ± 1.9) of the total ancestry profile). The smaller ancestral proportion of MLZ-like ancestry found in Iberian_EN_North also explains why the nested two-way model including only Anatolia_Neolithic and WHG was sufficient to obtain a good initial model fit with qpAdm in the iteration approach. Taken together, we found that the total amount of HG-like ancestry (WHG + MLZ) is larger in the South.

C) Iteration through sources to evaluate different Solutrean, Magdalenian or Mesolithic ancestries.
We aimed to distinguish between Solutrean, Magdalenian or Southern Iberian Mesolithic ancestries that contributed to the HG-like ancestry carried by the Iberian EN groups. As the three HG sources are very similar, we iterated through the three potential sources of ancestry as described below: By applying this model, we obtained a good model fit value using any of MLZ, El Mirón or Moita do Sebastião. We interpreted this result as having reached the limit of the resolution of the data due to the similarity of the three sources.
Chronologically it is clear that the only population who could have been in direct contact with the EN farmers was the Mesolithic population, but the fact that all three models worked well highlights the higher HG population continuity in Iberia, particularly in the South.

D) Rotating strategy to distinguish the best source (MLZ, El Mirón, Moita do Sebastião)
Finally, in order to gain more power of resolution to distinguish between Soutrean, Magdalenian or Southern Iberian Mesolithic ancestries as potential sources, we rotated each of them to the outgroups while keeping the others as sources. We obtained a good model fit for every source-outgroup combination. The fact that all three models were well supported highlights the higher proportion of HG contribution, and thus a continuation of this ancestry in Iberia, and in particular in the South.

Exploration of individual ADS007.
The presence of 'steppe-related ancestry' was confirmed for individuals ADS007 by