The Yamnaya Steppe Herders from Russia

MATERIALS AND METHODS
Processing sites for the newly reported individuals
Most (186 of 206, 90.2%) of the newly reported individuals were entirely processed at the Max Planck Institute for the Science of Human History in Jena, Germany, and the full details of sampling and ancient DNA wet laboratory work and bioinformatic processing are summarized in what follows. The individuals from the site Makotřasy were initially sampled and processed into powder at the University of Vienna, followed by subsequent laboratory work and bioinformatic and ancient DNA analysis at Harvard Medical School following previously described protocols (61).

Sampling
In total, 389 pars petrosa, teeth, and bones from 261 individuals were processed as part of this study. Upon introduction into the clean room facilities at the Max Planck Institute for the Science of Human History in Jena, Germany, all samples were wiped with 5% bleach and ultraviolet (UV) irradiated for 20 min on each side. Teeth were sampled by removing the crown followed by drilling into the pulp chamber to create bone powder. Pars petrosa were sampled by drilling into their dense region (78) to create bone powder. Between 50 and 100 mg of resulting bone powder from each sample were collected in different 2-ml Biopure tubes (one tube per sample) and used in subsequent DNA extraction.

DNA extraction
One milliliter of extraction buffer [containing 0.9 ml of 0.5 M EDTA, 0.025 ml of proteinase K (0.25 mg/ml), and 0.075 ml of UV high-performance liquid chromatography (HPLC) water] was added to Biopure tubes containing bone powder. Biopure tubes were then sealed with Parafilm and incubated overnight on a rotating wheel at 37°C. After incubation, Biopure tubes were spun for 2 min at 18,500 relative centrifugal force (rcf), separating the soluble from insoluble parts of the resulting solution. The soluble part was transferred to a 50-ml falcon tube containing 10 ml of binding buffer and 400 μl of sodium acetate (3 M, pH 5.2). Resulting mixture was transferred to a high pure extender assembly (HPEA) falcon tube, which was centrifuged at 1500 rpm in a 50-ml Thermo Fisher Scientific TX-400 Swinging Bucket Rotor for 8 min. The column from each HPEA tube was removed and inserted into a fresh collection tube and centrifuged at 18,500 rcf for 2 min. A total of 450 μl of wash buffer from the high pure viral nucleic acid kit (HPVNAK) was added to each column, which was then centrifuged at 8000 rcf for 1 min. Columns were then removed and placed into new collection tubes. Another round of washing was performed whereby 450 μl of wash buffer from the HPVNAK was added to each column and centrifuged at 8000 rcf for 1 min. Columns containing washed DNA were then transferred to 1.5-ml siliconized tubes. Fifty microliters of TET (tris-EDTA + Tween 20) buffer were added to the center of columns, and the columns were then incubated at room temperature for 3 min and centrifuged at 18,500 rcf for 1 min. Another 50 μl of TET buffer was added to the center of columns, after which they were centrifuged once more at 18,500 rcf for 1 min. The resulting 100 μl of DNA extracts was stored at −20°C until further processing.

DNA libraries and in-solution capture
Twenty-five microliters of DNA extract was used for the construction of (in most cases) double-stranded uracil DNA glycosylase (UDG) half-treated DNA libraries. UDG repair was performed by adding DNA extract to a 25-μl mastermix containing 6 μl of 10× Buffer Tango, 6 μl of 10 mM adenosine 5′-triphosphate, 0.5 μl of bovine serum albumin (BSA; 20 mg/ml), 0.2 μl of 25 mM each deoxynucleoside triphosphate (dNTP), 3.6 μl of 1-U USER enzyme, and 8.7 μl of UV HPLC water. Resulting mixture was incubated at 37°C for 30 min followed by 12°C for 1 min. Inhibition of UDG treatment was achieved through the addition of 3.6 μl of 2-U uracil glycosylase inhibitor (UGI) to each tube followed by incubation at 37°C for 30 min and again at 12°C for 1 min. Blunt-end repair was performed by adding 3 μl of 10-U T4 polynucleotide kinase and 1.65 μl of 3-U T4 DNA polymerase and incubating the resulting mixture at 25°C for 20 min and then at 12°C for 10 min. Blunt-end repaired mixture was purified using a MinElute kit and eluted in 20 μl of elution buffer (EB) containing 0.05% Tween. Illumina adapters were ligated onto DNA molecules through the mixture of 18 μl of eluate from the previous step with 20 μl of 2× Quick Ligase Buffer, 1 μl of 10 μM Adapter Mix, and 1 μl of 5-U Quick Ligase. Resulting mixture was incubated at 22°C for 20 min followed by purification with a MinElute kit and elution in 22 μl of EB containing 0.05% Tween. Adapter fill in reaction was performed by adding 20 μl of eluate from the previous step to 4 μl of 10× isothermal buffer, 0.2 μl of 25 mM each dNTP, 2 μl of 8-U Bst 2.0 polymerase, and 13.8 μl of UV HPLC water, followed by incubation at 37°C for 30 min and 80°C for 10 min. Resulting libraries were stored at −20°C until further processing. Unique library-specific indexes were added to the 5′ and 3′ ends of molecules in each library through an indexing polymerase chain reaction (PCR). Each library was split into four separate indexing PCR reactions, which were carried out using 10 μl of 10× Pfu Turbo buffer, 1.5 μl of BSA (20 mg/ml), 1 μl of 25 mM each dNTP, 1 μl of 2.5-U Pfu Turbo polymerase, 73.5 μl of UV HPLC water, 2 μl of 10 μM P5 index, 2 μl of 10 μM P7 index, and 9 μl of DNA library. Amplification was achieved through an initial denaturation at 95°C for 2 min, followed by 10 cycles of 95°C for 30 s, 58°C for 30 s, and 72°C for 1 min, followed by 72°C for 10 min. Resulting indexed libraries of the same sample were pooled and purified using a MinElute purification kit. Purified libraries were quantified using quantitative PCR and amplified to 1013 copies. Amplified libraries were shallow shotgun sequenced (~5 million reads) on an Illumina HiSeq or NextSeq platform to estimate the general human DNA content, presence of ancient DNA damage, and mitochondrial:nuclear coverage ratio. Libraries with >0.1% endogenous human DNA and >5% C-to-T misincorporations at the 5′ end were chosen for 1240k capture (8) and mitochondrial DNA (mtDNA) capture. In cases where more than one library from the same individual satisfied the criteria for capture, the better quality (higher endogenous DNA content) library was used for 1240k and mtDNA capture.

Sequencing
Postcapture libraries were single-end (75 cycles) or paired-end sequenced (2 × 50 or 2 × 75) on HiSeq or NextSeq Illumina platforms to a depth of 20 to 50 million reads per library. Resulting sequence data were processed through EAGER (v1.92.38) (79). Illumina adapters were removed using AdapterRemoval (v2.2.0) (80), and in case of paired-end sequencing, corresponding reads from the same template molecule with a minimum of 11 base pairs (bp) of overlap were merged. Fastq files of merged and unmerged reads were concatenated, and reads shorter than 30 bp were discarded. Processed reads were mapped to the human reference genome (hg19) using BWA-aln and BWA-samse (v0.7.12) (81) applying maxdiff (-n) 0.01 and seeding turned off (-l 10000). Resulting bam files were sorted, and duplicate reads were removed using DeDup (v0.12.1) (https://github.com/apeltzer/DeDup). Damageprofiler (v0.3.10) (https://github.com/Integrative-Transcriptomics/DamageProfiler) was used to calculate rates of misincorporation in read termini of DNA fragments in our captured libraries. BAM files had the last three bases from both 5′ and 3′ ends of reads and corresponding base quality scores masked for downstream analyses (82).

Sex determination and authentication
The genetic sex of each sample (bam file) was determined by calculating the normalized mean coverage on the X (mean X coverage/mean autosome coverage) and Y (mean Y coverage/mean autosomal coverage) chromosomes (83). Samples with normalized mean Y coverage values greater than 0.2 were assigned male. Contamination was estimated in males by calculating the rate of heterozygosity on their X chromosome (84). In addition, we used Schmutzi to estimate the mitochondrial contamination in all libraries (85). Schmutzi was run on BAM files resulting from mapping 1240K capture sequencing data to the human mitochondrial reference genome. In cases where 1240k data were not enough to give an mtDNA contamination, we ran Schmutzi on the mtDNA capture data mapped to the human mitochondrial reference genome.

Genotyping
We used SAMtools (v1.3) (86) mpileup and pileupCaller from the sequenceTools (v1.4.0.2) package (https://github.com/stschiff/sequenceTools) to call pseudo diploid genotypes by sampling a random high-quality allele (base quality, ≥30; mapping quality, ≥30) from each of the 1240k sites (8). Newly generated genotype data for this study were then merged to a compiled dataset of previously published ancient and modern worldwide populations (https://reich.hms.harvard.edu/allen-ancient-dna-resource-aadr-downloadable-genotypes-present-day-and-ancient-dna-data) (v42.2) using mergeit (v2450) from the EIGENSOFT package (https://github.com/DReichLab/EIG).

Mitochondrial and Y chromosome haplogroups
Mitochondrial haplogroups were called by mapping 1240k or mtDNA capture data to the human mitochondrial reference genome followed by creating pileups at each position (map quality and base quality filter of 30) and calling the most frequent base at each position. Resulting genotype information was converted to fasta files, and haplogroups were called using haplofind (87). Y chromosome haplogroups were called by mapping 1240k capture data to the whole human reference genome (hg19) followed by visual inspection of ancestral/derived alleles (after map quality and base quality filter of 30) at ISOGG (v15.58 April 2020) sites.

Principal components analysis
PCA was conducted using smartpca (v1600) from the Eigensoft package (https://github.com/DReichLab/EIG). Principal components were calculated on the genotype data of modern West Eurasian individuals (table S7) (88–90). Ancient individuals were projected (lsqproject: YES) onto the axes calculated from modern individuals. “shrinkmode: YES” was used to account for artificial stretching of principal component axes between projected (ancient) individuals and modern individuals. Fst values were also calculated in smartpca using “fsthiprecision: YES” and “inbreed: YES” parameters.

Ancestry decomposition and admixture modeling
F statistics, qpWave, and qpAdm runs were conducted in admixr v0.7.1 (91), a wrapper program around ADMIXTOOLS (88). Selection of outgroups for each analysis is indicated in the corresponding Supplementary Table.
Linear modeling of pre-CW HG ancestry (Fig. 3A) was performed using segmented linear regression as implemented in R using the segmented function for v.1.2-0 of the segmented library (92). To select the optimal number of breakpoints, we compared Akaike’s information criterion (AIC) for models with between zero and four breakpoints. The AIC is a score that considers how well a model fits the data while simultaneously penalizing models with additional parameters. In this way, model fit must be significantly improved for a more complicated model with additional parameters to be accepted over a simpler, nested model. A linear regression model with one breakpoint was found to have the minimum AIC and hence was selected (93).
DATES v753 (61) was used to estimate length distributions of ancestry tracts and infer admixture dates between Anatolia_Neolithic and WHG (here, Loschbour+Körös_HG+Germany_BDB). Parameters binsize 0.001, maxdis 1, seed 77, jackknife YES, qbin 10, runfit YES, afffit YES, lovalfit 0.45, minparentcount 1, and checkmap YES were used.

Y haplogroup frequency simulations
To investigate the process of Y-haplogroup inheritance in early and late CW groups, we simulated 106 realizations, assuming a generational time of 25 years, and analyzed the results using approximate Bayesian computation (ABC). For the ith realization, we assumed a constant population size of Ni ~ U(102,104), with a starting a proportion of R1a-M417(xZ645) of pi ~ TN(0, 1)(0.27,0.134) from a truncated normal distribution based on the observed proportion of R1a-M417(xZ645) of 3 of 11 (0.27). For each simulation, we also included a selection coefficient denoted si ~ U(−1,1). Under random mating, for generation j + 1, let the number of a male offspring carrying R1a-M417(xZ645) be Xij+1 ~ B(Ni, wj), where wj = Xj−1/Nj−1. However, if one includes a selection coefficient, then wj = min(1,(1 + si)Xj−1/Nj−1). Hence, one may interpret sj as the average increase in the proportion of male offspring that R1a-M417(xZ645) individuals were having over this period. We then compared our observed number of per-generation R1a-M417(xZ645) to our simulated realizations using the rejection method and keeping the top closest 0.05% realizations (selected via cross-validation) to form samples from the joint posterior distributions for our simulation parameters. All ABC and cross-validation analyses were performed in R using the abc package (94).

Supplementary Materials for
Dynamic changes in genomic and social structures in third millennium
BCE central Europe

www.science.org/action/downloadSupplement?doi=10.1126%2Fsciadv.abi6941&file=sciadv.abi6941_SM.pdf