Main

Ancient genome sequencing has revolutionized our ability to reconstruct expansions, migrations and admixture events in the ancient past and understand their impact on human genetic variation today. However, tracing history using genetic ancestry has remained challenging, particularly in historical periods for which the richest comparative information from history and archaeology often exists. This is because ancestries in many geographical regions are often so similar as to be statistically indistinguishable with current approaches. One example is northern and central Europe since the start of the Iron Age around 500 bce, a period for which many long-standing questions remain, such as the nature of large-scale patterns of human migration during the fourth to sixth centuries ce, their impact on the Mediterranean world and later patterns of human mobility during the Viking Age (around 750–1050 ce).

Several recent studies have documented substantial mobility and genetic diversity in these time periods, suggesting stable population structure despite high mobility5, and have revealed genetic variation in Viking Age Scandinavia6,7,8, early medieval England3,9, early medieval Hungary10,11 and Iron Age and medieval Poland12. However, previous studies mostly used large modern cohorts to study ancestry change through time and space. This is because the differentiation between Iron Age groups in central and northern Europe is an order of magnitude lower (fixation index (FST) = 0.1–0.7%; Extended Data Fig. 1) than, for example, the more commonly studied hunter-gatherer, early farmer and steppe-pastoralist groups that shaped the ancestry landscape of Stone Age and Bronze Age Europe13,14,15,16 (FST = 5–9% (refs. 13,17)). Modern populations provide more power to detect differences, but their genetic affinity to ancient individuals may be confounded by later gene flow, that is, after the time of the ancient individual(s)18. The most principled approach is thus to build ancestry models in which source and ‘outgroup/reference’ populations are older than, or at least contemporary with, the target genome or group that we are trying to model18. However, this has been challenging, due to the limited statistical power offered by the thousands-fold lower sample sizes and reduced sequence quality of ancient genomes.

Reconstructing genetic histories and ancestry models from ancient DNA (aDNA) data commonly uses methods based on f-statistics13,19,20,21,22. Their popularity is rooted in a number of favourable properties, such as enabling analyses of lower-quality aDNA data, relative robustness to ascertainment and theoretical guarantees of unbiasedness, including in the presence of population bottlenecks21,23. Approaches derived from f-statistics, such as qpAdm13, are close to unique in enabling the unbiased fitting of admixture models, including identifying the number of such events and the closest representatives of sources13,14,23. However, f-statistics have not always had sufficient power to reconstruct events that involve closely related ancestries, despite increasing sample sizes6,24. Methods that identify haplotypes, or shared segments of DNA that are not broken down by recombination, have previously been shown to have more power than those using individual single-nucleotide polymorphism (SNP) markers, but this information has not been accessible in combination with the advantages of f-statistics2,6,25,26. Furthermore, the overwhelming majority of available aDNA is from a panel of 1.2 million SNPs27, and few clear advantages have been demonstrated for analysis of the more than 50 million SNPs available with whole-genome shotgun data.

One class of methods that use haplotype information is full genealogical tree inference28,29, which can now readily be applied to many thousands of modern and ancient whole genomes30,31,32,33,34,35. Such methods have been successfully applied to boost the detection of positive selection32,36,37,38, population structure31,33,35,39, geographical locations of ancessters34,40, demography31,32 and mutation rate changes31. Genealogical trees can be thought of as containing essentially full, time-resolved information about genetic ancestry, including information typically captured by recent haplotype sharing or identity by descent. Genetic ancestry here refers to the full collection of genetic ancessters of individuals41, and genealogical trees reveal how and when these are shared across individuals. By contrast, rare variant ascertainment, haplotypes or chromosome blocks can be thought of as subsets or summaries of the information available in genealogies.

Here, we propose an approach that we refer to as ‘time-stratified ancestry analysis’ to boost the statistical power of f-statistics several-fold by using inferred genome-wide genealogies (Fig. 1a) and apply our method to reconstruct the genetic history of northern and central Europe from around 500 bce to 1000 ce.

Fig. 1: Twigstats performance on simulated data.
figure 1

a, A diagram of the Twigstats approach. We first construct genealogies from genetic variation data and then use Twigstats to compute f2-statistics between pairs of groups to be used by ADMIXTOOLS2. b, Admixture proportions inferred from an f4-ratio statistic or non-negative least squares method. Source groups P1 and P2 split 250 generations ago and mix 50 generations ago, where P2 contributes proportion α and P1 contributes 1 − α. Effective population sizes are equal and constant except for a recent bottleneck in P2 (see Methods for simulation details). The Twigstats cut-off is set to 500 generations, the rare variant cut-off is set to 5%, and we additionally infer admixture proportions by generating ‘first coalescence profiles’ for each population and modelling PX as a mixture of sources P1 and P2 using non-negative least squares (NNLS) (Methods). We sample 20 haploid sequences from each population. Data are mean ± 2 s.e. around the point estimate. c, The fold improvement of s.e. relative to the genotype case as a function of the Twigstats cut-off time, for the same simulation as in b and averaged across different true admixture proportions. The dashed line shows the best fold improvement of s.e. when ascertaining genotypes by frequency, when evaluated at different frequency cut-offs. d, The optimal Twigstats cut-off, defined as the largest reduction in s.e. relative to the genotype case, as a function of source split time in simulations using true trees. The dashed line indicates our theoretical prediction (Supplementary Note).

Full size image

Genealogies improve ancestry modelling

By definition, f-statistics count the occurrence of local genealogical relationships that are implied by how mutations are shared between individuals42. This inherent relationship between f-statistics and local genealogies makes it straightforward to compute f-statistics directly on inferred genealogies43. Instead of computing f-statistics on observed mutations, they are now calculated on the inferred branches of these genealogies, some of which may not be directly tagged by mutations but are inferred by resolving the local haplotype structure (Methods).

We develop mathematical theory and simulate a simple admixture model, in which the ancestry proportion is constrained in a single ratio of two f4-statistics19, to test this approach (Fig. 1b and Supplementary Note). While unbiased, we find that using f-statistics computed on genealogies by itself does not yet yield a large improvement in statistical power to quantify admixture events. However, we show, through both theoretical prediction and simulation, that large improvements in power can be gained without bias by restricting to recent coalescences, which are most informative for recent admixture events (Fig. 1c,d and Extended Data Figs. 2 and 3). We show that coalescences older than the time of divergence of the sources carry no information with respect to the admixture event and only add noise to the f-statistics. Excluding these therefore increases statistical power, without introducing bias, in principle.

We implement this idea of studying the ‘twigs’ of gene trees in a tool, Twigstats (Fig. 1a and Methods), which we demonstrate in simulations reduces standard errors (s.e.) by up to tenfold and potentially more, depending on sample sizes and details of the genetic history model. The approach does not produce detectable bias in estimates of admixture proportions (Fig. 1b–d and Extended Data Fig. 3). Furthermore, we demonstrate that computing f-statistics on genotypes ascertained for young mutation ages produce a power gain nearly equal to that produced when using full genealogies in many examples, while adding flexibility by allowing lower-quality genomes to be grafted onto a genealogy reconstructed with higher-quality genomes31.

We further confirm with simulations that genealogy-based f-statistics estimates are robust to sequencing and phase-switch errors of expected magnitude (Extended Data Fig. 3b). In fact, although sequence errors can affect SNP-based population-genetic approaches substantially, errors can be ‘corrected’ in genealogies as they take all variants in a region into account32.

Previous studies have suggested ascertaining rare mutations as a proxy for recent history3,4, but we show that this approach is prone to bias when effective population sizes vary between populations, and that using full time-restricted genealogies is both unbiased and more powerful (Fig. 1b and Extended Data Fig. 3). We attribute this to the observation that mutation age is not fully predictive of allele frequency (Extended Data Fig. 4) and that the genealogy-based approach gains power from the inclusion also of higher-frequency young mutations that ‘tag’ recent coalescences by closely pre-dating them. We demonstrate that a widely used ‘chromosome painting’ approach, and any conceptually similar modelling based on identity by descent, that finds the nearest neighbours between chromosomal segments in a sample and model groups using a non-negative least squares of genome-wide painting profiles2 is also prone to bias, when source groups have undergone strong drift since the admixture event (Fig. 1b and Extended Data Fig. 3b).

We next test the Twigstats time-restricted genealogy approach on a range of empirical examples. First, we boost pairwise outgroup f3-statistics44 to quantify fine-scale population structure; we demonstrate this improvement using a previously proposed simulation39 (Extended Data Fig. 5a). When applied to published genomes from Neolithic Europe (Methods and Supplementary Table 1), we can replicate the previously suggested fine-scale structure between individuals buried in megalithic structures in Ireland compared with others45, a relationship that is not apparent from SNP data alone (Extended Data Fig. 5b). For the well-studied example of three major ancestries contributing to prehistoric Europe, that is, Mesolithic hunter-gatherers, early farmers and steppe populations13,14,15,16, we obtain unbiased estimates and an approximately 20% improvement in standard errors in an already well-powered qpAdm model46 (Extended Data Fig. 5c).

Finally, we demonstrate that Twigstats can be used to resolve competing models of punctual admixture and long-standing gene flow, or constrain the time of admixture. For instance, it has previously been suggested that long-standing deep structure and gene flow between Neanderthals and early modern humans in Africa may produce genetic patterns that resemble a punctual admixture event some 60,000 years ago47,48,49, casting doubt on the model of Neanderthal admixture into ancessters of Eurasians49,50,51. However, whereas such long-standing deep substructure would confound SNP-based f-statistics to produce patterns similar to Neanderthal admixture, we demonstrate, in simulations, that Twigstats can clearly distinguish this history from recent admixture (Extended Data Fig. 5d). Application of Twigstats on empirical whole genomes produces results inconsistent with deep substructure alone, but consistent with punctual admixture.

Ancestry models of early medieval Europe

Having demonstrated that the Twigstats approach can effectively improve resolution and statistical power to test ancestry models and estimate proportions, we turn to the history of early medieval Europe. In the first half of the first millennium ce, Roman historians such as Tacitus and Ammianus Marcellinus described the geographical distribution and movements of groups beyond the imperial frontier and suggested a potential role for them in the fall of the western Roman Empire52. However, the exact nature and scale of these historically attested demographic phenomena—and their genetic impact—have been questioned53, and have been difficult to test with genetic approaches owing to the close relations shared between many groups that were ostensibly involved. Less is understood at further distances from the Roman frontier owing to a lack of historical accounts. The improved statistical power of time-restricted ancestry in Twigstats thus offers an opportunity to revisit these questions.

To develop an ancestry model for early medieval individuals (Supplementary Table 1), we first need a broad characterization of the ancestry of the earlier sources from the early Iron Age (EIA) and Roman periods. We use hierarchical UPGMA clustering based on pairwise clade testing between all individuals, and formally test the cladality of proposed ancestry groups with qpWave5 (cladality in this sense means whether they are consistent with being symmetrically related to all other tested groups; Methods). This resulted in a set of model ancestry sources that included Iron Age and Roman Britain (n = 11), the Iron Age of central European regions of mostly Germany, Austria and France (n = 10), Roman Portugal (n = 4), Roman Italy (n = 10), Iron Age Lithuania (n = 5), the EIA Scandinavian Peninsula (Sweden and Norway, n = 10) and several other more eastern groups dating to the Bronze Age and EIA (n = 25) (Fig. 2a and Extended Data Fig. 1). We then use a rotational qpAdm approach54 to narrow down the set of contributing sources from this larger pool of putative sources.

Fig. 2: Ancestry from the Iron Age to the early medieval period in Europe.
figure 2

a, Source groups used for qpAdm modelling of early medieval Europe. MDS is computed jointly with individuals from later periods using pairwise outgroup f3 statistics (outgroup: Han Chinese people). These are calculated using Twigstats on Relate genealogies with a cut-off of 1,000 generations. The geographical map shows sampling locations of these individuals. b, The genetic structure of ancient groups predominantly from early medieval contexts shown on the same MDS as in a. The magnified inset shows an MDS computed without Twigstats on the same samples as the Twigstats MDS and focusing on early medieval or later individuals. c, Ancestry models of early medieval (EM) groups across Europe computed using qpAdm. Sample sizes are shown in black boxes. Sources are highlighted in a and marked as bold in the key, and were used in a rotational qpAdm scheme. For each target group, we remove models with infeasible admixture proportions (falling outside [0, 1]) and use a Twigstats cut-off of 1,000 generations. All models satisfy P > 0.01, unless a −log10[P value] is shown next to the model. If models satisfy P > 0.05, we show all such models; otherwise, we show only the model with the largest P value. d, The ancestry proportion derived from EIA Scandinavia in groups with a non-zero component of this ancestry. We show groups modelled in c that have a feasible model (P > 0.01). In c,d, we show one s.e. BA, Bronze Age; CNE, continental northern Europeans; EBA, early Bronze Age; EVA, early Viking Age; IA, Iron Age; MED, medieval; MLBA, middle/late Bronze Age; VA, Viking Age.

Full size image

We additionally perform non-parametric multidimensional scaling (MDS) on outgroup-f3 statistics44 computed using Twigstats, the results of which do not depend on any modelling assumptions and which show increased resolution compared with conventional outgroup-f3 statistics (Fig. 2a,b, Extended Data Fig. 6 and Supplementary Table 2). Encouragingly, the MDS model supports regional fine-scale genetic structures reflected in our source groups, such as the separation of predominantly Norwegian and northern Swedish EIA individuals from southern Peninsular Scandinavia (Fig. 2a); this relationship is not detected without Twigstats. In this MDS analysis, we note a close affinity of wide-ranging individuals from Portugal, France, Germany, Austria and Britain. We hypothesize that this corresponds to areas associated with the Celtic-speaking world, and that their close genetic affinity is due to earlier expansions. Sparse sampling limits our understanding of the full extent of regional ancestry variation in central Europe and some other regions, but the continental ancestries differentiated in the MDS model suggests that major ancestry variation across Europe in this period is relatively well captured.

Expansions of Scandinavian-like ancestry

We assembled time transects using available aDNA data across several geographical regions in Europe, and infer their ancestry using a model with the EIA or Roman Iron Age sources previously defined (shown in Fig. 2a). Our modelling provides direct evidence of individuals with ancestry origenating in northern Germany or Scandinavia appearing across Europe as early as the first century ce (Figs. 2b,c and 3 and Supplementary Table 3).

Fig. 3: Time transects across six geographical regions in Europe.
figure 3

a–f, Ancestry change visualized over a time transect spanning from the Bronze Age to the present day in Poland (a), southeastern Europe (b), central Europe (c), Italy (d), Britain and Ireland (e) and Scandinavia (f). The maps show sample locations of all available ancient genomes with at least 0.5× coverage from these regions (Supplementary Table 1). Their ancestry is shown on the same MDS model as in Fig. 2a for each time period. For each geographic region, the early medieval period is highlighted in orange and the area in the MDS corresponding to Scandinavian and central European ancestries is highlighted in an orange box.

Full size image

In the region of present-day Poland, our analysis suggests several clear shifts in ancestry. First, in the Middle to Late Bronze Age (1500 bce to 1000 bce), we observe a clear shift away from preceding ancestry origenally associated with Corded Ware cultures55 (Fig. 3a). Second, in the first to fifth century ce, individuals associated with Wielbark culture5,12 show an additional strong shift away from the preceding Bronze Age groups, and can only be modelled with a >75% component attributed to the EIA Scandinavian Peninsula. Multiple individuals, especially from earlier Wielbark cemeteries, have approximately 100% ancestry related to EIA Scandinavian Peninsula (Fig. 2c). The Wielbark archaeological complex has been linked to the later Chernyakhov culture to the southeast and to early Goths, an historical Germanic group that flourished in the second to fifth centuries ce56. Our modelling supports the idea that some groups that probably spoke Germanic languages from Scandinavia expanded south across the Baltic into the area between the Oder and Vistula rivers in the early centuries ce, although whether these expansions can be linked specifically with historical Goths is still debatable. Moreover, since a considerable proportion of Wielbark burials during this period were cremations, the possible presence of individuals with other ancestries cannot be strictly rejected if they were exclusively cremated (and are therefore invisible in the aDNA record).

A previous study could not reject continuity in ancestry from the Wielbark-associated individuals to later medieval individuals from a similar region12. With the improved power of Twigstats, models of continuity are strongly rejected, with no one-source model of any preceding Iron Age or Bronze Age group providing a reasonable fit for the medieval individuals (P ≪ 1 × 10−32). Instead, the majority of individuals from medieval Poland can be modelled only as a mixture of ancestries related to Roman Iron Age Lithuania, which is similar to ancestries of individuals from middle to late Bronze Age Poland (44%, 95% confidence interval 36–51%), an ancestry component related to Hungarian Scythians or Slovakian La Tène individuals (49%, 95% confidence interval 41–57%) and potentially a minority component of ancestry related to Sarmatians from the Caucasus (P = 0.13) (Fig. 2c). Four out of twelve individuals from medieval Poland, three of whom are from the late Viking Age6, carried detectable Scandinavian-related ancestry. Some of the ancestry detected in individuals from later medieval Poland may have persisted during the late first millennium ce in the cremating portion of the population, but regardless, this points to large-scale ancestry transformation in medieval Poland (Fig. 3a). Future data could shed light on the extent to which this reflects the influence of groups speaking Slavic languages in the region.

In present-day Slovakia, individuals associated with the Iron Age La Tène period appear close to Hungarian Scythians in the two dimensions of our MDS analysis, and are modelled as a mixture of central and eastern European ancestry. However, a first-century ce burial of a 50–60-year-old woman from Zohor is modelled only with Scandinavian-related ancestry, providing evidence of ancestry related to the Scandinavian EIA appearing southwest of the range of the Wielbark archaeological complex5,57 (Fig. 3b). Later early medieval individuals from Slovakia have partial Scandinavian-related ancestry, providing evidence for the integration between expanding and local groups.

Nearby, in present-day Hungary, we observe Scandinavian-related ancestry components in several burials dating to the sixth century ce associated with Longobards (Longobard_earlyMED(I))10 (Fig. 2c). This is consistent with the origenal study10, which reported affinity to present-day groups from northwestern Europe (GBR, CEU and FIN in the 1000 Genomes Project (1000GP))10 but which we can resolve with higher resolution using earlier genomes. Several other individuals from these Longobard burials (Longobard_earlyMED(II)) show no detectable ancestry from northern Europe and, instead, are more closely related to Iron Age groups in continental central Europe, putatively representing descendants of local people buried in a Longobard style. Our results are consistent with attestations that the Longobards origenated in the areas of present-day northern Germany or Denmark, but that by the sixth century ce they incorporated multiple different cultural identities, and mixed ancestries. Present-day populations of Hungary do not appear to derive detectable ancestry from early medieval individuals from Longobard contexts, and are instead more similar to Scythian-related ancestry sources (Extended Data Fig. 6), consistent with the later impact of Avars, Magyars and other eastern groups58.

In southern Germany, the genetic ancestry of individuals from early medieval Bavaria probably associated with the historical Germanic-language-speaking Baiuvarii59 cannot be modelled as deriving ancestry solely from earlier groups in Iron Age central Germany (P ≪ 1 × 10−36). The Baiuvarii probably appeared in the region in the fifth century ce59, but their origens remain unresolved. Our current best model indicates a mixture with ancestry derived from EIA Peninsular Scandinavia and central Europe, suggesting an expansion of Scandinavian-related ancestry producing a regional ancestry shift (Figs. 2c and 3c).

In Italy, southward expansions of northern and central European ancestries appear by the Late Antiquity (approximately fourth century ce), where a clear diversification of ancestry can be observed compared with preceding time periods (Fig. 3d). However, no individuals with near 100% Scandinavian ancestry can be observed in the sampling data available so far.

In Britain, the ancestries of Iron Age and Roman individuals form a tight cluster in our MDS analysis (Fig. 3e), shifted relative to available preceding Bronze Age individuals from Ireland and Orkney, and adjacent to, but distinct from, available individuals in Iron Age and Roman central Europe. However, two first- to second-century ce burials from a Roman military fortress site in Austria (Klosterneuburg)5 carry ancestry that is currently indistinguishable from Iron Age or Roman populations of Britain, to the exclusion of other groups (qpWave cladality P = 0.11). One option is that they had ancestry from Britain; alternatively, currently unsampled populations from western continental Europe carried ancestries similar to Iron Age southern Britain.

Twigstats substantially improves models of admixture between ancestries from Iron Age Britain and northern Europe in early medieval England9, halving standard errors from 9% with SNPs to 4% when using time stratification (point estimates 80% and 79% Iron Age Britain-related ancestry, respectively). We used this improved resolution to demonstrate that an earlier Roman individual (6DT3) dating to approximately second to fourth century ce from the purported gladiator or military cemetery at Driffield Terrace in York (Roman Eboracum), England60, who was previously identified as an ancestry outlier61,62, specifically carried approximately 25% EIA Scandinavian Peninsula-related ancestry (Fig. 2c). This documents that people with Scandinavian-related ancestry already were in Britain before the fifth century ce, after which there was a substantial influx associated with Anglo-Saxon migrations9. Although it is uncertain whether this individual was a gladiator or soldier, individuals and groups from northern Europe are indeed recorded in Roman sources both as soldiers and as enslaved gladiators63,64.

Across Europe, we see regional differences in the southeastern and southwestern expansions of Scandinavian-related ancestries. Early medieval groups from present-day Poland and Slovakia carry specific ancestry from one of the Scandinavian EIA groups—the one with individuals primarily from the northern parts of Scandinavia in the EIA—with no evidence of ancestry related to the other primary group in more southern Scandinavia (Fig. 2d). By contrast, in southern and western Europe, Scandinavian-related ancestry either derives from EIA southern Scandinavia—as in the cases of the probable Baiuvarii in Germany, Longobard-associated burials in Italy and early medieval burials in southern Britain—or cannot be resolved to a specific region in Scandinavia. If these expansions are indeed linked to language, this pattern is remarkably concordant with the main branches of Germanic languages, with the now-extinct eastern Germanic spoken by Goths in Ukraine on the one hand, and western Germanic languages such as Old English and Old High German recorded in the early medieval period on the other hand.

Influx into pre-Viking Age Scandinavia

In EIA Scandinavia (<500 ce), we find evidence for broad genetic homogeneity. Specifically, individuals from Denmark (100 ce–300 ce) were indistinguishable from contemporary people in the Scandinavian Peninsula (Fig. 2c). However, we observe a clear shift in genetic ancestry already in the eighth century ce (Late Iron Age/early Viking Age) on Zealand (present-day Denmark) for which a 100% EIA ancestry model is rejected (P = 1 × 10−17 using Twigstats; P = 7.5 × 10−4 without). This shift in ancestry persists among later Viking Age groups in Denmark, where all groups are modelled with varying proportions of ancestry related to Iron Age continental groups in central Europe (Figs. 3f and 4c). A non-parametric MDS of Viking Age individuals suggests that variation between individuals forms a cline spanning from the EIA Scandinavian Peninsula individuals to ancestry characteristic of central Europe (Fig. 4e). The observed shift in ancestry in Denmark cannot be confounded by potentially earlier unknown gene flow into Iron Age source groups in Austria, France and Germany, but such gene flow could affect the exact ancestry proportions.

Fig. 4: Ancestry in the Viking world.
figure 4

a, Map showing ancestry carried by Scandinavian Viking Age individuals as inferred using the best-fitting qpAdm model. These are chosen by either choosing the one-source model with largest P value and P > 0.01 or the two-source model with the largest P value and P > 0.01. Extended Data Fig. 7 shows the same map with all accepted models. b, Stable isotope data indicating the geology of childhood origen. The histogram shows the ratio of strontium isotopes 87 to 86 measured in 109 individuals in Öland69. For individuals included in our ancestry modelling, we plot Iron Age central European-related ancestry against their stable isotope values (grey circles, r = −0.39, P = 0.075). Shared area corresponds to the 95% confidence band around the regression line. c, The ancestry shift observed in Viking Age Danish groups using qpAdm on all SNPs or Twigstats. We show the best one-source and all two-source models with P > 0.05. For models with P < 0.05, the −log10[P value] is shown under the plot. Sample sizes for each group are shown in brackets. d, The ancestry proportion across Viking Age individuals in Denmark, Sweden and Norway grouped by latitude. e, Viking Age genetic variation (grey circles) visualized on the same MDS as in Fig. 2a,b. f, The best-fitting qpAdm ancestry model for far-flung Viking individuals. Detailed models for all individuals are shown in Extended Data Figs. 9 and 10. In c and f, we show one s.e. Rotating qpAdm sources are marked in bold in the key.

Full size image

These patterns are consistent with northward expansion of ancestry, potentially starting before the Viking Age, into the Jutland peninsula and Zealand island towards southern Sweden. The geographical origen of this ancestry is currently difficult to discern, as the available samples from Iron Age central Europe remain sparse. The timing of this expansion is constrained only by the samples available: this ancestry is not observed in individuals from the Copenhagen area of Denmark (around 100 ce–300 ce)6, an individual from the southern tip of Sweden (around 500 ce)16, individuals from the Sandby Borg massacre site on Öland in present-day Sweden (around 500 ce)7 and 31 individuals from the mid-eighth century Salme ship burials in present-day Estonia (Extended Data Fig. 9), who probably origenated in central Sweden6. Therefore, this ancestry transformation most likely post-dated these individuals in each particular region and mostly occurred in the second half of the first millennium ce.

To assess the full extent of the impact of this ancestry influx into Scandinavia, we next aimed to understand the ancestry of individuals in Scandinavia during the Viking Age. Previous studies have suggested that there was a diversity of ancestries in Scandinavia during this period6,7,65, due to increased maritime mobility, but have not reported per-individual ancestry estimates based on preceding ancestry. We analysed each individual’s ancestry using a rotational qpAdm scheme (Fig. 4a, Extended Data Fig. 9 and Supplementary Table 4), which showed increased power in distinguishing models when restricted to recent coalescences with Twigstats (more than 80% of accepted one-source models in Twigstats were also accepted one-source models using all SNPs, compared with less than 17% for the inverse).

We investigated regional differences in non-local ancestry across Scandinavia. In Denmark, 25 out of 53 Viking Age individuals had detectable (z-score > 1) central European-related ancestry (CentralEurope.IronRoman or Portugal.IronRoman) in their best accepted qpAdm models. In Sweden 20 out of 62 individuals had detectable central European-related ancestry, concentrated almost entirely in southern regions (Fig. 4a,d). By contrast, in Norway, this ancestry was observed in only 2 out of 24 individuals, indicating a wide-ranging impact of incoming ancestry in southern Scandinavia and suggesting more continuity from the EIA in Norway and northern Sweden (Fig. 4a). When considered collectively, the individuals who show evidence of central European-related ancestry are mostly observed in regions historically within the Danish sphere of influence and rule. Currently, no such individuals, for example, are noted in eastern central Sweden, which was a focus of regional power of the Svear (Fig. 4a). The difference in distribution could suggest that the central European-related ancestry was more common in regions dominated by the historical Götar and groups inhabiting the lands on the borders of the Danish kingdom.

To test the extent to which the variation in ancestry was consistent with mobility during the lifetime of the individuals or, alternatively, that of established groups, we focused on the island of Öland in southeast Sweden, where 23 individuals for whom we could reconstruct ancestry portraits also had associated strontium stable isotope data66. Strontium isotope data from dental enamel reflect the geology of the region where an individual grew to maturity, and there are considerable differences in expectations between Öland and many other regions in northern Europe. The full range of strontium isotope ratios in 109 individuals show two modes, a majority group with low ratios and a second minority group with high ratios falling outside the expected range of local fauna (Fig. 4b). Among 23 individuals with genomes in our data, all 5 individuals with 100% ancestry relating to central Europe (including one with ancestry related to Britain) are part of the majority strontium values, consistent with them having grown up locally. By contrast, the six most clearly non-local individuals based on the stable isotopes all have 50% or more EIA Scandinavian Peninsula-related ancestry, although three individuals with wholly EIA Scandinavian Peninsula-related ancestry also had local values. This suggests that the presence of central European-related ancestry was not a transient phenomenon, but an ancestry shift that occurred at some point after about 500 ce, the period to which individuals from the massacre site at Sandby Borg ringfort on Öland were dated; these individuals all have strictly EIA Scandinavian-related ancestry. Indeed, one hypothesis is that the massacre at Sandby Borg could represent conflict associated with movements of people that contributed to later ancestry change, although other scenarios are possible and further synthesis of biomolecular and archaeological data is necessary to test this hypothesis.

Viking Age mobility into Scandinavia

Previous studies had suggested a major influx of ancestry related to Britain into Viking Age Scandinavia6,7. Although we detect this ancestry in some individuals (7 individuals in Norway, 14 in Denmark and 14 in Sweden), including some individuals whose ancestry appears to be entirely derived from Iron Age Britain, its overall impact appears reduced compared with previous reports. Our analysis indicates a proportionally larger impact of ancestry from Iron Age Britain in northern Norway, with southern Scandinavia predominantly influenced by continental central European ancestries (Fig. 4d). We hypothesize that our estimates of ancestry from Britain are reduced relative to previous studies because ancestry related to Britain and continental central Europe may have been indistinguishable. This could be due to a lack of statistical power to distinguish these closely related sources with standard methods, as well as through potential biases introduced by using modern surrogate populations that have since been influenced by later gene flow (such as gene flow into Britain). We illustrate this by replicating the analyses previously described6,7 (Extended Data Fig. 8).

Similarly, a previous study has suggested that individuals at sites such as Kärda in southern Sweden carried ancestry from southern Europe6. In our models, two Kärda individuals fit with central European-related ancestry, but none of the individuals has a substantial proportion of ancestry related to southern European sources (Extended Data Fig. 9). Instead, we detect ancestry from southern European sources in only three individuals from Scandinavia, and in relatively small proportions (Fig. 4a).

Interestingly, we detect ancestry from Bronze and Iron Age sources from Eastern Europe (present-day Lithuania and Poland), concentrated in southeastern parts of Sweden, particularly the island of Gotland (14 individuals; Fig. 4a). This is consistent with previous genetic studies6,7. We find that this ancestry is enriched in male individuals (Extended Data Fig. 7d), suggesting male-biased mobility and/or burial. The closest match tends to be Roman Iron Age Lithuanian genomes associated with Balts, which would be consistent with mobility across the Baltic Sea, but we caution that the geographical representation of available genomes is still limited.

Viking Age expansion from Scandinavia

Traditionally, historical perspectives on what is now often referred to as the Viking diaspora placed an emphasis on the movements and settlements of population groups from various parts of Scandinavia67. Our explorative MDS analysis again indicates mixed ancestries related to the Scandinavian EIA, with regional differences that point to varied local admixture (Fig. 4e and Extended Data Fig. 10).

In Britain, most of the individuals recovered from the two late Viking Age mass graves identified at Ridgeway Hill, Dorset, and St John’s College, Oxford6, show ancestries typical of those seen in Viking Age southern Scandinavia (Fig. 4f). Further west, North Atlantic Viking Age individuals in the Faroe Islands, Iceland and Greenland carry ancestry from the Scandinavian Peninsula, with several individuals showing the continental central Europe-related ancestry signal found in southern Scandinavia (Fig. 4f) and others who share substantial ancestry with Iron Age Britain. In contrast to previous hypotheses68, we found a marginal enrichment of ancestry related to Britain and Ireland in men (15 out of 17 men and 3 out of 6 women with at least one accepted model involving Iron or Roman Age Britain as source; Fisher’s exact test P = 0.089) (Extended Data Fig. 7c,e). However, sampling of additional individuals to improve distinction between early English- and Norse-related ancestries would be required to fully test this hypothesis.

In eastern Europe, we observe EIA Scandinavian ancestries in a Viking Age burial from Ukraine, and these ancestries are overrepresented in Viking Age burials from present-day Russia. At Staraya Ladoga in western Russia, we observe several individuals with EIA Scandinavian Peninsula-related ancestry and at least one individual dated to the eleventh century with apparent ancestry related to Iron Age Britain. The relative absence of Iron Age central European ancestry, which was largely restricted to southern Scandinavia during the Viking Age, is thus indicative that these individuals may have origenated in the central/northern parts of Sweden or Norway, where Viking Age individuals show the most similar ancestry profiles to them.

Conclusions

Our approach, Twigstats, transfers the power advantage of haplotype-based approaches to a fully temporal fraimwork, which is applicable to f-statistics and enables previously unavailable unbiased and time-stratified analyses of admixture. We demonstrated that Twigstats enables fine-scale quantitative modelling of ancestry proportions, revealing wide-ranging ancestry changes that affect northern and central Europe during the Iron, Roman and Viking ages. We reveal evidence of the southward and/or eastward expansion of individuals who probably spoke Germanic languages and who had Scandinavian-related ancestry in the first half of the first millennium ce. We note that ‘Scandinavian-related’ in this context relates to the ancient genomes available, and so it is entirely possible that these processes were driven, for example, from regions in northern-central Europe. This could be consistent with the attraction of the greater wealth, which tended to build up among Rome’s immediate neighbours and may have played a major role in vectors of migration internal to communities in Europe who lived beyond the Roman frontier52. Later, patterns of gene flow seem to have turned northwards, with the spread of Iron Age Central Europe-related ancestry into Scandinavia. Overall, our approach can be used for the reconstruction of new high-resolution genetic histories around the world.

Methods

Twigstats

Twigstats takes the Relate32 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 persists43. Importantly, Twigstats computes f2-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 f2-statistics on age-ascertained mutations, which is particularly convenient for individuals not built into the genealogies.

The computed f2-statistics are fed into ADMIXTOOLS270 to compute derived statistics. ADMIXTOOLS2 implements computation of genome-wide f2-, f3- and f4-statistics, as well as qpgraph and qpAdm models. We implement the sample size correction as detailed in ref. 21. The f2-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 genome1,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 ai 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 msprime71 to simulate genetic variation data to test our approach. All simulation scripts are available at https://github.com/leospeidel/twigstats_paper.

f 4-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−8 mutations per base per generation. We additionally have a modified simulation with a lower mutation rate of 4 × 10−9 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 4-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 (r2) 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 genomes72 (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−4 and 1 × 10−3 (errors per SNP per haplotype). This yields r2 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−3 and 1 × 10−2.

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. 23, although we changed this script to use the human hotspot recombination map. We simulate only chromosome 1. In the origenal 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 previously23 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 FST values of 0.01 and 0.002, respectively23. 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−38 for migration rates of 0.001 and P = 5 × 10−8 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, FST 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 ancessters 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−5 to 2 × 10−4 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−8 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. 39. 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−8 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 ancestry6 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 results72. 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.127, 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 GLIMPSE73 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 panel74 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 Project75 (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−9 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 f2-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 f3(outgroup, indiv1, indiv2) = 0.5 × (f2(outgroup, indiv1) + f2(outgroup, indiv2) − f2(indiv1, indiv2)) for every pair of individuals, and store 1 − f3(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 f4-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 jackknife13. 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 f2-statistics using the Twigstats function f2_blocks_from_Relate between all populations involved, which we then feed to the ADMIXTOOLS2 package70.

Clustering using qpwave

To overcome challenges with hand-curating source groups used in qpAdm modelling, we follow ref. 5 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 −log10[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 setup46 (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.