remove any duplicate individuals, such as individuals who were rese-
quenced, 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 popula-
tions (CEU, Utah residents with northern and western European ances-
try; CHB, Han Chinese in Bejing, China; FIN, Finnish in Finland; GBR, 
British in England and Scotland; BS, Iberian populations in Spain; TSI, 
T oscani 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 gen-
eration 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 indi-
viduals, 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 itera-
tively selected a subset of source groups and used all remaining puta-
tive 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 gene-
alogies of SGDP modern and imputed ancient genomes in our qpAdm 
modelling and first compute f2-statistics using the Twigstats function