Article
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 simula-
tion in Extended Data Fig. 5d emulates Neanderthal admixture, in which 
Neanderthals and ancestors of modern humans split 25,000 genera-
tions 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 sub-
group 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 muta-
tion 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 muta-
tion 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 pro-
vided by refs. 7,8 were only available aligned against the GRCh38 refer-
ence 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 rese-
quenced, choosing the file with the highest coverage. We then filter