Bronze Age population dynamics and the rise of dairy pastoralism on the eastern Eurasian steppe

We built and sequenced uracil-DNA-glycosylase–half (19), double-indexed Illumina libraries for genomic DNA extracted from teeth or femora from DSKC-associated burials in Khövsgöl, Mongolia. Twenty of 22 libraries exhibited good human DNA preservation, with a mean host endogenous content of 14.9% (range 0.2–70.0%); two libraries yielded very little human DNA (<0.05%) and were excluded from further analysis (SI Appendix, Table S2). Libraries were then enriched for 1.2 million variable sites in the human genome (1240K) using in-solution hybridization (2, 3). All individuals (12 males, 8 females) showed characteristic patterns of chemical modifications typical of ancient DNA (SI Appendix, Fig. S5), and 18 individuals yielded both low estimates of modern DNA contamination (≤1% mitochondrial and nuclear contamination) and sufficient genome coverage for subsequent analysis (0.11× to 4.87× mean coverage for target sites) (SI Appendix, Table S3). No close relative pairs were identified among the ancient individuals (SI Appendix, Fig. S6). Two individuals with high endogenous content on screening (ARS008, 70.0%; ARS026, 47.6%) were deeply sequenced to obtain whole genomes (∼3.3× coverage) (SI Appendix, Table S3). We intersected our ancient data with a published world-wide set of ancient and contemporary individuals (Dataset S1) whose genotypes are determined for 593,124 autosomal SNPs on the Affymetrix HumanOrigins 1 array (20).

To characterize the genetic profile of DSKC-associated LBA Khövsgöl individuals (Khövsgöls), we performed principal component analysis (PCA) of Eurasian populations (SI Appendix, Fig. S7). PC1 separates eastern and western Eurasian populations, with central and north Eurasian populations falling in an intermediate position (SI Appendix, Fig. S7). PC2 separates eastern Eurasian populations along a north–south cline, with northern Siberian Nganasans and the Ami and Atayal from Taiwan forming the northern and southern end points, respectively. Most LBA Khövsgöls are projected on top of modern Tuvinians or Altaians, who reside in neighboring regions. In comparison with other ancient individuals, they are also close to but slightly displaced from temporally earlier Neolithic and Early Bronze Age (EBA) populations from the Shamanka II cemetry (Shamanka_EN and Shamanka_EBA, respectively) from the Lake Baikal region (SI Appendix, Fig. S7) (4, 21). However, when Native Americans are added to PC calculation, we observe that LBA Khövsgöls are displaced from modern neighbors toward Native Americans along PC2, occupying a space not overlapping with any contemporary population (Fig. 2A and SI Appendix, Fig. S8). Such an upward shift on PC2 is also observed in the ancient Baikal populations from the Neolithic to EBA and in the Bronze Age individuals from the Altai associated with Okunevo and Karasuk cultures (1). These observations are consistent with LBA Khövsgöls and other ancient Siberians sharing more ancestry with Native American-related gene pools than modern populations in the region do.

Notably, two individuals fall on the PC space markedly separated from the others: ARS017 is placed close to ancient and modern northeast Asians, such as early Neolithic individuals from the Devil’s Gate archaeological site (22) and present-day Nivhs from the Russian far east, while ARS026 falls midway between the main cluster and western Eurasians (Fig. 2A). Genetic clustering with ADMIXTURE (23) further supports these patterns (Fig. 2B and SI Appendix, Fig. S9). We quantified the genetic heterogeneity between Khövsgöl individuals by calculating f₄ symmetry statistics (24) in the form of f₄(chimpanzee, outgroup; Khövsgöl₁, Khövsgöl₂) for all pairs against 18 outgroups representative of world-wide ancestries (SI Appendix, Fig. S10). As expected, the two outliers did not form a clade with the rest of individuals and therefore we treated each individual separately in subsequent analyses. For the remaining 16 individuals, 14 were merged into a single main cluster based on their minimal genetic heterogeneity. The other two individuals (ARS009 and ARS015) were excluded from this cluster because they broke symmetry with four and two individuals (maximum |Z| = 3.9 and 4.7 SE), respectively, and were also slightly displaced from the others in our PCA (Fig. 2A).

Next, we quantified the genetic affinity between our Khövsgöl clusters and world-wide populations by calculating outgroup-f₃ statistics with Central African Mbuti as an outgroup (25). For the main cluster, top signals were observed with earlier ancient populations from the Baikal region, such as the early Neolithic and EBA individuals from the Shamanka II cemetry (4), followed by present-day Siberian and northeast Asian populations, such as Negidals from the Amur River basin and Nganasans from the Taimyr peninsula (Fig. 3A and SI Appendix, Fig. S11 A and B). As expected based on their nonoverlapping positions on PCA, however, Khövsgöls do not form a cluster with these high-affinity groups, as shown by f₄ symmetry tests in the form of f₄(Mbuti, X; Siberian, Khövsgöls). Interestingly, Upper Paleolithic Siberians from nearby Afontova Gora and Mal’ta archaeological sites (AG3 and MA-1, respectively) (25, 26) have the highest extra affinity with the main cluster compared with other groups, including the eastern outlier ARS017, the early Neolithic Shamanka_EN, and present-day Nganasans and Tuvinians (Z > 6.7 SE for AG3) (red shades in Fig. 3B and SI Appendix, Fig. S11 C and D). This extra affinity with so-called “Ancient North Eurasian” (ANE) ancestry (27) may explain their attraction toward Native Americans in PCA, because Native Americans are known to have high proportion of ANE ancestry (20, 25). Main-cluster Khövsgöl individuals mostly belong to Siberian mitochondrial (A, B, C, D, and G) and Y (all Q1a but one N1c1a) haplogroups (SI Appendix, Table S4).

Previous studies show a close genetic relationship between WSH populations and ANE ancestry, as Yamnaya and Afanasievo are modeled as a roughly equal mixture of early Holocene Iranian/Caucasus ancestry (IRC) and Mesolithic Eastern European hunter-gatherers, the latter of which derive a large fraction of their ancestry from ANE (20, 28). It is therefore important to pinpoint the source of ANE-related ancestry in the Khövsgöl gene pool: that is, whether it derives from a pre-Bronze Age ANE population (such as the one represented by AG3) or from a Bronze Age WSH population that has both ANE and IRC ancestry. To test these competing hypotheses, we systematically compared various admixture models of the main cluster using the qpAdm program (20). Ancient Baikal populations were chosen as a proxy based on both their spatiotemporal and genetic similarities with the Khövsgöl main cluster (Figs. 2 and 3). When the early Neolithic Shamanka_EN is used as a proxy, we find that Baikal+ANE provides a better fit to the main cluster than Baikal+WSH, although no two-way admixture model provides a sufficient fit (P ≥ 0.05) (SI Appendix, Table S5). Adding a WSH population as the third source results in a sufficient three-way mixture model of Baikal+ANE+WSH with a small WSH contribution to the main cluster (e.g., P = 0.180 for Shamanka_EN+AG3+Sintashta with 3.7 ± 2.0% contribution from Sintashta) (Fig. 4 and SI Appendix, Table S6).

Using the temporally intermediate EBA population Shamanka_EBA, we can narrow down the time for the introduction of WSH ancestry into the main cluster. Shamanka_EBA is modeled well as a two-way mixture of Shamanka_EN and ANE (P = 0.158 for Shamanka_EN+AG3) (Fig. 4) but not as a mixture of Shamanka_EN and WSH (P ≤ 2.91 × 10⁻⁴) (SI Appendix, Table S5), suggesting no detectable WSH contribution through the early Bronze Age. Similar results are obtained for other Late Neolithic and EBA populations from the Baikal region (SI Appendix, Table S5). In contrast, the Khövsgöl main cluster is modeled well by Shamanka_EBA+WSH but not by Shamanka_EBA+ANE (P ≥ 0.073 and P ≤ 0.038, respectively) (SI Appendix, Table S5). A three-way model of Shamanka_EBA+ANE+WSH confirms this by providing the ANE contribution around zero (SI Appendix, Table S6). The amount of WSH contribution remains small (e.g., 6.4 ± 1.0% from Sintashta) (Fig. 4 and SI Appendix, Table S5). Assuming that the early Neolithic populations of the Khövsgöl region resembled those of the nearby Baikal region, we conclude that the Khövsgöl main cluster obtained ∼11% of their ancestry from an ANE source during the Neolithic period and a much smaller contribution of WSH ancestry (4–7%) beginning in the early Bronze Age.

Using the same approach, we obtained reasonable admixture models for the two outliers, ARS017 and ARS026. The eastern outlier ARS017, a female, shows an extra affinity with early Neolithic individuals from the Russian far east (Devil’s Gate) (22) and in general with contemporary East Asians (e.g., Han Chinese) compared with the Khövsgöl main cluster (Fig. 3B and SI Appendix, Fig. S12). ARS017 is also similar to Shamanka_EN in showing no significant difference in qpAdm (SI Appendix, Fig. S12 and Table S7). Using contemporary East Asian proxies, ARS017 is modeled as a mixture of predominantly Ulchi and a minor component (6.1–9.4%) that fits most ancient western Eurasian groups (P = 0.064–0.863) (SI Appendix, Table S7). This minor Western component may result from ANE ancestry; however, given the minimal western Eurasian contribution, we do not have sufficient power to accurately characterize this individual’s western Eurasian ancestry.

The Western outlier ARS026, a male dating to the end of the radiocarbon series, has the highest outgroup-f₃ with the main LBA Khövsgöl cluster, with extra affinity toward Middle Bronze Age (MBA) individuals from the Sintashta culture (Fig. 3B and SI Appendix, Fig. S13) (1). DNA recovered from this individual exhibited expected aDNA damage patterns (SI Appendix, Fig. S5) but was otherwise excellently preserved with >47% endogenous content and very low estimated contamination (1% mitochondrial; 0.01% nuclear). ARS026 is well modeled as a two-way mixture of Shamanka_EBA and Sintashta (P = 0.307; 48.6 ± 2.0% from Sintashta) (SI Appendix, Table S7). Similar to ARS026, contemporaneous LBA Karasuk individuals from the Altai (1400–900 BCE) (1, 29) also exhibit a strong extra genetic affinity with individuals associated with the earlier Sintashta and Andronovo cultures (SI Appendix, Fig. S14). Although two-way admixture models do not fit (P ≤ 0.045) (SI Appendix, Table S8), the Karasuk can be modeled as a three-way mixture of Shamanka_EBA/Khövsgöl and AG3 and Sintashta, suggesting an eastern Eurasian source with slightly higher ANE ancestry than those used in our modeling (P ≥ 0.186) (SI Appendix, Table S8). Like ARS026, admixture coefficients for the Karasuk suggest that MBA/LBA groups like the Sintashta or Srubnaya are a more likely source of their WSH ancestry than the EBA groups, like the Yamnaya or Afanasievo. Notably, Karasuk individuals are extremely heterogeneous in their genetic composition, with the genetically easternmost Eurasian individual nearly overlapping with the EBA Baikal groups (Fig. 2A and SI Appendix, Figs. S7 and S8). Earlier groups, such as the Afanasievo, Sintashta, and Andronovo, are mostly derived from WSH ancestries, and this may suggest that admixture in the Altai-Sayan region only began during the LBA following a long separation since the Eneolithic. Although ARS026 exhibits substantial WSH ancestry, strontium isotopic values obtained from his M3 enamel resemble local fauna and fall within the range of the main Khövsgöl cluster (SI Appendix, Fig. S15 and Table S9); however, because the enamel this individual also exhibited elevated manganese levels, postmortem trace element alteration from soil could not be excluded.

Contemporary Mongolia has a dairy- and meat-based subsistence economy, and to more precisely understand the role of dairy products in the diets of present-day mobile pastoralists in Khövsgöl aimag, we conducted a detailed nutritional investigation of summer and winter diets. We find that dairy-based foods contribute a mean of 35% total dietary energy, 36–40% total carbohydrate, 24–31% total protein, and 39–40% total fat to rural summer diets in Khövsgöl aimag, with liquid milk and dairy product consumption of 216–283 and 172–198 g/d, respectively (SI Appendix, Table S10 and Dataset S2).

Despite the importance of dairying today, its origins in Mongolia are poorly understood. Given the limited WSH ancestry of the main Khövsgöl cluster, we sought to determine if dairy pastoralism was practiced by this putatively pastoralist LBA population by testing for the presence of milk proteins (30) in the dental calculus of these individuals. We extracted proteins from 12 dental calculus samples representing 9 individuals (SI Appendix, Table S11) and analyzed tryptic peptides using LC-MS/MS (31). Observed modifications included deamidation (N, Q) and oxidation (P, M) (SI Appendix, Table S12). All protein identifications were supported by a minimum of two peptides across the dataset, and only peptides with an E value ≤ 0.001 were assigned; the estimated peptide false-discovery rate (FDR) across the full dataset was 1.0%, and protein FDR was 4.6%. Milk proteins were detected in seven of the nine individuals analyzed (SI Appendix, Table S13 and Dataset S3), confirming that dairy foods were consumed as early as 1456 BCE (1606–1298 BCE, 95% probability of the earliest directly dated individual) (SI Appendix, Fig. S4 and Table S2). Specifically, we detected the milk whey protein β-lactoglobulin (Fig. 5 A and B) and the curd protein α-S1-casein, with peptides matching specifically to sheep (Ovis), goat (Capra), Caprinae, Bovinae, and a subset of Bovidae (Ovis or Bovinae) (Fig. 5C, SI Appendix, Table S13, and Dataset S3). These peptides exhibited asparagine and glutamine deamidation, as expected for ancient proteins (32), and the frequency and distribution of recovered β-lactoglobulin (Fig. 5B) and α-S1-casein peptides closely matched that empirically observed for modern bovine milk (33), thereby providing additional protein identification support through appropriate proteotypic behavior.

Given the evidence for dairy consumption by the LBA Khövsgöl population, we sought to determine if the dairy-adaptive -13910*T (rs4988235) lactase persistence (LP) allele found today in Western steppe (34) and European (35) populations was present among LBA Khövsgöls dairy herders, and we examined this position in our SNP-enriched dataset. The -13910*T LP allele was not found in the LBA Khövsgöls (SI Appendix, Fig. S17 and Table S14), and additionally all observed flanking sequences in the lactase transcriptional enhancer region contained only ancestral alleles.

Based on an 850-km² archaeological survey of DSKC-associated burial mounds in Arbulag soum, Khövsgöl, Mongolia, we selected 22 burial mounds from 6 distinct burial mound groupings (A–F) for excavation and analysis (Fig. 1 and SI Appendix, sections 1 and 2 and Table S1). Bone and tooth samples from 22 individuals (11 femora, 11 teeth) were analyzed for ancient DNA, and 12 dental calculus samples from 9 individuals were analyzed for ancient proteins (SI Appendix, Table S2). Twenty-one individuals were successfully direct radiocarbon dated to ca. 1380–975 BCE (SI Appendix, section 3 and Table S2).

DNA extraction and library construction was performed in a dedicated clean room facility at the Max Planck Institute for the Science of Human History in Jena, Germany following previously published protocols (38), including partial uracil-DNA-glycosylase treatment (19). Following screening, 20 samples with ≥0.1% endogenous content were enriched for 1.2 million informative nuclear SNPs (1240K) by in-solution hybridization (2, 3). Additionally, preenrichment libraries for two well-preserved samples (ARS008 and ARS026) were deep-sequenced to generate ∼3.3× genomes. All sequencing was performed using single-end 75-bp (for screening and enriched libraries) or paired-end 50-bp (for whole-genome sequencing of two preenrichment libraries) sequencing on an Illumina HiSEq 4000 platform following the manufacturer’s protocols (SI Appendix, section 4).

DNA sequences were processed using the EAGER v1.92.50 pipeline (39). Adapter-trimmed reads ≥30 bp were aligned to the human reference genome using BWA aln/samse v0.7.12 (40) with the nondefault parameter “-n 0.01,” and PCR duplicates were removed using dedup v0.12.2 (39). The first and last three bases of each read were masked using the trimbam function in bamUtils v1.0.13 (41). For each target SNP, a single high-quality base (Phred-scaled quality score ≥30) from a high-quality read (Phred-scaled mapping quality score ≥30) was randomly chosen from the 3-bp masked BAM file to produce a pseudodiploid genotype for downstream population genetic analysis. DNA damage was assessed using mapDamage v2.0.6 (42), and mitochondrial DNA contamination was estimated using Schmutzi (43). For males, nuclear contamination was estimated using ANGSD v0.910 (44) (SI Appendix, section 4).

Mitochondrial haplogroups were determined by generating a consensus sequence using the log2fasta program in Schmutzi (43), followed by haplogroup assignment both by HaploGrep2 (45) and HaploFind (46). The Y haplogroup was determined using the yHaplo program (47). Genetic relatedness was estimated by calculating pairwise mismatch rate of pseudodiploid genotypes (48) (SI Appendix, section 4).

Khövsgöl SNP data were merged with published ancient genome-wide data for the 1240K panel (1, 3, 4, 20–22, 25–28, 49–59) (Dataset S1). A comparative dataset of present-day individuals was compiled from published datasets either genotyped on the Affymetrix Axiom Human Origins 1 array (HumanOrigins) or sequenced to high-coverage in the Simons Genome Diversity Project (20, 60–62) (SI Appendix, section 4). Intersecting with SNPs present in the HumanOrigins array, we obtain data for 593,124 autosomal SNPs across world-wide populations. Population structure was investigated by PCA as implemented in the smartpca v13050 in the Eigensoft v6.0.1 package (63) and by unsupervised genetic clustering using ADMIXTURE v1.3.0 (23) (SI Appendix, sections 4 and 5). The f₃ and f₄ statistics were calculated using the qp3Pop (v400) and qpDstat (v711) programs in the admixtools v3.0 package (24). For calculating the f₄ statistic, we added the “f4mode: YES” option to the parameter file. For admixture modeling, we used qpAdm v632 (20) in the admixtools v3.0 package (SI Appendix, sections 4 and 5).

Strontium isotopes (⁸⁷Sr/⁸⁶Sr) measured from human and faunal tooth enamel (n = 16) and bone (n = 5) were analyzed at the University of Georgia Center for Applied Isotope Studies (n = 17) and the University of Florida Department of Geological Sciences (n = 4) using a thermo-ionization mass spectrometer (SI Appendix, section 6).

Up to 6 d of weighed diet records were collected from 40 subjects (n = 231 total person-days) randomly sampled from the rural soum of Khatgal and the provincial center of Mörön in June 2012 and January 2013 by trained medical students from the Mongolian National University of Medical Sciences and Ach Medical Institute. Nutrient consumption was determined using a purpose-built food composition table (64), which we appended with unpublished food composition data from the Mongolian University of Science and Technology and the Mongolian Public Health Institute, as well as published data from the United States and Germany (65, 66) (SI Appendix, section 7). Contemporary dietary data were collected under Harvard Institutional Review Board Protocol #21002.

Ancient protein analysis was performed in a dedicated clean room facility at the Max Planck Institute for the Science of Human History following recommended guidelines (32). Dental calculus was decalcified in 0.5 M EDTA, and proteins were extracted and trypsin-digested using a modified low-volume Filter-Aided Sample Preparation protocol (67). The resulting peptides were analyzed by LC-MS/MS using a Q-Exactive HF mass spectrometer (Thermo Scientific) coupled to an ACQUITY UPLC M-Class system (Waters) at the Functional Genomics Center Zurich, according to previously published specifications (25). Extraction blanks and injection blanks were processed and analyzed alongside experimental samples (SI Appendix, section 8).

Raw spectra were converted to Mascot generic files using MSConvert using the 100 most intense peaks from each spectrum, and MS/MS ion database searching was performed using Mascot software (v2.6; Matrix Science) with the databases SwissProt (version 2017_07; 555,100 sequences) and a custom dairy database consisting of 244 dairy livestock milk protein sequences obtained from the National Center for Biotechnology Information GenBank. Before analysis, an error-tolerant search was performed (SI Appendix, Table S12) to identify common variable modifications (deamidation N, Q; oxidation M, P). Reversed sequences for each entry in both databases were added to perform downstream FDR calculations in R. Peptide tolerance was set at 10 ppm, with an MS/MS ion tolerance of 0.01 Da, and the data were filtered to only include peptides with an E-value ≤ 0.001 and proteins supported by a minimum of two peptides (SI Appendix, section 8). Peptides identified as matching milk proteins were tested for taxonomic specificity using BLASTp against the National Center for Biotechnology Information nr database and aligned to protein sequences of known dairy livestock. Modeling of β-lactoglobulin coverage was rendered using VMD v.1.9.4a7, and an additional level of protein identification confirmation was performed by comparing the observed ancient milk peptides to modern proteotypic peptides using the R package ggplot2 (68) with published data for bovine β-lactoglobulin obtained from the Peptide Atlas (33) (SI Appendix, section 8).

Genotype likelihoods for phenotype-associated SNPs were calculated using the UnifiedGenotyper program in the Genome Analysis Toolkit v3.5 (69) (SI Appendix, section 9).