Sampling

Written informed consent was obtained from the Saami individual whose genome was analysed in this study, which was approved of by The Hospital District of Helsinki and Uusimaa Ethical Committee (Decision 329/13/03/00/2013) and the ethics committee of the University of Leipzig (Approval reference number 398-13-16122013).

Sampling and extracting ancient DNA requires a strict procedure in order to avoid contamination introduced by contemporary genetic material. For the 13 Iron-Age individuals from Finland available to us, the sampling took place in a clean-room facility dedicated to ancient DNA work, at the Institute for Archaeological Sciences in Tübingen. The preliminary workflow included documenting, photographing and storing the samples in individual, ID-coded plastic tubes and plastic bags. As a result of an early pilot study, the tooth samples we used were fragmented, and some of the dentine was removed. The remaining dentine was collected by carefully separating it from the enamel with a dentist drill and cooled-down diamond drill heads, rotated at a speed below 15 rpm, to avoid possible heat-caused damaging to the ancient DNA.

For samples from the sites of Bolshoy and Chalmny Varre, we used leftover tooth powder that was originally processed at the Institute of Anthropology at the University of Mainz for replication purposes as described in ref. 24. In brief, the sample preparation steps included UV-irradiation for 30–45 min, followed by gentle wiping of the surface with diluted commercial bleaches. The teeth were then sand-blasted using aluminium oxide abrasive (Harnisch & Rieth) and ground to fine powder using a mixer-mill (Retsch).

Radiocarbon date calibration

We calibrated the radiocarbon date of Bolshoy, reported in refs 24,55 as 3473 ± 42 years BP, using Intcal1356 as the calibration curve, using OxCal 4.357.

DNA extraction and library preparation

DNA from the six Bolshoy and the two Chalmny Varre samples was extracted in the ancient DNA facilities of the Max Planck Institute for the Science of Human History (MPI-SHH) in Jena, Germany. Extraction for the Levänluhta samples was similarly conducted in the clean-room facilities of the Institute for Archaeological Sciences in Tübingen. For each specimen, ~ 50 mg of dentine powder was used for an extraction procedure specifically designed for ancient DNA retrieval58. Extraction buffer containing 0.45 M EDTA, pH 8.0 (Life Technologies) and 0.25 mg/ml Proteinase K (Sigma-Aldrich) was added to bone powder and incubated at 37 °C with rotation overnight. The supernatant was separated from the pellet of bone powder by centrifugation (14,000 rpm). A binding buffer consisting of 5 M GuHCL (Sigma Aldrich) and 40% Isopropanol (Merck), together with 400 μl of 1 M sodium acetate (pH 5.5) was added to the supernatant, and the solution purified by spinning it through a purification column attached to a High Pure Extender Assembly funnel (8 min. in 1500 rpm, with slow acceleration). The column was then spun into a collection tube (1 min 14,000 rpm) 1–2 times to maximise the yield. This was followed by two subsequent washing steps of 450 μl of wash buffer (High Pure Viral Nucleic Acid Large Volume Kit) and two dry spin steps of 1 min centrifugation at 14,000 rpm. The final total volume of 100 μl eluate was reached by two separate elution rounds of 50 μl of TET (10 mM Tris-HCL, 1 mM EDTA pH 8.0, 0.1% Tween20), each spun for 1 min at 14,000 rpm into a fresh Eppendorf 1.5 ml tube. Negative controls (buffer instead of sample) were processed in parallel at a ratio of 1 control per 7 samples.

Of the 100 μl extract, 20 μl was used to immortalize the sample DNA as a double-stranded library. The procedure included a blunt-end repair, adapter ligation and adapter fill-in steps, as described by Meyer and Kircher59. During the blunt-end repair step, a mixture of 0.4 U/μl T4 PNK (polynucleotide kinase) and 0.024 U/μl T4 DNA polymerase, 1× NEB buffer 2 (NEB), 100 μM dNTP mix (Thermo Scientific), 1 mM ATP (NEB) and 0.8 mg/ml BSA (NEB) was added to the template DNA, followed by incubation in a thermocycler (15 min 15 °C, 15 min 25 °C) and purification with a MinElute kit (QIAGEN). The product was eluted in 18 μl TET buffer. The adapter ligation step included a mixture of 1× Quick Ligase Buffer (NEB), 250 nM Illumina Adapters (Sigma-Aldrich) and 0.125 U/μl Quick Ligase (NEB), added to the 18 μl eluate, followed by a 20 min incubation, and second purification step with MinElute columns, this time in 20 μl eluate. For the fill-in step, a mixture of 0.4 U/μl Bst-polymerase and 125 μM dNTP mix was added and the mixture then incubated in a thermocycler (30 min 37 °C, 10 min 80 °C). Libraries without Uracil-DNA-glycosylase (UDG) treatment were produced for all of the 13 extracts from Levänluhta. In addition, UDG-half treated libraries were produced for seven of the original 13 extracts from Levänluhta, and for all Bolshoy and Chalmny Varre extracts. To introduce the UDG-half treatment, an initial stage was included in the library preparation, in which 250 U USER enzyme (NEB) was added into the 20 μl of extract, followed by an incubation at 37 °C for 30 min, and then 12 °C for 1 min. This again was followed by the addition of 200 U UGI (Uracil Glycosylase inhibitor, by NEB) and another identical incubation to stop the enzymatic excision of deaminated sites, as described in60. For each library, a unique pair of eight-bp-long indexes was incorporated using a Pfu Turbo Cx Hotstart DNA Polymerase and a thermocycling program with the temperature profile as follows: initial denaturation (98 °C for 30 sec), cycle of denaturation/annealing/elongation (98 °C for 10 sec/ 60 °C for 20 sec/ 72 °C for 20 sec) and final extension at 72 °C for 10 min61. Bone powder from a cave bear was processed in parallel serving as a positive control. Negative controls for both extraction and library preparation stages were kept alongside the samples throughout the entire workflow.

Experiment efficiency was ensured by quantifying the concentration of the libraries on qPCR (Roche) using aliquots from libraries before and after indexing. The molecular copy number in pre-indexed libraries ranged from ~ 10E8 to ~ 10E9 copies/μl, indicating a successful library preparation, whereas the indexed libraries ranged from ~ 10E10 to ~ 10E12 copies/μl, stating an admissible indexing efficiency. The negative controls showed 4–5 orders of magnitude lower concentration than the samples, indicating low contamination levels from the laboratory processing stages.

The libraries were amplified with PCR, for the amount of cycles corresponding to the concentrations of the indexed libraries, using AccuPrime Pfx polymerase (5 μl of library template, 2 U AccuPrime Pfx DNA polymerase by Invitrogen, 1 U of readymade 10× PCR mastermix, and 0.3 μM of primers IS5 and IS6, for each 100 μl reaction) with thermal profile of 2 min denaturation at 95 °C, 3–9 cycles consisting of 15 sec denaturation at 95 °C, 30 sec annealing at 60 °C, 2 min elongation at 68 °C and 5 min elongation at 68 °C. The amplified libraries were purified using MinElute spin columns with the standard protocol provided by the manufacturer (Qiagen), and quantified for sequencing using an Agilent 2100 Bioanalyzer DNA 1000 chip.

For the modern Saami individual, total DNA was phenol-chloroform extracted and physically sheared using COVARIS fragmentation. A modified Illumina library preparation was performed using blunt-end repair followed by A-tailing of the 3’-end and ligation of forked adapters. Indexing PCR was followed by excision of fragments ranging from 500 to 600 bp from a 2% agarose gel.

Capture & sequencing

We used the in-solution capture procedure from ref. 62 to enrich our libraries for DNA fragments overlapping with 1,237,207 variable positions in the human genome4. The sequences used as bait, attached to magnetic beads, were mixed in with the DNA sample template in solution, and left to hybridize with the target DNA during a 24-hour incubation at 60 °C in a rotation oven. 4–6 samples were pooled in equal mass ratios at a total of 2000 ng of DNA. The captured libraries were sequenced (75 bp single-end, plus additional paired-end for the three non-UDG libraries of the Levänluhta individuals) on an Illumina HiSeq 4000 platform at the Max Planck Institute for the Science of Human History in Jena. Out of the 13 originally processed Iron-Age samples from Finland, seven proved to be of an adequate quality to be used in downstream analyses. The modern Saami genome was sequenced in on a Genome Analyser II (8 lanes, 125 bp paired-end) at the Max Planck Institute for Evolutionary Anthropology in Leipzig.

Processing of sequenced reads

We used EAGER63 (version 1.92.50) to process the sequenced reads, using default parameters (see below) for human-originated, UDG-half treated, single-end sequencing data, when processing the UDG-half libraries for all individuals. Specifically, AdapterRemoval was used to trim the sequencing adapters from our reads, with a minimum overlap of 1 bp, and using a minimum base quality of 20 and minimum sequence length of 30 bp. BWA aln (version 0.7.12-r1039, https://sourceforge.net/projects/bio-bwa/files)64 was used to map the reads to the hs37d5 human reference sequence, with a seed length (-l) of 32, max number of differences (-n) of 0.01 while doing no quality filtering. Duplicate removal was carried out using DeDup v0.12.1. Terminal base deamination damage calculation was done using mapDamage, specifying a length (-l) of 100 bp (Supplementary Table 1).

For downstream analyses, we used bamutils (version 1.0.13, https://github.com/statgen/bamUtil.git) TrimBam to trim two bases at the start and end of all reads. This procedure eliminates the positions that are affected by deamination, thus removing genotyping errors that could arise due to ancient DNA damage.

For three Levänluhta individuals exceeding the threshold coverage of 1% in the preliminary screening, we used the non-UDG treated libraries to confirm the authenticity of the ancient data. For these untreated libraries, two rounds of sequencing were carried out, which were processed using EAGER with the above parameters, but specifying a non-UDG treatment mode and setting the correct sequencing type between the libraries. The merged reads were extracted from the resulting bam files, and merged with the bam file containing reads from the single end sequence run using samtools merge (version 1.3)65.

The modern Saami genome was generated using Ibis for base calling and an in-house adapter trimming script. The resulting reads were then aligned to the hs37d5 human reference genome using bwa 0.5.9-r16 (parameters -e 20 -o 2 -n 0.01).

Genotyping

We used a custom program (pileupCaller) to genotype the 15 ancient individuals. A pileup file was generated using samtools mpileup with parameters -q 30 -Q 30 -B containing only sites overlapping with our capture panel. From this file, for each individual and each SNP on the 1240K panel, one read covering the SNP was drawn at random, and a pseudohaploid call was made, i.e., the ancient individual was assumed homozygous for the allele on the randomly drawn read for the SNP in question. PileupCaller is available at https://github.com/stschiff/sequenceTools.git.

For the three Levänluhta libraries that did not undergo UDG treatment, we only genotyped transversions, thus eliminating artefacts of post-mortem C- > T damage from further analyses.

The shotgun genome of the modern Saami individual was genotyped using GATK (version 1.3-25-g32cdef9) Unified Genotyper after indel realignment. The variant calls were filtered for variants with a quality score above 30, and a custom script was used to convert the variants into EigenStrat format.

The data were merged with a large dataset consisting of 3871 ancient and modern individuals genotyped on the Human Origins and/or 1240K SNP arrays, using mergeit.

Sex determination

To determine the genetic sex of each ancient individual we calculated the coverage on the autosomes as well as on each sex chromosome. We used a custom script (https://github.com/TCLamnidis/Sex.DetERRmine) for the calculation of each relative coverage as well as their associated error bars (Supplementary Figure 1, Supplementary Note 3 for more information on error calculation). Females are expected to have an x-rate of 1 and a y-rate of 0, while males are expected to have both x- and y-rate of 0.5 (ref. 49).

Authentication

We first confirmed that the deamination pattern at the terminal bases of DNA reads were characteristic of ancient DNA (1.04–4.5% for non-UDG libraries, and 4.7–9.5% for non-UDG libraries, see Supplementary Table 1), using mapDamage (version 2.0.6)66. We performed a number of different tests to ensure the authenticity of our ancient data. For male individuals, we investigated polymorphisms on the X chromosome27 using the ANGSD software package (version 0.910)67. This revealed robust contamination estimates for 2 male Bolshoy individuals, and 1 male Chalmny-Varre individual. All of these were below 1.6% contamination (Table 1). For the female individuals from these two sites, we note that they are projected close to the males in PCA space (Fig. 2a, Supplementary Figure 3), suggesting limited effects of potential contamination. In addition, we generated a PMD-filtered dataset for all individuals using pmdtools (version 0.60)30. PMD-filtering was done using a reference genome masked for all positions on the 1240K capture panel, to avoid systematic allelic biases on the analysed SNP positions. We set a pmd-threshold of 3, which, according to the original publication30, effectively eliminates potential modern contaminants based on the absence of base modifications consistent with deamination.

To provide a more quantitative estimate of possible contamination in females, we used the ContamMix program (version 1.0-10)29 for estimating mitochondrial contamination. We extracted reads mapping to the mitochondrial reference for each of the ancient individuals using samtools (version 1.3)65. We then generated a mitochondrial consensus sequence for each of the ancient individuals using Geneious (version 10.0.9, http://www.geneious.com,68), and calling N for all sites with a coverage lower than 5. Finally, all mitochondrial reads were aligned to their respective consensus sequence, using bwa aln (version 0.7.12-r1039)64 with a maximum number of differences in the seed (-k) set to 5 and the maximum number of differences (-n) to 10, and bwa samse (version 0.7.12-r1039)64. A multiple alignment of the consensus sequence and a reference set of 311 mitochondrial genomes69 was generated, using mafft (version v7.305b)70,71,72 with the --auto parameter. The read alignment, as well as the multiple alignment of the consensus and the 311 reference mitochondrial genomes were then provided to ContamMix. We report here the a posteriori mode of contamination, along with the upper and lower bounds of the 95% posterior interval (Table 1).

For additional authentication, we ran supervised ADMIXTURE28 (version 1.3.0) for all samples using the six present-day populations (Atayal, French, Kalash, Karitiana, Mbuti and Papuan) as defined genetic clusters, to locate any large differences in genetic clustering among individuals from the same site (Supplementary Figure 2). We tested the power of this method to detect contamination and find that it can detect contamination that is distantly related to the ancestries present within the test individuals already at rates of 5–8%, but lacks the power to identify contamination closely related to the test individuals (see Supplementary Note 2). We did not observe significant differences (within our resolution) in the ancestry patterns between the ancient individuals from the same site, with the exception of Levänluhta, where the individual sample JK2065 seems to derive from a different ancestry. We therefore assigned it a separate population label, Levänluhta_B in this study.

Finally, using smartpca, we projected PMD-filtered and non-filtered datasets on the same set of principal components constructed on modern European populations, to ensure that the ancient individuals remain projected in the roughly equivalent positions regardless of PMD-filtering. This was possible for all samples with UDG-half treatment, except for the individuals from Levänluhta, which represented too little damage for the PMD-filtering to be applied. Regarding this site, we therefore relied on the non-UDG libraries (using transversions only) that were generated for the three individuals used in the main analysis. We found that within expected noise due to a low number of SNPs, all samples show consistency between the filtered and non-filtered datasets, suggesting a low amount of contamination in all of the samples (Supplementary Figure 3a, b). Four additional individuals from Levänluhta were excluded from the main analysis and from this authentication test because of low coverage (< 15,000 covered SNPs) and lack of non-UDG libraries.

F statistics

All programmes used for calculating F statistics in this study can be found as part of the Admixtools package (https://github.com/DReichLab/AdmixTools)2,42.

We used qp3Pop (version 412) for all F 3 calculations.

F 4 statistics were calculated using qpDstat (version 711), and qpAdm (version 632)2 was used to estimate mixture proportions using the following: Sources (Left Populations): Nganasan; WHG; EHG; Yamnaya_Samara; LBK_EN. Outgroups (Right Populations): OG1: Mbuti; CHG; Israel_Natufian; Onge; Villabruna; Ami; Mixe. OG2: Mbuti; CHG; Onge; Villabruna; Ami; Mixe. OG3: Mixe; CHG; Israel_Natufian; Villabruna; Onge; Ami. OG4: Mbuti; Israel_Natufian; Onge; Villabruna; Ami; Mixe. OG5: Mbuti; Samara_HG; CHG; Israel_Natufian; Villabruna; Ami.

To ensure the outgroup sets had enough power to distinguish the ancestries present in the sources, we ran qpWave (version 410) using only the sources as left populations and each outgroup set as rights. All such qpWave runs were consistent only with maximum rank, meaning all outgroup sets had enough power to distinguish between the five different sources. All qpWave and qpAdm models were run using the option allsnps: YES. When the five-way admixture models provided by qpAdm had p-values above 0.05, but included infeasible mixture proportions and one of the sources was assigned a negative mixture proportion, we ran the model again with that source was excluded. For each Test population, if outgroup set OG1 did not produce a working full model (p < 0.05), we tried alternative outgroup sets with one right population removed. This resulted in outgroup sets OG2-4. In the case of Levänluhta, multiple outgroup sets produced working models, which are listed in Supplementary Data 4. The admixture model needing the minimum number of sources while still providing feasible admixture proportions is always shown. In the case of PWC from Sweden where none of the outgroup sets OG1-4 produced a working model, a revised set of right populations was used (OG5) which includes Samara_HG to provide more power to distinguish hunter-gatherer ancestries. We preferred models with OG1-4 over OG5 in general, because OG5 contains more ancient genomes with potential biases in the right populations, which more often causes failing models for modern Test populations. The excluded sources in the minimal models were specified as N/A (Supplementary Data 4). If either Yamnaya or EHG could be dropped (as is the case for Levänluhta), we show the model which is more consistent with previous publications3,7,8,45 in Fig. 4, but show both models in Supplementary Data 4.

Principal component analysis

We used smartpca (version #16000)73 (https://github.com/DReichLab/EIG) to carry out Principal Component Analysis (PCA), using the lsqproject: YES and shrinkmode: YES parameters.

For the Eurasian PCA (Fig. 2a), the following populations were used to construct principal components: Abkhasian, Adygei, Albanian, Altaian, Ami, Armenian, Atayal, Avar.SG, Azeri_WGA, Balkar, Balochi, Basque, BedouinA, BedouinB, Belarusian, Borneo, Brahui, Bulgarian, Buryat.SG, Cambodian, Canary_Islanders, Chechen, Chuvash, Croatian, Cypriot, Czech, Dai, Daur, Dolgan, Druze, English, Estonian, Even, Finnish, French, Georgian, Greek, Han, Hazara, Hezhen, Hungarian, Icelandic, Iranian, Italian_North, Italian_South, Japanese, Jew_Ashkenazi, Jew_Georgian, Jew_Iranian, Jew_Iraqi, Jew_Libyan, Jew_Moroccan, Jew_Tunisian, Jew_Turkish, Jew_Yemenite, Jordanian, Kalash, Kalmyk, Kinh, Korean, Kumyk, Kurd_WGA, Kyrgyz, Lahu, Lebanese, Lezgin, Lithuanian, Makrani, Mala, Maltese, Mansi, Miao, Mongola, Mordovian, Naxi, Nganasan, Nogai, North_Ossetian.DG, Norwegian, Orcadian, Oroqen, Palestinian, Pathan, Russian, Saami.DG, Saami_WGA, Sardinian, Saudi, Scottish, Selkup, Semende, She, Sherpa.DG, Sicilian, Spanish, Spanish_North, Syrian, Tajik, Thai, Tibetan.DG, Tu, Tubalar, Tujia, Turkish, Turkmen, Tuvinian, Ukrainian, Ulchi, Uygur, Uzbek, Xibo, Yakut, Yi, Yukagir.

For the West Eurasian PCA (Supplementary Figure 3a, b), the following populations were used to construct principal components: Abkhasian, Adygei, Albanian, Armenian, Balkar, Basque, BedouinA, BedouinB, Belarusian, Bulgarian, Canary_Islander, Chechen, Chuvash, Croatian, Cypriot, Czech, Druze, English, Estonian, Finnish, French, Georgian, Greek, Hungarian, Icelandic, Iranian, Italian_North, Italian_South, Jew_Ashkenazi, Jew_Georgian, Jew_Iranian, Jew_Iraqi, Jew_Libyan, Jew_Moroccan, Jew_Tunisian, Jew_Turkish, Jew_Yemenite, Jordanian, Kumyk, Lebanese, Lezgin, Lithuanian, Maltese, Mordovian, North_Ossetian, Norwegian, Orcadian, Palestinian, Polish, Russian, Sardinian, Saudi, Scottish, Sicilian, Spanish, Spanish_North, Syrian, Turkish, Ukrainian.

ADMIXTURE analysis

ADMIXTURE28 was run with version 1.3.0, following exclusion of variants with minor allele frequency of 0.01 and after LD pruning using plink (version 1.90b3.29)74 with a window size of 200, a step size of 5 and an R2 threshold of 0.5 (https://www.genetics.ucla.edu/software/admixture/download.html). Five replicates were run for each K value, with K values ranging between 2 and 15. The populations used were: Ami, Ami.DG, Armenian, Atayal, Atayal.DG, Balochi, Basque, BedouinB, Belarusian, Brahmin_Tiwari, Brahui, Chuvash, Croatian, Cypriot, Czech, English, Estonian, Even, Finnish, Finnish.DG, French, Greek, GujaratiB, Hadza, Han, Hungarian, Icelandic, Kalash, Karitiana, Lithuanian, Makrani, Mala, Mansi, Mansi.DG, Mari.SG, Mbuti, Mbuti.DG, Mixe, Mordovian, Nganasan, Norwegian, Onge, Orcadian, Papuan, Pima, Russian, Saami.DG, ModernSaami, Sardinian, Scottish, Selkup, Spanish, Ukrainian, Ulchi, Yoruba, ALPC_Hungary_MN, Baalberge_MN, Baltic_BA, Baltic_CCC, Baltic_CWC, Baltic_LN, BolshoyOleniOstrov, Bu_kk_Culture_Hungary_MN, ChalmnyVarre, CHG, EHG, Esperstedt_MN, Ganj_Dareh_Iran_Neolithic, Hungary_MN, Hungary_Neolithic, Iran_Chalcolithic, JK2065, Koros_Hungary_EN, Kunda, Latvia_HG3, Latvia_MN1, LBK_EN, LBK_Hungary_EN, Levanluhta, Narva, PWC_Sweden_NHG.SG, Scandinavia_LNBA, SHG, Sweden_HG.SG, TRB, Ukraine_HG1, Ukraine_N1, WHG, Yamnaya_Samara.

We find that K = 11 results in the lowest Cross-Validation error, as shown in Supplementary Figure 4b.

Y-chromosomal haplotyping

We assigned ancient males to Y haplogroups using the yHaplo program (https://github.com/23andMe/yhaplo)75. In short, this program provides an automated search through the Y haplogroup tree (as provided within yHaplo, as accessed from ISOGG on 04 Jan 2016) from the root to the downstream branch based on the presence of derived alleles and assigns the most downstream haplogroup with derived alleles. For about 15,000 Y-chromosomal SNPs present both in our capture panel and in two published datasets76,77, we randomly sampled a single base and used it as a haploid genotype. We used a custom script to convert EigenStrat genotypes to the yHaplo format. We report the haplogroup assigned furthest downstream provided by the program (Table 1). We also manually checked derived status and absence of mutations defining the designated haplogroup because missing information might lead to a premature stop in its automated search.

Mitochondrial haplotyping

We imported the trimmed mitochondrial reads for each individual with mapping quality >30 into Geneious (version 10.0.9, https://www.geneious.com)68 and reassembled these reads to the reference genome RSRS78, using the Geneious mapper, with medium sensitivity and 5 iterations. We used the in-built automated variant caller within Geneious to find mitochondrial polymorphisms with a minimum coverage of 3 and a minimum Variant Frequency of 0.67. The resulting variants were exported to Excel and manually compared to the SNPs reported in the online mtDNA phylogeny (mtDNA tree Build 17, 18 Feb 2016, http://www.phylotree.org/). Nucleotide positions 309.1 C(C), 315.1C, AC indels at 515-522, 16182C, 16183C, 16193.1C(C) and 16519 were masked and not included in our haplotype calls.

Phenotypic SNPs

We used samtools mpileup (version 1.3)65, filtering for map- (-Q) and base- (-q) quality of 30, deactivating per-Base Alignment Quality (-B), on the trimmed bam files, to generate a pileup of reads mapping to a set of 43 phenotypic SNPs4,40,41,79 that are part of our genome capture panel. A custom python script was used to parse the pileup into a table containing the number of reads supporting each allele (Supplementary Data 2).

Code availability

All software first described in this study is freely available from online repositories. Sex.DetERRmine: https://github.com/TCLamnidis/Sex.DetERRmine

ContaminateGenotypes: https://github.com/TCLamnidis/ContaminateGenotypes