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

Recent paleogenomic studies have shown that migrations of Western steppe herders (WSH) beginning in the Eneolithic (ca. 3300-2700 BCE) profoundly transformed the genes and cultures of Europe and central Asia. Compared with Europe, however, the eastern extent of this WSH expansion is not well defined. Here we present genomic and proteomic data from 22 directly dated Late Bronze Age burials putatively associated with early pastoralism in northern Mongolia (ca. 1380-975 BCE). Genome-wide analysis reveals that they are largely descended from a population represented by Early Bronze Age hunter-gatherers in the Baikal region, with only a limited contribution (∼7%) of WSH ancestry. At the same time, however, mass spectrometry analysis of dental calculus provides direct protein evidence of bovine, sheep, and goat milk consumption in seven of nine individuals. No individuals showed molecular evidence of lactase persistence, and only one individual exhibited evidence of >10% WSH ancestry, despite the presence of WSH populations in the nearby Altai-Sayan region for more than a millennium. Unlike the spread of Neolithic farming in Europe and the expansion of Bronze Age pastoralism on the Western steppe, our results indicate that ruminant dairy pastoralism was adopted on the Eastern steppe by local hunter-gatherers through a process of cultural transmission and minimal genetic exchange with outside groups.

Archaeogenetic studies provide evidence that the Eurasian Eneolithic–Bronze Age transition was associated with major genetic turnovers by migrations of peoples from the Pontic-Caspian steppe both in Europe and in central Asia (1–5). The migration of these Western steppe herders (WSH), with the Yamnaya horizon (ca. 3300–2700 BCE) as their earliest representative, contributed not only to the European Corded Ware culture (ca. 2500–2200 BCE) but also to steppe cultures located between the Caspian Sea and the Altai-Sayan mountain region, such as the Afanasievo (ca. 3300–2500 BCE) and later Sintashta (2100–1800 BCE) and Andronovo (1800–1300 BCE) cultures. Although burials typologically linked to the Afanasievo culture have been occasionally reported in Mongolia (6), the genetic profile of Eastern steppe populations, as well as the timing and nature of WSH population expansion and the rise of dairy pastoralism in Mongolia, remain unclear.

The remarkable demographic success of WSH populations has been linked to mobile pastoralism with dairying (7), a system that efficiently converts cellulose-rich wild grasses into protein- and fat-rich dairy products. Dairy foods provide a rich source of nutrients and fresh water, and function as an adaptive subsistence strategy in cold, dry steppe environments where most crop cultivation is highly challenging. Dairy pastoralism became widely practiced in the eastern Eurasian steppe, the homeland from which subsequent historical nomadic dairying empires, such as the Xiongnu (ca. 200 BCE to 100 CE) and the Mongols (ca. 1200–1400 CE) expanded; however, it is not fully understood when, where, and how this subsistence strategy developed. At Botai, in central Kazakhstan, evidence for Eneolithic dairying has been reported through the presence of ruminant and equine dairy lipids in ceramic residues as early as 3500 BCE (8, 9). In the Altai and Tarim basin, where WSH populations have left strong genetic footprints (1, 3, 10, 11), archaeological evidence supports the presence of dairy products by the early second Millennium BCE and later (8, 12, 13). In the Eastern steppe, however, no direct observations of dairy consumption have been made for a comparable time period, despite the fact that skeletal remains of domestic livestock (such as sheep, goats, cattle, and horses) have been found at Mongolian ritual sites and in midden contexts as early as the 14th century BCE (14–17). In the absence of direct evidence for Bronze Age milk production or consumption on the Eastern steppe, it remains unclear whether these animals are merely ritual in nature or signify a major shift in dietary ecology toward dairy pastoralism, and whether their appearance is connected to possible WSH migrations onto the Eastern steppe.

To understand the population history and context of dairy pastoralism in the eastern Eurasian steppe, we applied genomic and proteomic analyses to individuals buried in Late Bronze Age (LBA) burial mounds associated with the Deer Stone-Khirigsuur Complex (DSKC) in northern Mongolia (). To date, DSKC sites contain the clearest and most direct evidence for animal pastoralism in the Eastern steppe before ca. 1200 BCE (18). Focusing on six distinct burial clusters in Arbulag soum, Khövsgöl aimag, Mongolia (Fig. 1 and ), we produced genome-wide sequencing data targeting ∼1.2M single nucleotide polymorphisms (SNPs) for 22 DSKC-associated individuals directly dated to ca. 1380–975 calibrated BCE (), as well as sequenced whole genomes for two individuals (>3× coverage). Nine of the individuals in this group yielded sufficient dental calculus for proteomic analysis, and we tested these deposits for the presence of milk proteins using liquid chromatography-tandem mass spectrometry (LC-MS/MS). Overall, our results find that DSKC subsistence strategy included dairying of Western domesticated ruminants, but that there was minimal gene flow between analyzed DSKC populations and WSH groups during the LBA. Thus, in contrast to patterns observed in western Europe where, for example, the arrival of WSH is associated with population replacement and continental-level genetic turnover (5), contact between WSH and Eastern steppe populations is characterized by transcultural transmission of dairy pastoralism in the near absence of demic diffusion.

Fig. 1.

Map of the Eurasian steppes. (A) Distribution of the Western (brown) and Eastern (green) steppes and the locations of ancient (red) and modern (black) populations discussed in the text. Population codes are provided in the Dataset S1. A box indicates the location of the LBA burial mounds surveyed in the Arbulag soum of Khövsgöl aimag. (B) Enhanced view of LBA burial mounds (white circles) and burial clusters selected for excavation (boxes a–f) with the number of analyzed individuals in parentheses (). (C) Photograph of burial 2009-52 containing the remains of ARS026, a genetic outlier with Western steppe ancestry.

Results

Ancient DNA Sequencing and Quality Assessment.

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 (). 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 (), 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) (). No close relative pairs were identified among the ancient individuals (). Two individuals with high endogenous content on screening (ARS008, 70.0%; ARS026, 47.6%) were deeply sequenced to obtain whole genomes (∼3.3× coverage) (). 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).

Characterization of the Genetic Profile of the Khövsgöl Gene Pool.

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 (). PC1 separates eastern and western Eurasian populations, with central and north Eurasian populations falling in an intermediate position (). 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 () (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. 2 and ). 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.

Fig. 2.

The genetic profile of LBA Khövsgöl individuals summarized by PCA and ADMIXTURE. (A) Khövsgöl (Kvs, ARS017, and ARS026) and other ancient individuals (colored symbols) are projected onto the top PCs of modern Eurasian and Native American individuals. Contemporary individuals are marked by gray circles. Mean coordinates for each of the contemporary populations are marked by three-letter codes and by colors assigned to the associated geographic regions. Population codes are provided in Dataset S1 and . (B) ADMIXTURE results for Khövsgöl and other ancient individuals with K values 9 and 17. In K = 17, the Khövsgöls main cluster is mainly modeled as a mixture of components most enriched in modern northeast Asians (e.g., Nivh) and ancient Siberians (e.g., AG3, Botai, and Okunevo).

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. 2). Genetic clustering with ADMIXTURE (23) further supports these patterns (Fig. 2 and ). 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öl1, Khövsgöl2) for all pairs against 18 outgroups representative of world-wide ancestries (). 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. 2).

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. 3 and ). 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. 3 and ). 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 ().

Fig. 3.

The genetic affinity of the Khövsgöl clusters measured by outgroup-f and -f statistics. (A) The top 20 populations sharing the highest amount of genetic drift with the Khövsgöl main cluster measured by f(Mbuti; Khövsgöl, X). (B) The top 15 populations with the most extra affinity with each of the three Khövsgöl clusters in contrast to Tuvinian (for the main cluster) or to the main cluster (for the two outliers), measured by f(Mbuti, X; Tuvinian/Khövsgöl, Khövsgöl/ARS017/ARS026). Ancient and contemporary groups are marked by squares and circles, respectively. Darker shades represent a larger f statistic. Population codes are provided in Dataset S1; see also for further details.

Source of ANE Ancestry in the LBA Khövsgöl Population.

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) (). 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 ).

Fig. 4.

Admixture modeling of Altai populations and the Khövsgöl main cluster using qpAdm. For the archaeological populations, (A) Shamanka_EBA and (B and C) Khövsgöl, each colored block represents the proportion of ancestry derived from a corresponding ancestry source in the legend. Error bars show 1 SE. (A) Shamanka_EBA is modeled as a mixture of Shamanka_EN and AG3. The Khövsgöl main cluster is modeled as (B) a two-way admixture of Shamanka_EBA+Sintashta and (C) a three-way admixture Shamanka_EN+AG3+Sintashta. Details of the admixture models are provided in .

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−4) (), 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 (). 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) (). A three-way model of Shamanka_EBA+ANE+WSH confirms this by providing the ANE contribution around zero (). The amount of WSH contribution remains small (e.g., 6.4 ± 1.0% from Sintashta) (Fig. 4 and ). 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.

Admixture Testing of Genetic Outliers.

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. 3 and ). ARS017 is also similar to Shamanka_EN in showing no significant difference in qpAdm (). 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) (). 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. 3 and ) (1). DNA recovered from this individual exhibited expected aDNA damage patterns () 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) (). 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 (). Although two-way admixture models do not fit (P ≤ 0.045) (), 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) (). 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. 2 and ). 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 (); however, because the enamel this individual also exhibited elevated manganese levels, postmortem trace element alteration from soil could not be excluded.

Dairy Subsistence and Lactase Persistence.

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 ( 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 () and analyzed tryptic peptides using LC-MS/MS (31). Observed modifications included deamidation (N, Q) and oxidation (P, M) (). 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 ( 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) (). Specifically, we detected the milk whey protein β-lactoglobulin (Fig. 5 ) 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. 5, , 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. 5) and α-S1-casein peptides closely matched that empirically observed for modern bovine milk (33), thereby providing additional protein identification support through appropriate proteotypic behavior.

Fig. 5.

Presence of ruminant β-lactoglobulin and α-S1-casein milk protein in LBA Khövsgöl dental calculus. (A) B- and Y-ion series for one of the most frequently observed β-lactoglobulin peptides, TPEVD(D/N/K)EALEKFDK, which contains a genus-specific polymorphic residue: D, Bos; N, Ovis; K, Capra. See for peptide and fragment ion error distribution graphs. (B) Alignment of observed peptides to the 178 amino acid β-lactoglobulin protein, with peptide taxonomic source indicated by color. Trypsin cut sites are indicated by gray ticks. The position and empirically determined observation frequency of BLG peptides for bovine milk are shown as a heatmap scaled from least observed peptides (light gray) to most frequently observed peptides (dark red), as reported in the Bovine PeptideAtlas (34). Inset displays a 3D model of the β-lactoglobulin protein with observed peptide positions highlighted in black. (C) Taxonomically assigned β-lactoglobulin (black) and α-S1-casein (gray) peptides presented as scaled pie charts on a cladogram of Mongolian dairy domesticates. Bracketed numbers represent the number of peptides assigned to each node. Ruminant milk proteins were well supported, but no cervid, camelid, or equid milk proteins were identified.

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 (), and additionally all observed flanking sequences in the lactase transcriptional enhancer region contained only ancestral alleles.

Discussion

In this study, we find a clear genetic separation between WSH populations and LBA Mongolians more than a millennium after the arrival of WSH at the furthest edges of the Western steppe and the earliest appearance of the WSH Afanasievo cultural elements east of the Altai-Sayan mountain range. This genetic separation between Western and Eastern steppe populations appears to be maintained with very limited gene flow until the end of the LBA, when admixed populations, such as the Karasuk (1200–800 BCE), first appear in the Altai (1) and we observe the first individual with substantial WSH ancestry in the Khövsgöl population, ARS026, directly dated to 1130–900 BCE. Consistent with these observations, we find that the WSH ancestry introduced during these admixture events is more consistent with MBA and LBA steppe populations, such as the Sintashta (2100–1800 BCE), than with earlier EBA populations, such as the Afanasievo (3300–2500 BCE), who do not seem to have genetically contributed to subsequent populations.

Despite the limited gene flow between the Western and Eastern steppes, dairy pastoralism was nevertheless adopted by local non-WSH populations on the Eastern steppe and established as a subsistence strategy by 1300 BCE. Ruminant milk proteins were identified in the dental calculus of most of the tested LBA Khövsgöl individuals, and all identified milk proteins originated from ruminants, specifically the Western dairy domesticates sheep, goat, and Bovinae. These findings suggest that neighboring WSH populations directly or indirectly introduced dairy pastoralism to local indigenous populations through a process of cultural exchange. Further research on other regional cultures in Mongolia, such as Chemurchek, Hemsteg, and Ulaanzuukh, is needed to determine if this pattern of cultural adoption observed among DSKC sites is broadly shared across other Bronze Age cultures throughout the Eastern steppe.

Bronze Age trade and cultural exchange are difficult to observe on the Eastern steppe, where mobile lifestyles and ephemeral habitation sites combine to make household archaeology highly challenging. Burial mounds are typically the most conspicuous features on the landscape, and thus much of Mongolian archaeology is dominated by mortuary archaeology. However, unlike WSH, whose kurgans typically contain a range of grave goods, many LBA mortuary traditions on the Eastern steppe did not include grave goods of any kind other than ritually deposited animal bones from horse, deer, and bovids. Given that Mongolian archaeological collections are typically dominated by human remains with limited occupational materials, the ability to reconstruct technological exchange, human–animal interaction, and secondary product utilization through the analysis of proteins preserved in dental calculus represents an important advance.

The 3,000-y legacy of dairy pastoralism in Mongolia poses challenging questions to grand narratives of human adaptation and natural selection (36). For example, despite evidence of being under strong natural selection (36), LP was not detected among LBA Khövsgöls, and it remains rare (<5%) in contemporary Mongolia even though levels of fresh and fermented dairy product consumption are high (35). Recent studies in Europe and the Near East have found that dairying preceded LP in these regions by at least 5,000 y, suggesting that LP may be irrelevant to the origins and early history of dairying (36). As a non-LP dairying society with a rich prehistory, Mongolia can serve as a model for understanding how other adaptations, such as cultural practices or microbiome alterations (37), may be involved in enabling the adoption and long-term maintenance of a dairy-based subsistence economy. Early herding groups in Mongolia present a historical counter-example to Europe in which WSH migrations resulted in cultural exchange rather than population replacement, and dairying was maintained for millennia without the introgression or selection of LP alleles.

Materials and Methods

Experimental Design.

Based on an 850-km2 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 ). 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 (). Twenty-one individuals were successfully direct radiocarbon dated to ca. 1380–975 BCE ().

Ancient DNA Extraction, Library Construction, and Sequencing.

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 ().

DNA Sequence Data Filtering and Quality Assessment.

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) ().

Uniparental Haplogroup and Kinship Analysis.

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) ().

Population Genetic Analysis.

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) (). 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) (). 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 ().

Strontium Isotope Analysis.

Strontium isotopes (87Sr/86Sr) 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 ().

Dietary Analysis in Contemporary Khövsgöl, Mongolia.

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) (). Contemporary dietary data were collected under Harvard Institutional Review Board Protocol #21002.

Protein Extraction, Digestion, and LC-MS/MS.

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 ().

Spectrum Analysis, Data Filtering, and Authentication.

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 () 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 (). 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) ().

Phenotype-Associated SNPs.

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