High-resolution genomic history of early medieval Europe

Twigstats

Twigstats takes the Relate³² output format as input and allows the computation of f-statistics directly on genealogies, by using the inferred expected number of mutations on each branch as input, which is computed as the product of a prespecified average mutation rate per base per generation, the branch length and the number of bases each tree persists⁴³. Importantly, Twigstats computes f₂-statistics ascertained by an upper date threshold, such that only branches younger than this threshold are used. If a branch crosses the threshold, we use only the proportion of the branch underneath the threshold. Twigstats additionally enables us to specify a minimum derived allele frequency and lower date threshold. Twigstats can also compute f₂-statistics on age-ascertained mutations, which is particularly convenient for individuals not built into the genealogies.

The computed f₂-statistics are fed into ADMIXTOOLS2⁷⁰ to compute derived statistics. ADMIXTOOLS2 implements computation of genome-wide f₂-, f₃– and f₄-statistics, as well as qpgraph and qpAdm models. We implement the sample size correction as detailed in ref. ²¹. The f₂-statistics are computed in blocks, typically of prespecified centimorgan size or of prespecified physical distance. These blocks are used downstream in ADMIXTOOLS2 to compute standard errors using a block-jackknife approach. By default, we compute f-statistics only on internal branches and exclude singleton tip branches to increase robustness against sample age.

The optimal Twigstats time cut-off is a priori unknown; however, we develop a theory that predicts the optimal choice in a simple two-way admixture as a function of the admixture date, source split time and admixture proportion (Supplementary Note). In this case, the optimal cut-off equals approximately 1.4 times the split time between admixing source groups, depending on exact parameters in the model (Fig. 1b,c and Extended Data Fig. 2).

Non-negative least squares ancestry modelling

We implement an approach that uses genealogies to emulate the chromosome painting technique of identifying closest genetic relatives along the genome^1,2 to fit admixture weights. When applied to true genealogies in simulations, this approach represents an idealized version of this idea.

We implement this function in Twigstats, which, given known assignment of each sample to a population, identifies, at each position in the genome, the population with which a sample coalesces first. Our implementation takes a list of reference populations as input, such that any coalescences that do not involve these reference populations are ignored when traversing back in time through genealogical trees. If the first coalescence involves multiple different reference populations, this coalescence event will be assigned to each population with a weight proportional to the number of samples in each population involved in that event.

We then implement a second function in Twigstats to compute, for each target population and putative source populations, the proportion of the genome ‘painted’ by each of the reference populations. Given k reference populations, we denote by a_i the vector of length k storing these proportions for population i. We fitted our target population as a mixture of putative source populations using a non-negative least squares approach that finds a solution to the optimization problem \({\text{min}}_{0\le \varSigma {{\boldsymbol{\beta }}}_{{\mathcal{l}}}\le 1}{||{\bf{a}}}_{{\rm{t}}{\rm{a}}{\rm{r}}{\rm{g}}{\rm{e}}{\rm{t}}}-{\bf{A}}{\boldsymbol{\beta }}|{|}_{2}\), where A is a matrix storing \({{\bf{a}}}_{{\mathcal{l}}}\) for putative source populations as its column vectors with \(\ell \) indexing source populations and β are non-negative mixture weights.

Admixture simulations

We use msprime⁷¹ to simulate genetic variation data to test our approach. All simulation scripts are available at https://github.com/leospeidel/twigstats_paper.

f
₄-ratio admixture simulation

Our simulation in Fig. 1b and Extended Data Fig. 3b simulates five populations named PI, PO, P1, P2 and PX, where PO splits from all other populations 10,000 generations ago, P1 and P2 represent two proxy source groups that split from each other at 250 generations or 500 generations ago, PI splits from P1 100 generations ago and PX emerges from a pulse admixture between P1 and P2 50 generations ago. All populations have a constant diploid population size of 5,000, a variable human-like recombination map, in which our simulation only covers chromosome 1, and a human-like mutation rate of 1.25 × 10⁻⁸ mutations per base per generation. We additionally have a modified simulation with a lower mutation rate of 4 × 10⁻⁹ mutation per base per generation, emulating a transversions-only dataset, and a simulation in which P2 has a diploid population size of 1,000 in the last 50 generations, emulating a recent bottleneck in this population. We sample 20 haploid sequences from all populations. The ‘large sample size’ simulation samples 100 haploid sequences from all populations.

f
₄-ratio admixture simulation with genotype and phasing errors

We emulate the data quality we expect in imputed ancient genomes (Extended Data Fig. 3b). We implement a simple error model in which every haploid genotype at any segregating site can switch with a certain error probability. We can theoretically compute the predicted squared correlation coefficient (r²) between the true simulated genotypes and the genotypes that include error, stratified by minor allele frequency, to generate a plot similar to those used for evaluating imputation accuracy using downsampled high-coverage ancient genomes⁷² (Extended Data Fig. 3a). As imputation accuracy varies for each individual in real settings, we randomly sample the error probability for each individual uniformly between 1 × 10⁻⁴ and 1 × 10⁻³ (errors per SNP per haplotype). This yields r² curves that are comparable to those observed in real data. We additionally simulate a high error case, for which we sample error probabilities between 1 × 10⁻³ and 1 × 10⁻².

In real settings, we are additionally required to computationally phase genomes. We emulate this by combining two haploid sequences to construct a diploid individual. We then computationally rephase these diploid individuals without a reference panel. This approach is expected to result in suboptimal phasing and should therefore be well suited to test robustness to phase-switch errors.

qpAdm simulation

Our simulation in Extended Data Fig. 3c uses the simulation model and script provided with ref. ²³, although we changed this script to use the human hotspot recombination map. We simulate only chromosome 1. In the original simulation model, admixing sources split 1,200 generations ago, with admixture occurring 40 generations ago. We additionally simulate a version in which all population split times and admixture times are reduced by a factor of 5. We sample 20 haploid sequences per population.

Stepping-stone separation by distance simulation

We adapt the simulation model provided previously²³ to simulate a stepping-stone model of nine populations organized on a 1D grid, in which individuals are able to migrate between adjacent populations (Extended Data Fig. 3d). We changed this script to use the human hotspot recombination map and simulate only chromosome 1. We simulate under migration rates of 0.001 and 0.005, corresponding to average F_ST values of 0.01 and 0.002, respectively²³. We sample 20 haploid sequences per population. We then fitted population 4 using pairs of other populations as sources in a rotational qpAdm scheme such that unused populations are assigned to the reference set.

We expect that this simulation model violates qpAdm assumptions of no (or limited) gene flow after admixture between sources and reference groups. Consistent with this idea, qpAdm models are rejected (P = 4 × 10⁻³⁸ for migration rates of 0.001 and P = 5 × 10⁻⁸ for migration rates of 0.005) when using Twigstats with a cut-off of 1,000 generations. However, these are not rejected using regular qpAdm, including when migration rates are high (and, therefore, F_ST is low), indicating that Twigstats is better powered to detect such scenarios of continued migration. Encouragingly, a model that involves the two immediately adjacent populations is selected in all replicates as the ‘best’ model (highest qpAdm P value) using Twigstats, whereas this is the case in only 80% (migration rate of 0.001) and 30% (migration rate of 0.005) of all replicates using regular qpAdm.

Neanderthal admixture and deep structure simulation

Our simulation in Extended Data Fig. 5d emulates Neanderthal admixture, in which Neanderthals and ancestors of modern humans split 25,000 generations ago and admixture occurs 2,000 generations ago. The resulting admixed non-African-like population coexists with the non-admixed African-like population until the present day. Furthermore, two Neanderthal populations split from each other 7,000 generations ago, which can be interpreted as emulating the Altai and Vindija Neanderthal populations, with Vindija being closer to the admixing source.

We simulate an alternative model with two subgroups emulating ancestral modern humans in Africa that have a non-zero symmetric migration rate, ranging from 4 × 10⁻⁵ to 2 × 10⁻⁴ per generation, up until 3,000 generations before present. One of these subgroups gives rise to a present-day African-like population, while the other gives rise to a present-day non-African-like population. We further sample two Neanderthal populations that split 7,000 generations ago and merge 25,000 generations ago with the same ancestral modern human subgroup that will eventually give rise to a non-African-like population.

We simulate whole genomes with human-like recombination rates and a mutation rate of 1.25 × 10⁻⁸ mutations per base per generation. Diploid effective population sizes are set to 10,000 except on the Neanderthal lineage, in which it is set to 3,000. We sample 2 haploid sequences for each Neanderthal population and 20 haploid sequences for the target admixed population and African non-admixed population.

Fine-scale structure simulation

Our simulation in Extended Data Fig. 5a emulates the emergence of a fine-scale population structure and is adapted from ref. ³⁹. In this simulation, populations split 100 generations ago into 25 subpopulations followed by a period in which individuals are allowed to migrate at a rate of 0.01 between adjacent populations in a 5 × 5 grid. The diploid effective population size is 500 in each of the 25 populations, and 10,000 in the ancestral population. We simulate ten replicates of chromosome 10, with a human-like mutation rate of 1.25 × 10⁻⁸ and hotspot recombination map. We sample two diploid individuals from each population. Furthermore, we sample 100 individuals from an ancestral population that splits from the 25 target populations 100 generations ago, before the emergence of structure in these 25 populations. Relate trees are inferred assuming true mutation rates, recombination rates and average coalescence rates across all samples.

Ancient sample selection

A full list of ancient genomes can be found in Supplementary Table 1. Published ancient shotgun genomes provided by refs. ^7,8 were only available aligned against the GRCh38 reference sequence. These data were realigned to the GRCh37d5 reference sequence using bwa aln (v. 0.7.17-r1188).

We select genomes with average autosomal coverage above 0.5×, except for VK518, which has previously been suggested to be of Saami ancestry⁶ and which had a coverage of 0.438. We included VK518 in our panel to capture this ancestry. Genomes above a coverage cut-off of 0.5× have previously been shown to result in reliable imputation results⁷². We exclude samples with evidence of contamination. We remove any duplicate individuals, such as individuals who were resequenced, choosing the file with the highest coverage. We then filter out any relatives annotated in the Allen Ancient DNA Resource v. 54.1²⁷, retaining the individual with the highest coverage in each family clade.

Our final dataset includes 1,556 ancient genomes.

Imputation of ancient genomes

We follow the recommended pipeline of GLIMPSE⁷³ and first call genotype likelihoods for each genome in the 1000GP, segregating sites using bcftools mpileup with filter -q 20, -Q 20 and -C 50. We subsequently impute each genome separately using GLIMPSE v. 1.1.1 using the 1000GP phase 3 reference panel⁷⁴ downloaded from https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/. These imputed genomes are merged into a single VCF (variant call format) for further downstream processing.

We filter any site for which more than 2% of sites have an imputation posterior of less than 0.8 and retain all remaining sites so as not to have any missing genotypes at individual SNPs.

Relate-inferred genealogies

We merge imputed ancient genomes with a subset of the 1000GP dataset, including all European populations (CEU, Utah residents with northern and western European ancestry; CHB, Han Chinese in Bejing, China; FIN, Finnish in Finland; GBR, British in England and Scotland; BS, Iberian populations in Spain; TSI, Toscani in Italy, YRI, Yoruba in Ibadan, Nigeria). We create a second dataset in which we merge imputed genomes with the Simons Genome Diversity Project⁷⁵ (SGDP) downloaded from https://sharehost.hms.harvard.edu/genetics/reich_lab/sgdp/phased_data2021/. These two datasets contain, respectively, a total of 2,270 and 1,834 modern and ancient individuals.

We then infer genealogies for the joint dataset of ancient and modern genomes using Relate v. 1.2.1. We restrict our analysis to transversions only and assume a mutation rate of 4 × 10⁻⁹ mutations per base per generation and input sample dates as shown in Supplementary Table 1. We use coalescences rates pre-inferred for the 1000GP and SGDP datasets.

MDS analysis

We compute f₂-statistics using the Twigstats function f2_blocks_from_Relate between all pairs of individuals and between all individuals and an outgroup (Han Chinese people in SGDP) using the Relate genealogies of SGDP modern and imputed ancient genomes. We set the argument t to specify a time cut-off and set the argument use_muts to FALSE to compute these f-statistics on branches of the genealogy and to TRUE to compute these only on the mutations. We use these to compute f₃(outgroup, indiv1, indiv2) = 0.5 × (f₂(outgroup, indiv1) + f₂(outgroup, indiv2) − f₂(indiv1, indiv2)) for every pair of individuals, and store 1 − f₃(outgroup, indiv1, indiv2) in a symmetric N × N matrix (where N is the number of individuals) for which we then compute an MDS using the R function cmdscale.

qpAdm modelling

In brief, qpAdm models are a generalization of f₄-ratios, for which one-, two- and three-source models can be tested as hypotheses and admixture components and their s.e. obtained with a block jackknife¹³. A qpAadm model is fully specified by a set of putative source groups and additional ‘outgroups’ that are used to distinguish source ancestries. We used a rotating approach in which we iteratively selected a subset of source groups and used all remaining putative sources as outgroups. This approach penalizes models where true contributing sources are used as outgroups. With sufficient statistical power, qpAdm models will be statistically rejected if true contributing sources are used as outgroups. If statistical power is more limited, several models will fit the data, but the correct model is expected to be preferred over wrong models. Throughout, we use the Relate genealogies of SGDP modern and imputed ancient genomes in our qpAdm modelling and first compute f₂-statistics using the Twigstats function f2_blocks_from_Relate between all populations involved, which we then feed to the ADMIXTOOLS2 package⁷⁰.

Clustering using qpwave

To overcome challenges with hand-curating source groups used in qpAdm modelling, we follow ref. ⁵ and run qpwave using Twigstats between pairs of ancient individuals. We use Han Chinese individuals from Beijing and five European populations from the 1000GP as reference groups. This approach tests whether two individuals form a clade with respect to reference groups. The reason why this is a principled approach despite the 1000GP groups post-dating the ancient individuals is that if a group of ancient individuals are truly homogeneous, they will be so also with respect to later individuals.

We then define clusters by running UPGMA (unweighted pair group method with arithmetic mean) on −log₁₀[P values] obtained from qpwave between all pairs of individuals and cut the resulting dendrogram at a height corresponding to a P value of 0.01. We then further subdivide clusters by requiring all samples to be within 500 years of the mean cluster age.

To choose the source groups shown in Fig. 2a and Extended Data Fig. 1d, we run this algorithm on samples from Iron and Roman Age Europe (Supplementary Table 1). We retain groups that have at least three individuals and, therefore, exclude clusters of size one or two.

This approach results in two clusters in the Scandinavian Peninsula, approximately separating northern from southern Scandinavia, three clusters in Poland and Ukraine that separate samples temporally between the early and later Bronze Age, a cluster combining the Hungarian Scythian and Slovakian La Tène-associated individuals, and a cluster each for Iron and Roman Age Portugal, Italy and Lithuania. In present-day Austria, Germany and France, this approach identifies three clusters, with each cluster spanning multiple archaeological sites in different countries, indicating genetic diversity in this region in the first millennium ce. Encouragingly, these clusters separate in our non-parametric MDS analysis (Fig. 2a), indicating that we are capturing real genetic differences between groups using this approach.

Fine-scale structure in Neolithic Europe

To quantify fine-scale structure in Neolithic Europe (Extended Data Fig. 5b), we aimed to select individuals in Neolithic Europe who have not yet been affected by the arrival of Steppe ancestry and do not show excess hunter-gatherer ancestry. We infer distal ancestry sources using Balkan_N, Yamnaya and Western Hunter-gatherers as source groups and reference groups according to a previously proposed qpAdm setup⁴⁶ (Supplementary Table 1). For this analysis, we infer ancestry using qpAdm applied to 1.2 million SNP sites of imputed genomes. We retain only Neolithic individuals with P > 0.01, z < 2 for Yamnaya ancestry, and z < 2 or proportion <0.25 for Western Hunter-gatherer ancestry.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.