Ancient genomes from present-day France unveil 7,000 years of its demographic history

The French Mesolithic Substrate.

We generated genomic data for five Mesolithic individuals from Les Perrats cave (Agris, Charente) dated from 7177 to 7057 calibrated years BCE (calBCE) including low-coverage shotgun sequencing for three of them (0.111 to 0.134×). All possess haplogroup U5b, demonstrating in the territory of present-day France the presence of a Mesolithic substrate akin to that described elsewhere, where it was characterized by haplogroups U5b (85%) and, to a lesser extent, U5a (15%) (3, 23).

Mesolithic hunter-gatherers (HGs) discovered in France are located near the edge of a principal-component analysis (PCA) plot obtained from genome-wide genotypes, falling near the Mesolithic individuals from western Europe (5, 6, 9) (Fig. 1B and SI Appendix, Fig. S4-3). While the ancestry from an ∼14,000-y-old individual from Villabruna in Italy was found to be dominant in Holocene HGs across western and central Europe (3), recently published studies conducted on ancient Iberians highlighted the survival of Magdalenian-associated ancestry (11). This dual ancestry is found in all Iberian HGs (ranging from 23.7 to 75.3%), while in most other regions only the Villabruna-associated ancestry remained. In order to investigate the proportions of ancestry in French HGs derived from these two Late Pleistocene lineages, we used qpAdm to model it as a mixture of two sources: ∼15,000-y-old GoyetQ2 from Belgium, the least-admixed individual associated with the Magdalenian complex in Europe sequenced so far (11), and Villabruna. We found that French HGs from Les Perrats harbored relatively high proportions of GoyetQ2 ancestry ranging from 31.3 to 45.6% (Fig. 1C), comparable to proportions described in the La Braña or Canes1 Mesolithic individuals from Spain, suggesting a late survival of Magdalenian-associated ancestry in HGs outside the Iberian Peninsula.

Successive Waves of Migration and Admixture over the Course of the Neolithic in France.

The arrival of an Anatolian Neolithic-associated ancestry component at the onset of the Neolithic is clearly visible in our dataset from both uniparentally inherited markers and at the genome-wide level. Maternal lineages from Neolithic farmers in France (from ∼5300 BCE onward) are more diverse than those from HGs, and display variable frequencies across the time transect (SI Appendix, Fig. S3-2). Early and Middle Neolithic individuals share more affinity with present-day southern Europeans (SI Appendix, Fig. S4-2), and their genetic variation is encompassed within that of contemporaneous European populations (Fig. 1B). Three of our individuals associated with the LBK culture, which stems from the Danubian Neolithic current, cluster genetically with other central European Early Neolithic (EN) individuals and share drift with other early farmers associated with the LBK cultural complex (SI Appendix, Fig. S4-6). Values of D (Mbuti, France_EN; western hunter-gatherers [WHGs], Anatolia_Neolithic) are consistent with those observed in LBK-associated individuals, with the exception of Mor6 (∼7,100 calB.P.) in northeastern France (Fig. 2B). Despite a consistent dating and cultural assignment to the LBK, this individual falls within the genetic diversity of Iberian_EN individuals and harbors the highest proportion of shared alleles with WHGs. After testing the potential origin of this ancestry with qpAdm, we found the best fit to be with a two-source ancestry model between Anatolia_Neolithic and GoyetQ2 (P = 0.0911155), a situation so far only described in Early Neolithic individuals from the Iberian Peninsula. This is the northern-most description of such ancestry. To date, the distribution of the Magdalenian-associated ancestry in the HG population in western Europe that admixed with the Neolithic farmers in Europe is not clear because of the scarcity of genetic data from Mesolithic HGs. The observation of this ancestry in the Iberian Peninsula has been interpreted as a feature distinguishing the Cardial from the LBK Neolithic migrants (12).

The observation that the ancestry of GoyetQ2 was found alongside Villabruna ancestry in HGs from Les Perrats in western France (∼9,100 calB.P.), as well as in HGs from Rigney1 (central-eastern France; ∼15,500 calB.P.) and from two caves in southwestern Germany, Hohlefels (∼15,000 calB.P.) and Burkhardtshöhle (∼14,600 calB.P.), rather indicates that this mixed ancestry is characteristic of western European HGs (11). Thus, the admixture between HGs with GoyetQ2 ancestry and Neolithic migrants is not necessarily characteristic of the Cardial Neolithic migration wave, as it seems to have occurred as well in more northern locations in western Europe.

While the only male Mesolithic individual from our dataset was assigned to Y-haplogroup I, LBK individuals from France showed relative diversity, belonging to Y-haplogroups C1a2, G2a, and H2. G2a and H2 have been described in various European Early Neolithic contexts in both central Europe and Iberia (3, 6, 24). I2 haplogroups are likely to have been introduced into the Neolithic pool through admixture with hunter-gatherers.

In eastern France, the transition into the Middle Neolithic is represented by individuals from the Grossgartach culture, that derives from the Danubian sphere (∼4700 to 4500 BCE; SI Appendix, Fig. S1-1). Their proximity to central European contemporaneous individuals on PCA plots generated from both mitochondrial haplogroup frequencies and genome-wide markers suggests genetic homogeneity within the Danubian current (Fig. 1B and SI Appendix, Fig. S3-3). Interestingly, the BUCH2 individual harbors mixed GoyetQ2 and Villabruna HG ancestry (Fig. 1C), suggesting again that the GoyetQ2 ancestry is not a specific signature of southwestern Europe. However, we can observe a subtle shift away from LBK individuals and toward WHGs, reflecting an increase in shared alleles with the latter. This echoes observations from other ancient populations across Europe (Figs. 1C and 2A), yet with a higher proportion. This phenomenon seems to appear in France very early and becomes stronger as the Neolithic progresses, since individuals from French Middle Neolithic sites share a higher proportion of ancestry with HGs (Fig. 2A) and display Mesolithic signature haplogroups (e.g., U5b) at variable frequencies across the territory (SI Appendix, Fig. S3-2). Predominant in the Champagne region of northeastern France, these haplogroups suggest variable proportions of admixture between autochthonous hunter-gatherers and farmers. During the second half of the Middle Neolithic (MN2, ∼4500 BCE), individuals from northern and southern France are poorly distinguished on the PCA plot (Fig. 1B).

Unlike in Iberia, the migration wave from which the Neolithic ancestry in ancient populations from France originated could not clearly be established using D statistics comparing the ancestry proportion between Iberia_EN and LBK_EN, as most individuals display nonsignificant values (SI Appendix, Fig. S4-6) (12). Using qpAdm, we investigated the ancestral source of the HG ancestry in individuals from French MN2 sites. Only individuals from southern France harbor varying levels of GoyetQ2-like ancestry, suggesting admixture with different populations of hunter-gatherers (Fig. 1C).

The Bell Beaker Period in France.

At the end of the Neolithic, the Bell Beaker complex (named after the particular shape of its ceramics) appears and becomes widespread across Europe. The Beaker phenomenon is peculiar because of its large geographic distribution and because of its coexistence with local Late Neolithic and Copper Age cultures.

We report two Bell Beaker-associated individuals (CBV95 and PEI2), that we coanalyzed with previously reported contemporaneous individuals from Europe, including France (10). French Beaker-associated individuals display a wide range of steppe-ancestry proportions (Figs. 1C and 2D). CBV95 in northern France derives the highest proportion of alleles from the Yamnaya in our dataset, and belongs to Y-chromosome haplogroup R1b (Figs. 1C and 2D and SI Appendix, Fig. S4-5), providing the earliest clear evidence of the presence of this haplogroup in France around 2500 BCE (Dataset S10). This lineage was associated with the arrival of migrants from the steppe in central Europe during the Late Neolithic, and was described in other parts of Europe and in Bell Beaker-associated individuals from southern France, while being almost absent in Iberia prior to the Bronze Age (10, 13). PEI2, a male unearthed from a collective burial site near Carcassonne in southwestern France with artifacts of the Bell Beaker complex, falls within the genetic diversity of Neolithic individuals in the PCA. Modeling admixture proportions between three source populations, Anatolia_Neolithic, Villabruna, and Yamnaya_Samara, we could, however, detect 28.3% of steppe ancestry in PEI2 (Fig. 1C). These observations are consistent with previous findings and confirm that steppe ancestry appeared later and with a lower impact in southwestern Europe than in other parts of the continent (10, 25). Surprisingly, the admixture model did not include Villabruna as an ancestry source for either CBV95 or PEI2, which differs from previously known Late Neolithic individuals.

Relative Continuity between the Bronze Age and the Iron Age.

Bronze Age individuals in France display new mitochondrial haplogroups such as U2, U4, and I, although at a low frequency. First described in eastern European populations and in Late Neolithic central Europe (13), these haplogroups suggest gene flow with a population deriving maternal ancestry from the Pontic steppe herders at the onset of the Bronze Age. As foretold by its early occurrence in Beaker-associated CBV95, a drastic Y-chromosome turnover occurs during the Bronze Age, where R1b replaces the preexisting diversity of Neolithic lineages in our sampling (Bronze Age, 11/13 and Iron Age, 7/10 individuals having R1b; SI Appendix, Fig. S5-1 and Datasets S10 and S11).

Bronze and Iron Age France share a common space in the PCA plot, both shifted toward modern central Europe and falling within the genetic diversity of Bronze Age Britain and central Europe (Fig. 1B), with a homogenization of the steppe component (Figs. 1C and 2C). In contrast to what was described for central Europe (9), there is no further shift toward eastern Eurasian genotypes during the Iron Age (9). Instead, the steppe component, heterogeneously distributed between individuals during the Bronze Age (ranging from 30 to 70%), becomes homogeneous (Figs. 1C and 2C), and individuals from the Hallstatt and La Tène culture in the French territory display similar affinities toward both modern and ancient populations. This could indicate that the transition from the Bronze Age to the Iron Age in France was mostly driven by cultural diffusion, without major gene flow from an external population. This would be consistent with an archeological and linguistic hypothesis proposing that the Celts from the second Iron Age descended from populations already established in western Europe, within the boundaries of the Bell Beaker cultural complex (26–29). It is important to mention, however, that due to the relative genetic homogeneity among European populations by the Bronze Age, subsequent migrations between different parts of Europe could easily remain unnoticed at this level of coverage (10, 25, 30).

Analysis of Genetic Markers Associated with Phenotypes.

To follow the evolutionary timelines of genetic adaptation to changing lifestyles and environments, we analyzed 73 autosomal loci, associated with both physical and physiological traits, genotyped across 149 individuals (Fig. 3; see Dataset S3 for details on the targeted positions). Among the phenotypic traits of interest, we focused on genetic variants involved in eye and skin pigmentation (such as SLC24A5, SLC45A2, GRM5, HERC2, IRF4, and TYR). It has been previously shown that Mesolithic HGs and Anatolian Neolithic farmers differ in their pigmentation alleles, with the latter carrying derived alleles responsible for a light skin that are near fixation in Europe today (22), while HGs are described as having a rather dark skin and light eye color (6, 31). The Mesolithic individuals from whom we retrieved reliable genotype information after stringent filtering (one to three individuals) possessed the ancestral pigmentation variant for loci SLC45A2 and GRM5, which had been associated with darker skin. In contrast, the derived alleles of these two loci are found at significantly higher frequency in the French Neolithic population (SLC45A2: 24.9%; GRM5: 38.4%), although they remain far from the frequency reached by these alleles in present-day Europe, where they are near fixation in the population (93.8 and 68.9%, respectively). SLC24A5, described as the principal mutation associated with light skin pigmentation, is found in 96% of our Neolithic population, close to the 99% found in present-day Europe. These results suggest different evolutionary timelines for these mutations in France, where a similar phenotype (light skin pigmentation) possibly resulted from several selection events, although the genetic determinism of these traits is likely more complex and not encoded by these three loci alone (32). The Mesolithic individual who displayed coverage of variants involved in the pigmentation of eye color (HERC2 and IRF4) carried the alleles associated with blue eye color. At the time of the Neolithic, the frequency of these two variants reaches 29.6 and 24.4%, respectively. Interestingly, these two variants seemingly have different trajectories after the Neolithic. While HERC2 rises in frequency again from 37.5% during the Bronze Age to finally reach more than 63% in present-day Europe, IRF4 further decreases from 23.81% during the Bronze Age to reach a frequency of 11% at present. None of the genotyped Neolithic individuals carried the mutation responsible for the persistence of lactase in Europeans, consistent with a later origin for this mutation (13, 30). We also analyzed several innate immunity-related variants that show signs of recent positive selection in modern European populations (6, 33–35). For a number of genes involved in the immune response (NOS2A, TOLLIP, CCL18, STAT6, and IL3), the frequency of derived alleles was comparable to that found in present-day Europeans, indicating that selection on these loci predates the Neolithic, whereas frequencies of variants in the TLR1–6–10 locus, the Toll-like receptor family involved in the innate immune response, are much lower than those found in present-day Europeans (36). The TLR1–6–10 gene cluster is possibly associated with resistance to leprosy, tuberculosis, or other mycobacteria (36). Our data suggest that the onset of agriculture and the proximity to domesticated animals did not exert a selective pressure on this locus. Likewise, variants associated with celiac disease in SH2B3 or SLC22A4 that could have been the target for positive selection owing to the role they play in inflammation and the elimination of a wide array of environmental toxins (33) did not yet reach the frequency of modern populations.

Processing of Archaeological Samples.

A total of 243 individuals from 54 archeological sites across France passed screening for the preservation of ancient DNA and were processed for sequencing (SI Appendix, Fig. S1-1 and Dataset S1). These individuals covered a time period ranging from the Mesolithic to the Iron Age, with the majority of samples covering the Neolithic Period (∼5500 to 3500 BCE). As direct dates were not available for all samples, they were labeled and grouped for analyses based on geographical proximity and cultural assignment (SI Appendix, Fig. S1-1).

Targeted DNA Enrichment and Sequencing.

DNA was extracted from petrous bones and teeth and used to build double-stranded Illumina libraries with thorough decontamination procedures and optimized protocols in a dedicated facility at Institut Jacques Monod (SI Appendix, Text S2) (37). Sequence enrichment for both the whole mitochondrial genome and nuclear regions was performed through a capture approach using biotinylated RNA baits obtained through in vitro transcription of PCR-amplified target regions (37) (see SI Appendix, Text S2 for details on the procedure). Both enriched and whole-genome shotgun libraries were sequenced on an Illumina MiSeq or Illumina NextSeq at Institut Jacques Monod and Institut de Recherche Biomédicale des Armées, respectively.

Data Processing and Authentication.

Adapters were trimmed and overlapping paired-end reads were merged using leeHom (38). Merged reads were then quality-trimmed with cutadapt and only reads longer than 30 bp were retained. Mapping was performed with BWA version 1.2.3 (11) against the Cambridge reference sequence (39) for human mitochondrial DNA with a duplication of its first 100 bases at the end to ensure mapping of the reads overlapping the junction resulting from the virtual linearization of the circular mitogenome and to the human reference genome hs37d5. In both cases, PCR duplicates were removed using Picard’s MarkDuplicate (https://broadinstitute.github.io/picard/). Ancient DNA authentication was assessed using mapDamage (40) and contamination was estimated using both the mitochondrial genome and X chromosome (41).

Uniparental Haplogroups.

Consensus sequences of the mitochondrial genome were obtained using the software ANGSD version 0.910, relying only on reads with base quality above 20 and mapping quality above 30, and with more than 3× coverage. Mitochondrial haplotypes were determined based on PhyloTree phylogenetic tree (42) build 17 and by using the HAPLOFIND web application (43) and Phy-Mer (44). Unexpected or missing mutations were visually inspected in Geneious version R6 (45) to check if they could be the results of misincorporations in low-coverage regions. Sequences with more than 5% undetermined sites were removed. Y-chromosome haplogroups were retrieved using Yleaf software (46).

Nuclear Markers.

Target regions were summarized in a file, allowing their extraction using the mpileup tool from the SAMTools toolkit (47). Proper variant-calling files were obtained using bcftools call and filtered for a base quality of 20 and a minimum coverage of five reads using the SelectVariants tool from the Genome Analysis toolkit (48). Allelic count and frequency estimation was carried out with PLINK (49).

Reference Datasets.

We merged our dataset with published ancient individuals from western Eurasia as well as genotyping data from the modern populations of the Human Origins panel (50). We only used published data obtained from UDG-treated libraries. Genotype was called for the 1.2 million SNPs broadly used in ancient DNA in-solution capture procedures (22). A single allele was drawn at random for each position (minimum mapping and base quality of 30) using pileupCaller (https://github.com/stschiff/sequenceTools), rendering the individuals from the dataset homozygous for each locus. For genotyping data from modern individuals of the Human Origins panel, heterozygous sites in the eigenstrat file were randomly recoded as homozygous for the reference or the alternative allele using an awk script. The two files were then merged using PLINK version 1.9 (49).

Population Genetic Analyses.

We conducted PCA on ancient individuals, projecting them onto the genetic variation of present-day Eurasians using the lsqproject option from smartpca (51). We computed D and f statistics using Admixtools to estimate gene flow and shared drift between populations and to test tree topologies (52). We used qpAdm to estimate the ancestry proportions attributable to early Anatolian farmers, hunter-gatherers, and steppe herders in the individuals of this study (13).

Data Availability.

Sequencing data can be found in the European Nucleotide Archive (ENA) at EMBL-EBI under accession no. PRJEB38152.