Genome flux and stasis in a five millennium transect of European prehistory

The Great Hungarian Plain was a crossroads of cultural transformations that have shaped European prehistory. Here we analyse a 5,000-year transect of human genomes, sampled from petrous bones giving consistently excellent endogenous DNA yields, from 13 Hungarian Neolithic, Copper, Bronze and Iron Age burials including two to high (~22 × ) and seven to ~1 × coverage, to investigate the impact of these on Europe’s genetic landscape. These data suggest genomic shifts with the advent of the Neolithic, Bronze and Iron Ages, with interleaved periods of genome stability. The earliest Neolithic context genome shows a European hunter-gatherer genetic signature and a restricted ancestral population size, suggesting direct contact between cultures after the arrival of the first farmers into Europe. The latest, Iron Age, sample reveals an eastern genomic influence concordant with introduced Steppe burial rites. We observe transition towards lighter pigmentation and surprisingly, no Neolithic presence of lactase persistence.


Supplementary Table 2. Alignment statistics for all ancient samples sequenced on an
Ilumina Miseq platform. Fifty base pair (bp) single-end sequencing was used. Adapters were trimmed using cutadapt as in Supplementary Method 3.2, retaining sequences greater than 16bp in length. Sequences were aligned to the GRCh37 build of the human genome with the mitochondrial sequence replaced by the revised Cambridge reference sequence and duplicates were removed as outlined in Supplementary Method 3.2. N/A indicates that the library did not pass quality assessment (Supplementary Method 1.4) Depth of  coverage  KO1  8590752  C  T  T  IJ  2  14144593  A  G  G  IJ  1  14484379  A  C  C  I  1  14847792  A  C  C  I  1  15023364  T  C  C  I  1  16354708  G  A  A  I  1  19048602  G  A  A  I  1  16638804  A  G  G  I2  4  18700150  C  T  T  I2  2  7879415 A G  IJK  2  15023364  T  C  C  I  2  19038302  A  C  C  IJ  1  19166861  T  C  C  IJ  1  7879415  A  C  C  I2a  1  BR2  6753519  G  A  G  IJK  9  7173143  A  G  A  IJK  13  8590752  C  T  T  IJ  15  14144593  A  G  G  IJ  10  19038302  A  C  C  IJ  13  19166861  T  C  C  IJ  10  14237131  C  T  T  J  9  22749853  A  C  C  J  6  19179335  T  C  C  J  14  14969634  T  G  G  J2  7  22243566  C  T  T  J2a  6  22942897  T  C  C  J2a1  10  IR1  6753519  G  A  G  IJK  1  15472863  T  C  T  K  1  20837553  T  C  T  K  1  Reference allele is with reference to the GRCh37 build of the human genome. TT  TT  TT  TT  TT  TT  TT  TT  TT  TT  TT  TT  TT

SUPPLEMENTARY NOTES 1 -Archaeological context and Samples Information
The Great Hungarian Plain (Supplementary Fig. 1) is a major lowland in the Carpathian basin area which is connected with the Mediterranean, the Pontic steppe, and the loess-lands of central Europe 14,43 . This area represents a crossroads of human prehistoric cultures.
The earliest agricultural communities belong to the Körös culture (Criş in Romania), which is part of the Balkan Early Neolithic complex, the "First Temperate Neolithic" [10][11][12] . Early Neolithic burials across eastern Europe are characterized by a lack of cemeteries, in contrast to those of later Neolithic periods. Furthermore, the majority of burials are found inside what seem to be refuse pits. In most cases it is not possible to discern whether the artifacts recovered from these pits were simply deposited or discarded into the pit or were rather interred with the person as grave goods 44,45 . In the specific case of skeletons from Körös culture sites from Hungary only 13 of the 184 burials were accompanied by objects or fragments of objects that could be identified to belong to a grave assemblage 46 .
Around 5,600/5,500 cal. BC (calibrated years BC) major changes took place in the Carpathian Basin. In the Great Hungarian Plain the Körös culture was followed by the Linearbandkeramik (LBK) cultural complex, which consisted of two synchronous regional groups: the Alföld Linear Pottery (ALP) culture -also known in the literature as Alföld Linearbandkeramik Culture and Alföldi Vonaldíszes Kerámia (AVK) 14,15in the Great Hungarian Plain and the Eastern LBK variant in west Hungary 15 . The Western LBK group dispersed agriculture into the Central European loess belt and became the dominant LBK culture of Europe. During the period between 5,500-5,000 BC ALP presence in the Hungarian Plains expanded into the Tisza lake region (eastern Hungary) and was accompanied by a synchronous differentiation into local groups. Archaeological research indicates that in this region the Late Neolithic (ca. 5,000-4,500 BC) comprises the tell settlement culture of the Tisza-Herpály-Csőszhalom complex which is contemporaneous with the Lengyel Culture of Transdanubia. These two cultural complexes met in the foreland of the northern Hungarian Range. This cultural division had its roots in the Middle Neolithic development.
In the last phase of the Neolithic, around 4,500 BC, at the end of the Atlantic Period, the climate favorable for crop cultivation was replaced by a cooler sub-boreal climate. This change had a major impact on the economy of Late Neolithic societies, and resulted in a shift from crop agriculture to an increased focus on animal husbandry. In the course of complex processes new network-systems emerged with the formation of the Copper Age cultural entities on the Great Hungarian Plain 47 . While there is no apparent discontinuity in the cultural sequences of the Late Neolithic and Early Copper age, there is change in settlement structures and burial customs, which correspond to differences in social status and wealth 48 . The Copper Age period in the Great Hungarian Plain was divided traditionally into 3 phases: Tiszapolgár, Bodrogkeresztúr/Hunyadihalom and Baden. Based on radiocarbon dates the following chronological framework was proposed: Early Copper Age (4,500/4,400-4,000 BC), Middle Copper Age (4,000-3,600/3,500 BC), and Late Copper Age (3,600/3,500-2,600/2,500 BC) 49,50 .
During the period from the end of the fourth to the middle of third millennium BC North Pontic kurgans (burial mounds), characteristic for steppe, became distributed into the Tisza region. Supporting the archaeological evidence, isotopic analyses also indicate the appearance of migrant groups 51 . During the course of the Bronze Age between 2,700/2,500 and 900/800 tell settlements occurred repeatedly and with the enclosured sites, hillforts represented special cultural attitudes of the new sociopolitical entities. The deposition of metal "hoards" in special places -a general phenomenon of European Bronze and Iron Age societies -was a new custom of these times 52 .
Iron metallurgy first appeared in central Europe during the 1 st millennium BC. During the early phase of the Iron Age the regions east and west of the Danube were parts of two separate cultural provinces 53

Archaeological sites
Tiszaszőlős-Domaháza is a small site at the northernmost settlement of the Körös culture along the Middle Tisza River. The Körös cultural phase includes a pit, a house plan and some graves 21 . The house plan was very regular in shape, forming an oblong of 12×5.5m. It is directed NW-SE along its longer axis and consists of a thick layer of ceramic and mussel shell fragments and, in the centre of its southern half, the burnt remains of a hearth. Apart from the Körös culture pit and house other later occupational phases were excavated at the site including early LBK, Early Bronze Age and Medieval phases. A total of 7 individuals were exhumed from Neolithic contexts: two complete skeletons, two nearly complete skeletons identifiable by their dispersed body parts, and three single bones. One skeleton (grave 2-3), separated into two or three parts probably by later disturbance, was found in the lower 10-13 th layers of Pit 6. The two most complete skeletons (Grave 1 and 6) were found in the eastern and northern sections of Pit 6. These belonged to regular burials in flexed positions. The other two sets of human remains (graves 4 and 5) were exhumed from the same pit: a skull (Grave 4, specimen KO1, recovered at 140-155 cm depth, Supplementary Table 1) and a mandible (Grave 5). Other human remains were found in the area of the house plan and belonged to two separate individuals (Grave 7 and 8). The skull of Grave 4 was found in clear Körös context in the lower part of Pit 6. Radiocarbon dates provide a comprehensive chronological framework for the site. Polgár-Ferenci hát is an ALP site, which is located on the former bank of the Tisza floodplain. A contiguous surface measuring approximately 4 ha was investigated during two seasons 2001/2002. The site consists of a section of an elongated oval double ditch system [59][60][61][62] . Traces of intensive human activity, including the remains of several burnt houses and a thick layer of refuse, were found within this enclosure. A total of 113 burials were recovered. The graves were clustered in smaller and larger concentrations within the excavated settlement area. Their density is higher in the proximity of ditches and within the enclosure. The grave goods recovered included obsidian cores, beads, armrings made of spondylus shells and vessels containing red ochre. These point to long-range trade contacts. A special group of burials is represented by regular inhumation graves in which the skull was missing 62 . Recovered pottery found within the circular ditch system is of typical Middle Neolithic Tiszadob-Bükk-Esztár and Szakálhát-Vinča ceramic styles. Grave 325/500 (specimen NE1, Supplementary Table 1) is an oval pit, orientation SE-NW, of a child, which was laid on the left side in flexed position. Red ochre was found on the skull. The only grave good is a deep bowl with incised decoration, which was placed in front of the skull. Grave 839/1198 (specimen NE4, Supplementary  69,70 . This settlement had a monumental round structure (a ditch bordered by a palisade), which probably served as a ceremonial center on its eastern side and had an urn cemetery at its northern part. The inner part of the settlement was full of pits, but house structures were not preserved. According to the radiocarbon dates the Late Bronze Age settlement spans the period between 1,540-1,000 cal BC. More than twenty pits yielded human burials. Pit 2161 contained two skeletons: a 2-3 year old child (specimen BR2, Supplementary Table  1) and an adult female with a similar absolute date as the child (Deb11101, 3130±70 uncal. BP). The ceramics in the pit are typical for the Late Bronze Age Kyjatice Culture, which is also in good agreement with the absolute radiocarbon dates.
At the southern edge of the settlement a pre-Scythian cemetery was excavated that contained two groups of skeletons. The first, with 7 graves, was very poorly preserved, while the second (with 12 graves) was better preserved. Grave 2638 (specimen IR1, a 4-6 years old child, Supplementary Table 1) belonged to the second group 71 . The skeleton of the child was lying on its back in an extended position but the burial was disturbed by recent ploughing. The only excavated grave good was one broken bronze pin, but in other graves within the cemetery both ceramic and iron objects were excavated and these confirm the attribution of the cemetery to the Early Iron Age Mezőcsát Culture.

Radiocarbon Dating
Each of the 13 individuals reported in this work was directly radiocarbon dated at the Oxford Radiocarbon Accelerator Unit Research Laboratory for Archaeology, Oxford, UK. Details about radiocarbon dating and grave numbers of each sample are reported in Supplementary Table 1. Supplementary Fig. 1 shows the location of the archaeological sites and radiocarbon dates in a timeline.

Sample preparation
All bones/teeth were photographed (Supplementary Fig. 2) and UV-irradiated (10 minutes on each of two sides) prior to any manipulation. Their surface was then carefully removed, inside a UV cabinet, using a dremel engraving cutter attached to a dental drill. When necessary, a dremel diamond wheel was used for subsampling. Samples were subsequently exposed to UV light for 10 minutes on each of two sides in a cross-linker (Biometra Biolink, 5 lamps at 254 nm). Fine powder was obtained by grinding the samples with a Mixer Mill (MM 400, Retsch).
All tools used in this phase were consecutively cleaned with the following products: bleach, DNA-ExitusPlus™, Ethanol and then exposed to UV for 30 minutes per side.

Petrous bones versus non-petrous bones
In order to evaluate DNA preservation in petrous bone samples versus non-petrous bones, we investigated both types of bones for seven of the 13 individuals investigated (Fig. 1, Supplementary Table 2). We compared dental crowns (IDs 8.  (Fig. 1, Supplementary Table 2). For individual NE6 we investigated two different areas of the same petrous bone (IDs 14.2, 14.3). We found that the denser part gave higher endogenous DNA yields. Therefore, this part was sampled in all other specimens. Dental portions gave the highest percentages of human DNA within the non-petrous bone samples, but the petrous bone from the same individual gave from 4-to 16-fold higher percentages of human DNA. For other non-petrous bone elements this ratio could reach up to 183-fold (Fig. 1, Supplementary Table 2).

DNA extraction
DNA extraction was carried out with around 300 mg of bone powder (Supplementary Table 3) and following a silica-column based protocol based on 72 , as modified by 73 . Bone powder was incubated for 24 h at 55°C followed by 24 h at 37°C in 2 ml tubes with 1 ml of freshly prepared lysis buffer (final concentrations: Tris HCl pH 7.4 -20 mM; Proteinase K, recombinant, PCR Grade -0.65 U/ml; Sarkosyl ® NL 30 -0.7 %; EDTA pH 8 -47.5 mM) shaking at 300 rpm in a thermo mixer (Thermomixer comfort Eppendorf ® ). The samples were then centrifuged for 10 min at 13,000 rpm and the supernatant was removed. Fresh lysis buffer (1 ml) was added to the pellet and the incubation and centrifugation steps were repeated. Contaminant DNA is likely to be more easily accessible during the first lysis and the second lysis buffer is more likely to include a higher relative concentration of endogenous DNA 74 . Therefore supernatant from the second lysis step was used for subsequent analyses.
The supernatant was transferred to an Amicon ® Ultra-4 Centrifugal Filter Unit 30K, diluted with 3 ml of 10 mM Tris-EDTA Buffer and centrifuged for 20 min at 2,500 rpm or until a final volume of 250 μl was obtained. The flow through was discarded and 3 ml of 10 mM Tris-EDTA Buffer was added to the sample and centrifuged again for 20 -30 min at 2,500 rpm or until a final volume of 100 μl was reached. This volume was then transferred to a silica column (MinElute PCR Purification Kit, QIAGEN) and purification was carried out following manufacturers' instructions, except for the last step in which TWEEN ® 20 (at 0.05 % final concentration, Sigma-Aldrich) was added to the elution buffer (EB). Then 60 μl pre-heated EB (at 60ºC) was added to each column for the final elution. A 1 μl aliquot was used for DNA quantification with Qubit, Invitrogen, following the manufacturer's instructions (see Supplementary Table 3).

Library construction and amplification
Libraries were constructed from 50 μl of DNA extract, based on 75 , with the following modifications. Firstly, blunt end repair was performed using NEBNext ® End Repair Module (New England BioLabs Inc.), following manufacturer's instructions. Secondly, Bst activity was arrested through heat inactivation (20 min at 80°C).
Indexing PCRs were performed using Accuprime™ Pfx Supermix (Life Technology), with primer IS4 (0.2 μM) and an indexing primer (0.2 μM; 75 ); 3 μl of library was added to the PCR mix to a total volume of 25 μl. For amplification the following temperature cycling profile was used: 5 min at 95°C, 8-12 cycles of 15 sec at 95°C, 30 sec at 60°C and 30 sec at 68°C followed by a final extension of 5 min at 68°C. PCR reactions were purified using MinElute silica columns as described above.
Quality assessment of the amplified library was performed on an Agilent 2100 Bioanalyzer, using an Agilent DNA 1000 Kit and following manufacturer's instructions. The length of the amplified fragments, excluding the adapters, ranged from 20-270 bp, with a peak around 40 bp. Details about indices and PCR cycles used are shown in Supplementary Table 3.

Raw reads processing of Hiseq Data
Raw reads were filtered based on the indices used in library preparation, allowing a mismatch in the first base of the index sequence and a single mismatch at any other position 76,77 . The majority of template molecules were shorter than the length of the sequence reads (100 bp), so part of the P7 adapter was usually found at the 3′ ends of reads. These adapter sequences were trimmed using cutadapt v1.3 78 , allowing a minimum length of 34 bp after trimming (-m 34) and a minimum overlap of 1 bp between the adapter and the read (-O 1).
Considering the high frequency of misincorporation sites at 5' ends of the reads (C to T) in ancient DNA, as well as the complementary G to A at 3' ends 76,79,80 , two bases at the 5' and 3' ends were trimmed before mapping, following 77 , using the program seqtk (trimfq with options -b 2 -e 2; https://github.com/lh3/seqtk), for a final minimum length of 30bp.

Mapping
Sequence reads were aligned independently to the GRCh37 build of the human nuclear genome and the revised Cambridge reference sequence for the mitochondrial genome (NCBI accession number NC_012920.1) using Burrows-Wheeler Aligner (BWA) version 0.7.5a-r405 81 . Apart from disabling the seed option (-l 1000), the program was run with default parameters. Duplicate reads were removed using Samtools version 0.1.19-44428cd 82 and indels were realigned using GATK's RealignerTargetCreator and IndelRealigner 39 . Genomic depth of coverage was calculated using depth-cover (https://github.com/jalvz/depth-cover).
BAM files obtained from different sequencing lanes were merged using the MergeSamFiles Picard tool (http://picard.sourceforge.net). Duplicates were further removed and reads with a mapping quality greater than 30 were retained and used for all subsequent analyses. Alignment statistics and mean depth of coverage are presented in Supplementary Table 3 (details on mitochondrial DNA are reported in Supplementary Method 4).

Preparation of modern datasets
We created a single modern SNP reference panel, labeled HGDP+, incorporating selected individuals from the Human Genome Diversity Project (HGDP) 18 , a dataset used to elucidate minor the Jewish diaspora 19 and a dataset enriched with populations from the Caucasus 83 .
Autosomal SNPs with rs identifiers present in all three datasets were mapped to the GRCh37 build of the human genome and genotypes were assigned to the forward strand. Datasets were merged and filtered to exclude SNPs with a minor allele frequency (MAF) ≤ 0.5% and SNPs lacking a genotype call in ≥ 10% of individuals. Populations relevant to the analyses were selected from this dataset (Supplementary Table 4) and outliers were identified using EIGENSOFT's smartpca program (pairwise r 2 > 0.2 84 ) and removed. The final dataset contained 517,623 autosomal SNPs genotyped in 552 individuals.

SNP calling in ancient samples
To explore the ancient data in the context of modern variation we called genotypes at all positions that overlapped with the HGDP+ dataset. Different SNP calling procedures were followed for the high coverage and low coverage data.
For the high coverage genomes NE1 and BR2, diploid genotypes were called for all sites in the HGDP+ dataset using GATK's UnifiedGenotyper tool 39 . Genotypes were only called for alleles observed in the HGDP+ dataset and using only sequence data with a base quality ≥ 20 and depth ≥ 15. To exclude false genotype calls due to molecular damage, a minimum allele count of 4 was required for the T or A alleles of C/T or G/A SNPs, respectively. Filtered VCF files were transformed to PED format using VCFtools 85 and merged with the modern dataset using PLINK, retaining only SNPs present in the ancient sample (Supplementary Table 5). For principal component and admixture analyses one allele was randomly chosen at each called locus, creating homozygous genotype calls.
For the low coverage genomes, genotypes were obtained for principal components analysis (PCA) following 3,86 using the GATK 39 tool Pileup. The resulting pileup files were filtered for a minimum base quality of 15. Triallelic genotypes and genotypes that could be explained by molecular damage were removed (A where the alternative state is G and T where the alternative state is C). For sites with > 1× coverage, one base call at that site was randomly selected. Each hemizygous allele call was duplicated to form a homozygous diploid genotype and PED format files were created which were merged with the modern datasets as above. For each SNP genotyped in the ancient sample, one allele was randomly selected from each individual's diploid genotype in the modern dataset, simulating hemizygous data. This allele was then duplicated to generate a homozygous diploid genotype. Supplementary Table 5. These SNPs were used in downstream analysis such as PCA for all samples and runs of homozygosity for the high coverage genomes (Supplementary Method 6 and 5, respectively).

Independent replication
Results obtained in Dublin were partially replicated at the Department of Biology, University of York for two portions of the temporal bones from individuals NE1 and BR2, using around 250 mg each. DNA was extracted following 87 . Libraries were built following 75 , with the following modifications: (a) the P5 adapter included an index sequence of 8 nucleotides at the 3' end (NE1 -ACCAACCG, BR2 -AACGCATC), and (b) the filtration step between the blunt end repair and the adapter ligation was substituted by heat inactivation of the enzymes. Libraries were amplified in two rounds using AmpliTaq Gold DNA Polymerase. Eight parallel amplifications of 15 cycles each were carried out using the primers IS7-indPCR and IS8 75 . Subsequently, 2.5 µl of the library were used as template for a second amplification of 5 cycles with primers IS4 and P7-indexing (NE1 -ACCAACCG, BR2 -AACGCATC, the same indexes used for the P5 adapter).
These libraries were subsequently sequenced on a MiSeq platform and reads were processed following the procedures described in Supplementary Method 2. We obtained 679,400 (18.49%) and 1,027,012 (26.96%) reads aligned (duplicates removed) to the human autosomal genome for NE1 and BR2 respectively. GATK's UnifiedGenotyper 39 was used to call HGDP+ variants from quality trimmed bam files, restricting the panel of variants to those called in the high coverage genomes (Supplementary Method 3). As a measure of result reproducibility, SNPs called in the replication analysis were then compared to corresponding SNPs in the high coverage genomes. Due to the low coverage data obtained in the replication, only homozygous calls in the high coverage genomes were considered, resulting in a comparison of 2,249 genotypes for NE1 and 3,609 genotypes for BR2. Our results were highly congruent. A concordance rate of 99.69% was found for NE1 and a similarly high rate of 99.08% was recorded for BR2. Only one and four SNPs in the low coverage data could not be potentially explained by molecular damage for NE1 and BR2 respectively (0.04% and 0.11%) and might therefore be due to contamination or genotyping error. This test supports high reproducibility of our results (over 99%) and an extremely low laboratory contamination rate (below 0.11%).
Mitochondrial DNA reads were also separately aligned to the rCRS as described in Supplementary Method 2, giving mean coverages of 0.95x for NE1 and 1.42x for BR2. The same haplogroup determined using the high coverage data was also predicted for these low coverage replicated samples, albeit with a lower associated probability (NE1, U5b2c, probability 51.1%; BR2, K1a1a, probability 61.6%).

Negative controls analysis
During the screening MiSeq run (Supplementary Method 2.5), contamination controls were sequenced in parallel with the samples. One extraction control, one library control and two PCR negative controls were prepared for each seven samples. Moreover, air controls (1.5 ml / 2 ml tubes were opened in the drilling working area for 30 minutes, prior to initiation of sample preparation, to check the quality of the air) and water controls (1 ml of water was added to clean grinding jars and shaken, to check they were contamination free) were periodically prepared as well.
To estimate the laboratory contamination at each stage, we compared the number of reads alignable to the human genome obtained for the controls with the number of reads obtained for the samples analyzed here, normalized by the total number of reads obtained in each MiSeq run. The contamination fraction (C) was calculated as following: where R b denotes the number of control reads aligning to the human genome, R s the number of reads aligned to the human genome for the sample analyzed in parallel with the control, T the total number of reads obtained in the specific Miseq run and V the volume added to the pool loaded on the sequencer for controls (V b ) and samples (V s ). Results are reported in Supplementary Table 6. For each sample, detailed information about relevant controls is provided.

Estimation of molecular damage and sequence length
Sequence length and molecular damage were calculated by considering a subsample of 500,000 reads per sample. Adapter sequences were trimmed and reads with a minimum length of 30 bp were retained (2.4). These sequences were mapped to the human genome as described in Supplementary Method 2.2; however no mapping quality filter was applied. The resulting bam file was used as input for mapDamage 2.0 41,88 with a minimum base quality of 20 (option -Q20, Phred score) specified.
The mean and standard deviation of sequence lengths are reported in Supplementary  Table 7, while the C to T patterns at 5' end are reported in Supplementary Fig. 3.
Both results point to authenticity of ancient molecules, with an average length of 50-70 bp, typical for ancient DNA 89 , and an upsurge of C to T transitions at 5' ends. A reciprocal pattern of increased G to A transitions is observed at 3' ends. The values for 5' C to T transitions lie between 0.147 and 0.284, typical for samples > 500 years and up to 50,000 years old 90 .

Contamination estimates on the mtDNA chromosome
Mitochondrial DNA contamination was estimated using high quality reads that aligned to the rCRS sequence (Supplementary Method 2.2). Only bases with a quality greater or equal to 23, mapping quality above 30 and sequence length higher than 30 were analyzed. Bases that matched the consensus base call (Supplementary Method 4.2) at a position were considered primary bases. Bases that diverged from the consensus base were considered secondary bases. The contamination rate (as percentage of contamination, % C, Supplementary Table 8) was calculated as the average of the number of secondary bases at haplotype defining positions that could not be explained by molecular damage (C to T or G to A, >82.35% of the total number of mismatches), following 86 . Given the high transition / transversion ratio of mtDNA, also the estimates based on all sites (including C to T and G to A possible due to damage) are provided (as percentage of contamination including possible molecular damage, % C+MD, Supplementary Table 8).
The contamination fraction (C) was calculated as following: where D i is the depth of secondary base at a haplotype defining site where the secondary base cannot be attributed to damage; D t the total depth of coverage at that site and P t the total number of haplotype defining positions evaluated.
Our contamination estimate, which might include sequencing errors as well as possible heteroplasmy and NUMTs (mtDNA transpositions into the nuclear DNA), span from 0 to 0.81%, suggesting a very low mtDNA contamination rate. See Supplementary Table 8 for details.

Contamination estimates on the X chromosome
For the six male individuals studied here it was possible to estimate X chromosome contamination. The X chromosome is a haploid marker in males and thus any polymorphic sites identified may be due to contamination, DNA damage, sequencing errors or mis-mapping. Contamination is expected to be more transparent at polymorphic sites than at monomorphic sites as there is greater opportunity for variability in contaminating population(s). In contrast, sequencing error should affect all positions of the genome in the same way.
In order to assess the contamination rate we examined discordance in the rate of heterozygous calls between known polymorphic sites and their adjacent sites. We used a set of polymorphic markers, reported in the 1000 Genomes panel 28 , on the non-pseudoautosomal region of the X chromosome. The frequency of the minor allele in Europe was used to isolate sites, which exhibited polymorphism in European populations. As the samples originate from Hungary and have been mainly handled by European individuals, we investigated the probability of contamination from a European source.
The aforementioned 1000 Genomes dataset (957,967 SNPs) was filtered as follows: 1. Only bases with quality ≥ 20 were considered 2. Sites were required to have a minimum depth of coverage of 5 (for the ~20× sample BR2) or 3 (for the ~1× samples KO1, NE5, NE6, NE7, IR1) and a maximum depth of 40 3. Tri-allelic and tetra-allelic sites were removed 4. Sites in which the reference allele frequency was equal to the alternate allele frequency were discarded We performed Test 1 and Test 2 as described in 91 . For Test 1 we computed the number of primary (most common) and secondary (least common) alleles present at each SNP site as well as for the four sites upstream and downstream of this site (Supplementary Table 9) as a proxy for background error. A contingency table (Supplementary Table 10) was prepared for each sample.
Using Test 2 we remove the assumption of independent error rates at each site by randomly sampling a single read at each site. We computed again the number of primary (most common) and secondary (least common) alleles as well as the four sites upstream and downstream (Supplementary Table 9) and prepared a contingency table (Supplementary Table 10) for each sample.

Maximum Likelihood Contamination estimate on the X chromosome
We took the average observed probability of error (e) at adjacent sites as the background error rate. The probability of observing a secondary allele (P s ) at this SNP site is given by: where c is the contamination rate and f is the frequency of the secondary allele in Europe.
Assuming that the number of secondary alleles follows a binomial distribution, the likelihood of observing the distribution of our data is proportional to the probability of observing a secondary allele and the probability of observing a primary allele at each site: where s and p are the secondary and primary allele count at each site. The likelihood to observe the distribution of our data is given by the product of all these probabilities at each site: We transformed this product to the sum of its logarithm and by examining the maximum likelihood surface determined the contamination rate to be between 0.21 -1.20% (Supplementary Table 11).

Conclusions
Authenticity of the sequences obtained is of key importance in any ancient DNA study. However, the probability of contamination differs with the types of samples investigated, with ancient human samples being particularly prone to contamination 92 . Despite this, and the fact that early NGS studies on ancient human remains revealed problems of contamination 93 , many projects rely solely on a posteriori insilico assessment of contamination levels. We have applied both a priori (extraction and library construction controls) and a posteriori (bioinformatics assessment of contamination levels) controls. We also performed classical replication of results in a secondary laboratory, a widely accepted measure to control for contamination in pre-NGS studies, which has thus far not been applied to NGS studies. All our analyses reveal consistent and extremely low levels of contamination of ~ 1%, or often much less.

Molecular sex determination
Sex was determined by analyzing the ratio of X to Y chromosome reads following 42 . All samples could be clearly assigned to one of the two sexes (Supplementary Fig.  4).

Mitochondrial DNA
All samples were mapped to the mitochondrial DNA (mtDNA) genome and filtered as described in Supplementary Method 4.2. BAM files were analyzed with the online tool MitoBamAnnotator 94 , using default parameters (minimal base quality of 23, minimal base quality per read of 5, maximal mismatch fraction of 3/36, maximal indel fraction of 3/36, maximal number of N's in a read of 0) except for the minimal mapping quality (30) and minimal read length (30 bp), already filtered in the input BAM file (Supplementary Method 3.2).
The minimum depth of coverage for establishing the consensus haplotype was 3, after filtering for molecular damage (Supplementary Table 12). Corresponding haplogroups were obtained through the online tool Haplogrep 95 based on the build 15 phylogeny of PhyloTree 96 .
KO1 shows the rare maternal lineage R3, which has only been described in one Finn 97 , one Armenian (Family Tree database, http://www.familytreedna.com/) two Yakuts and one Bengali individual 98 . New data support the placement of this haplogroup within the R1 branch and it has been recently suggested to rename this haplogroup as R1b 98 . Little is known about this haplogroup from modern populations. The finding of an Early Neolithic specimen from this branch suggests a molecular dating closer to the upper boundary of the range proposed by 97 , 12,245.8+-6,768.4 BP.
KO2 belongs to haplogroup K, the most frequent haplogroup found in the time series, together with NE6, BR1 and BR2. This haplogroup is present in around 6% of Europeans and Near Eastern populations (88), reaching higher frequencies in the Druze (16%) 99,100 and the Ashkenazi Jewish (32%) populations 101 . It has been already described in ancient European Neolithic specimens 23,24,102-106 and the Late Neolithic/Copper Age specimen from the Alps, Ötzi, belongs to the same sub-clade K1 6 . The KO2 specific haplogroup is K1 and the sub-clades K1a and K1c were found also in later samples from the Hungarian timeseries: the Neolithic specimen NE6 (K1a3a3) and the Bronze Age BR1 (K1c1) and BR2 (K1a1a). This might suggest a certain maternal continuity through time on the Great Hungarian Plain.
Individual NE1, dated to the Middle Neolithic ALP culture, belongs to the U5b2c mtDNA haplogroup. The clade U5 has been associated with the Palaeolithic background of Europeans and is currently widely distributed in Eurasia, showing a frequency of around 9% in Europeans 99 . It shows frequency peaks in the Basque country 107 and in Scandinavia, with the sub-clade U5b being present in 41% of the Saami population 108 . It has been suggested that this haplogroup might have spread into Europe after the Last Glacial Maximum from Palaeolithic refugia 109 . The presence, at low frequencies, of this clade outside Europe, in northern Africa and the Near East is probably due to back-migration from the European continent 99,108 . These hypotheses emerging from modern population genetics seem to be supported by ancient DNA results: haplogroup U5 assignments have been described in all published European pools of Mesolithic specimens 5,8,86,110,111 . However, it has also been described in Neolithic populations at lower frequencies 102,103,105,106 . These data, together with our detection of this haplogroup in the Middle Neolithic Hungarian specimen NE1, might suggest population admixture at the beginning of the Neolithic.
NE2 together with CO1 show haplogroup H, the most common clade among current Europeans (40-60%), with decreasing frequencies towards the Caucasus and the Near East (88). However, it has not been observed in Mesolithic hunter-gatherers and shows a lower frequency among the first Neolithic farmers (~19%) 112 . It has been suggested that its frequency started to increase during the Middle and Late Neolithic 112 and this is consistent with our observations here.
The sample NE3 belongs to haplogroup X2, which is a widespread haplogroup around the globe. The majority of X-clade individuals in Europe, the Caucasus and the Near East belongs to the sub-lineage X2 113 . Modern population genetics suggest an association of this clade with a Palaeolithic background 113 , but it has been described only in post-Neolithic samples 106,114,115 .
NE4 and NE5 belong to haplogroup J (J1c and J1c1 respectively), a widespread haplogroup in Eurasia showing a frequency cline from the Near East (13%) towards Europe (9%) 99 118 . The authors of 24 suggested that this was the result of a Neolithic diffusion by pioneer farming groups. Considering that it has never been typed in pre-Neolithic specimens, it seems that this haplogroup had a strong association with the Linear pottery cultures, both LBK and AVK.
The youngest sample of our time series, IR1, belongs to haplogroup G, specifically the sub-clade G2a1. Haplogroup G is typically East Asian 119 . G2 is found in northern China and central Asia, reaching peak frequencies in South Siberian populations 119 . This points clearly to an Asian maternal ancestry of the specimen IR1. The archaeological context supports this hypothesis; IR1 has been associated with the pre-Scythians, a culture linked to immigration from the East 120 .  Supplementary Table 13 we report haplogroup assignments using between 6 and 12 informative positions per individual.

Y-chromosome
KO1 was assigned to haplogroup I2a. Y-chromosome clade I is widespread in Europe with frequency peaks of 40%-50% in Scandinavia and southern Europe 121 . We can see parallels with his mtDNA assignment, U5, as the Y-chromosome I haplogroup is also virtually absent outside Europe 122 . Modern population genetics suggests a European Upper Palaeolithic origin of haplogroup I 121 , which is currently supported by ancient DNA, as five Mesolithic specimens from central and northern Europe have been typed as belonging to the I lineage 5 . Our data are consistent with the above and suggest that KO1 had a Mesolithic paternal ancestry.
Haplogroup I was probably introduced there during the Neolithic by farmers with a Mesolithic paternal ancestry, and possibly reached high frequency due to genetic drift. NE7 also belongs to haplogroup I2a. However, it is interesting to note the contrast in maternal and paternal ancestry of this specimen, which might suggest patrilocality practices. NE7 showed a Y-chromosome lineage associated with the Palaeolithic background and a mtDNA clade (N1a1a1a) associated with first Neolithic groups in central Europe 24 .
Both NE5 and NE6 belong to the same haplogroup C6 (also called C1a2 or C7). Haplogroup C has been described mainly in southern and eastern Asia, but some branches have also been found in Oceania, Malaysia, the Americas and Europe. The distribution of this clade has been suggested to support an African origin followed by a single migration 'Out-of-Africa' firstly towards India and then to southeast Asia and the other continents 124 . Only recently a new European branch has been identified, C6 (initially called C7), in southern Europe 125 . Further investigation will be needed to determine if it represents an ancient European lineage or it constitutes a recent Asian signal 125 . The first hypothesis is supported by the identification of this subclade in an ancient Mesolithic specimen 7 from northern Spain. These data suggest: (a) this currently uncommon and not well known haplogroup was probably much more common in the past, (b) it was widespread in Europe, as it has been found in eastern and Western European ancient specimens and (c) its distribution and frequencies have been significantly changing through time.  132,133 . It has been suggested that this haplogroup originated 12,000-14,000 years ago in southeast Asia followed by a subsequent spread toward northeastern Europe and Siberia 133 . It is widespread in Siberia but absent in Native Americans, dating the spread of this haplogroup after the Bering strait crossing towards the Americas 133 . These data might suggest an Eastern paternal lineage of IR1, in agreement with the archaeological link with Steppe groups and also an Asian maternal ancestry based on its mitochondrial haplogroup G2a1.

-Imputation of ancient genomes
Imputation is a method commonly used in genomics to infer missing genotypes based on comparison of their surrounding haplotypes with those of a larger dataset (often a phased reference panel). Imputation software can typically take genotype probabilities or likelihoods as input data, and it has been shown recently by Pasaniuc et al. 27 that even low-density variant likelihood calls made with ultra low-throughput sequence data are sufficient for an imputation engine such as Beagle 40 to make accurate estimates of both missing and low coverage genotypes. Output data represent the probabilities of the possible genotypes; suitable thresholds can therefore be chosen to control the error rate in the resulting data.
Imputation may be an effective tool for the analysis of low-coverage ancient genomes, as a dense set of accurate diploid calls may be made using the available low-coverage sequence data.

Assessment of imputation using subsampled data
Imputation was evaluated using high-coverage sequence data from individuals NE1 and BR2, subsampled using SAMtools 82 Fig. 5 b and d). Consistent patterns were observed for both high-coverage genomes.
Taking NE1 1× data as an exemplar, imputed data were marginally more accurate if the SNP site had spanning sequence data (Supplementary Fig. 5 e). When percentage of genotypes meeting threshold was plotted against coverage for two representative imputed genotype probability thresholds (Supplementary Fig. 5 f), a plateau was observed in the increase of retained genotypes commencing around 1×. This corroborated our target depth of 1× for several of the samples sequenced; extra sequencing would yield a diminishingly small increase in data.
To assess whether bias is introduced during imputation, genome-wide imputed genotypes (genotype probability threshold 0.99) from subsampled NE1 and BR2 datasets were merged with the HGDP+ modern reference dataset (Supplementary Method 2) and assessed using principal components analysis (PCA, Supplementary Method 6) and runs of homozygosity (ROH, Supplementary Method 7). PCA plots created using the subsampled datasets and high coverage genotypes from the same sample were merged using Procrustes transformation (Supplementary Method 6). Imputed data fall within a range of genome-wide variation defined by principal components 1 and 2 that is virtually indistinguishable from that of their corresponding high-coverage genotypes (Supplementary Fig. 6). The total genome length under ROH using imputed and high coverage data were compared to the distribution of this statistic for the HGDP+ panel measured using the same loci. The percentile position of the ancient samples was found to be reasonably robust to the process of imputation from low coverage sequence. The range in each was within 2.7% (Supplementary Table 14, Supplementary Fig. 7)

Imputation of low-coverage genomes
For genomes sequenced at 1× and 0.1×, genotype likelihoods were called using GATK's UnifiedGenotyper 39 Method 9), a genotype probability threshold of 0.85 was used to determine the genotype.

Assessment of Ancient Samples in the Context of the 1000 Genomes Data
As imputation relies on the presence of representative haplotypes in the reference panel, a potential weakness of this approach could be that haplotypes found in ancient samples may represent extinct lineages that are absent in modern populations, resulting in decreased accuracy of imputation for divergent samples. As KO1 and IR1 appear as outliers when examined in the context of the HGDP+ dataset (Fig. 2) we analysed these samples, as well as NE1 and BR2, in the context of the 1000 Genomes data using PCA (Supplementary Fig. 8) We found that KO1, NE1 and IR1 all fall outside clusters of the 1000 Genomes individuals, suggesting that this dataset is as representative for KO1 and IR1 as it is for NE1. As ancient high coverage data becomes more accessible it will be possible to directly assess imputation accuracy for diverse samples.

-Principal Components Analysis
To examine the samples in the context of the variation found in the HGDP+ dataset we performed individual principal component analysis (PCA) on all 13 samples (Supplementary Method 2.3) using SMARTPCA software 3,84 . Individual PCAs were plotted with R 134 with PC1 and PC2 vectors transformed to match the configuration of the HGDP+ dataset (data not shown).

Procrustes transformation
In order to visualize non-overlapping data from all samples jointly without compromising the resolution of the picture, we used an approach based on Procrustes analysis 3 using the R package 'vegan' 135 . We plotted the average of the transformed coordinates for individuals in the HGDP+ dataset with each of the ancient samples transformed coordinates. The final plot, Fig. 2 Fig. 9) closely reflects the original PCA presented in Fig. 2, allowing cross validation of the alternative approaches.

-Runs of Homozygosity
Runs of homozygosity (ROH) are extended homozygous haplotypes that are identical by descent, arising due to consanguinity, restricted population size and selection events 137 . The portion of the genome under long ROH increases with inbreeding. A predominance of shorter ROH is indicative of long-term reduced population size 138 . Pemberton et al. 30 detected superimposed distributions of ROH length in human populations, with a cut-off of ~1.6 Mb distinguishing longer and shorter runs.

Length of ROH analysis
The ROH analysis was carried out with PLINK 136 with the following specifications:  A sliding window of 5000 kb (--homozyg-window-kb 5000)  A sliding window of 50 SNPs (--homozyg-window-snp 50)  One heterozygote allowed in a window (--homozyg-window-het 1)  A proportion of 0.05 of overlapping windows that must be called homozygous to define any given SNP as 'in a homozygous segment' (-homozyg-window-threshold 0.05)  A minimum number of 50 SNPs to be called as homozygous (-homozyg-snp 50)  A minimum number of 500 Kb to be called as homozygous (-homozyg-kb 500)  If two SNPs within a segment were > 100 kb apart, that segment was split in two (--homozyg-gap 100)  A required minimum density of 50 SNPs (--homozyg-density 50) Short (< 1.6 Mb) versus long (≥ 1.6 Mb) ROH were plotted in Fig. 5 b.
The extreme position of KO1 in Fig. 5 b suggests a small population size and higher degree of shared ancestry in its population of origin. The remaining samples are progressively closer to modern population variation, with the most recent samples dating to the Iron, Bronze and Copper Age (IR1, BR2, CO1) giving values typical of modern individuals.
Supplementary Table 14 shows the comparison between high coverage and imputed data for NE1 and BR2. When the total genome length under ROH in each is compared to the distribution of this statistic for the HGDP+ panel measured using the same loci, the percentile position of the ancient samples was found to be reasonably robust to the process of imputation from low coverage sequence. The range in each was within 2.7%.

Correlation between ROH and Age
In Fig. 5 b samples seem to be distributed according to their age. We investigated the correlation between sample age and the total length of ROH and we estimated the coefficient of correlation r 2 (0.4042) and the p-value (0.06573) fitting the points into a linear model (Fig. 5 a). Even though the p-value is borderline significant, we can see a trend for decreasing total Mb of ROH through time. These results probably reflect an increase in mobility, population size and outbreeding over time 139 .

ADMIXTURE software
The unsupervised clustering algorithm implemented in ADMIXTURE 29 was used to estimate the ancestral genetic components of the seven ~1× imputed samples plus the two 1× sub-sampled ~20× genomes along with 552 modern reference samples (Supplementary Method 7.1). Prior to the ADMIXTURE analysis, the SNP dataset was filtered in PLINK 136 to remove SNPs in pairwise LD using a r 2 value of 0.2. This analysis resulted in a final dataset of 60,824 SNPs. An exploratory ADMIXTURE analysis was first completed using 10 replicates for each value of K({1..10}), this analysis returned K=4 to be the number of clusters that provided the optimal minimisation of the CV error. The K=4 analysis was then repeated 100 times and the results analyzed in CLUMPP using the Greedy algorithm with 1,000 random permutations 3 . The 100 runs produced similar results with a maximum similarity coefficient H' of 0.9933177. The clustering results were then visualized using R 134 (Fig. 4).

NgsAdmix software
NgsAdmix allows for the production of admixture analyses from low coverage genomes 140 , instead genotype likelihoods are used. To assess the impact of imputation on the ADMIXTURE analysis, the NgsAdmix analysis was completed on the same samples.
Prior to the NgsAdmix analysis, the ~1× bam files used for the imputation were converted to genotypes likelihoods using ANGSD (http://www.popgen.dk/angsd/index.php/Main_Page), specifically to genotype likelihood beagle input file format. These files were then merged with the HGDP+ dataset in the same genotype likelihood format and NgsAdmix run with default settings 140 . The clustering results (Supplementary Fig. 10) reflect the ones obtained using ADMIXTURE (Fig. 4).

Pigmentation
Skin reflectance increases with distance from the equator 141 suggesting a selective response to restricted environmental ultraviolet radiation, required to synthesize vitamin D. Several genes involved in pigmentation have been investigated [142][143][144] and the available results suggest a complex genetic history of these genes, with different selective pressure in space and time 32,[145][146][147][148] .
Here we report the imputed genotype of three pigmentation genes (SLC24A5, SLC45A2, associated with skin color, and TYRP1, associated with iris and hair pigmentation; Fig. 3) in which three SNPs (rs1426654, rs16891982 and rs2733831 respectively) have been specifically selected in Europeans. The selected alleles of the former two are close to fixation in Europe. The time depths for selective sweeps of these three SNPs have been proposed as 11,000-19,000 years ago based on modern population data 33 , while recently it has been suggested that positive selection on these variants were still ongoing later, based on ancient samples from the Pontic-Caspian steppe dating 6,500-4,000 years ago 32 .
Within the Hungarian time series, we observe the appearance of the SLC45A2 selected allele in a homozygous state during the Middle Neolithic period, whereas the selected allele in the gene SLC45A2 appears even later, towards the end of the Neolithic (Fig. 3). This suggests in both cases that the selective sweep occurred later than previously hypothesized 33 , as suggested by 32 . The selective sweep at the gene TYRP1 is more difficult to identify, given that the current frequency of the selected allele is not close to fixation. However its frequency in Neolithic samples is 0.375 (6/16 alleles, Fig. 3) while in subsequent periods it doubles to 0.75 (6/8 alleles, Fig.  3), suggesting positive selection on this allele.

Lactase persistence
In mammals the ability to digest raw milk usually declines after weaning, but several human populations are lactase (LCT) persistent, with an ability to digest milk during adulthood 149 . Their current geographical distribution is not uniform and correlates with historical milk consumption 150 . In Europeans, a selective sweep is strongly associated with the T allele at SNP rs4988235 151,152 , which regulates LCT expression persistence. Its presence in either homozygous or heterozygous form allows milk digestion into adulthood. It has been suggested that the T allele was selected during the Neolithic, around 7,500 years ago in a region between the central Balkans and central Europe, possibly within the Neolithic LBK culture in central Europe 150 . However, this allele has so far not been detected in Neolithic LBK samples 5,153 , nor in Late Neolithic samples from northeastern Iberia 38 and only in 5% of the later Pitted-Ware Culture hunter-gatherers from Sweden 154 . Our samples confirm the