High mitochondrial diversity of domesticated goats persisted among Bronze and Iron Age pastoralists in the Inner Asian Mountain Corridor

of in the precise pathways by which across Asia during the early to 2500 and later unclear. We analyzed sequences of hypervariable region 1 and cytochrome b gene in the mitochondrial genome (mtDNA) of goats from archaeological sites along two proposed transmission pathways as well as geographically intermediary sites. Unexpectedly high genetic diversity was present in the Inner Asian Mountain Corridor (IAMC), indicated by mtDNA haplotypes representing common A lineages and rarer C and D lineages. High mtDNA diversity was also present in central Kazakhstan, while only mtDNA haplotypes of lineage A were observed from sites in the Northern Eurasian Steppe (NES). These findings suggest that herding communities living in montane ecosystems were drawing from genetically diverse goat populations, likely sourced from communities in the Iranian Plateau, that were sustained by repeated interaction and exchange. Notably, the mitochondrial genetic diversity associated with goats of the

A second route for the spread of pastoralist subsistence to Central Asia reflects a vector of cultural transmission from the Iranian Plateau via a near continuous series of foothill ecozones to the Altai Mountains, dubbed the "Inner Asian Mountain Corridor" (IAMC) (Fig 1) [20]. Accordingly, pastoralist subsistence (perhaps through population movement or cultural transmission) is thought to have spread into Central Asia during the fourth millennium BC along montane elevational gradients as a strategy to provide graze for livestock by exploiting seasonal rhythms of vegetation growth [45]. Recent archaeological research in the IAMC has documented pastoralist occupations dating to the early third millennium BC that coincided with the initial trans-Eurasian transmissions of domesticated cereals of wheat and barley from the Near East and broomcorn millet from northern China [46][47][48]. Human paleogenetic research suggests that the IAMC in the second millennium BC was a distinct conduit for the northward human gene flow from agro-pastoralist communities from the region of Turan, roughly spanning present-day Turkmenistan and Uzbekistan [31]. In subsequent millennia, pastoralist interactions throughout the IAMC are associated with vibrant cultural exchanges, linking economies and cultures across Eurasia, which scholars commonly refer to as the ancient Silk Roads [45, [49][50][51][52][53][54].
While the NES and IAMC served as major cultural interaction spheres over the past six thousand years, there has been very little systematic evaluation of the dispersals of pastoralist subsistence or animal exchanges within and between these geographic loci. Here, we examine the mitochondrial (mtDNA) diversity of domesticated goats from sites located in the NES and IAMC with Bronze Age and/or Iron Age occupations in order to reconstruct the spread of this key pastoralist animal to Central Asia and subsequent genetic patterns (Fig 1) [55]. Key Bronze Age sites located between the NES and IAMC provide a test for the geographic extent of  [14], and supposed dispersal routes of domesticated goats are shown with dashed arrows. Made with Natural Earth. Free vector and raster map data @ http://www.naturalearthdata.com.
https://doi.org/10.1371/journal.pone.0233333.g001 maternally inherited goat genetic variation potentially characteristic of these macro regions (Fig 1). We analyzed the hypervariable region 1 (HVR1) of the control region and cytochrome b gene (MT-CYB) recovered from ancient goat skeletal remains. HVR1 sequences contain high frequencies of polymorphisms that characterize particular haplogroups, while MT-CYB sequences are highly conserved and generally reflect inter-species variation, which we use here to differentiate between domesticated sheep and goat for ambiguously identified archaeological specimens, in addition to wild Capra species [56], which are likely present in the zooarchaeological assemblages of the sampled sites. Prior to a major study across the Near East [14], a very limited number of ancient mtDNA sequences have been analyzed, which only represent Neolithic Europe [57] and the Near East [58][59][60] and also Iron Age China [61,62].

Sample selection, preparation, and DNA extraction
All necessary permits were obtained for the described study, which complied with all relevant regulations. Bones and teeth from 198 probable or possible goats, identified by morphological features [63,64] were selected from 11 archaeological sites with Bronze and/or Iron Age occupations in Kazakhstan, Kyrgyzstan, Russia, and Mongolia (Fig 1; Table 1, S1 Table). Wear stages of mandibular premolars and molars were scored following Payne [65], in order to differentiate individuals in case of recovering identical haplotypes from multiple specimens at each site. Laboratory work was conducted in the Ancient DNA Laboratory of Kiel University following established protocols [66][67][68][69]. Negative controls were employed during each set of DNA extractions and PCR preparations. Samples were soaked in consumer bleach solution for 5 min in order to remove external contamination, followed by a soaking for 5 min and two rinses with ddH 2 O. Samples were ground to a fine powder with a ball mill and digested in 960 μL 0.45 M EDTA (pH 8.0) and 40 μL proteinase K (25 mg/mL) overnight at 37˚C. DNA was extracted with a large-volume extraction protocol using a guanidinium-silica based method [70]. Two independent DNA extractions from sample powder were performed for each specimen yielding an MT-CYB sequence of Capra hircus, and sample amplicons were each sequenced twice, in forward and reverse directions.
Amplifications were performed as singleplex PCRs in a final volume of 25 μL (S2 Table). Reactions were performed in a Mastercycler (Eppendorf) with the following conditions: initialization at 94˚C for 10 min, 42 cycles of 94˚C for 30 s, 50-60˚C for 30 s (primer dependent, see S3 Table), 72˚C for 30 s, and a final elongation step at 72˚C for 10 min. The quality of PCR products was checked on a QIAxcel system, and those with single bands of expected size were subjected to direct sequencing with the BigDye Terminator v3.1 Kit (Applied Biosystems), following the instructions of the manufacturer. Sequencing products were purified applying the DyeEx 2.0 Spin Kit (Qiagen) and analyzed using the ABI Prism 310 and 3130 Genetic Analyzers (Applied Biosystems).
Ancient Capra hircus sequences of the control region and cytochrome b gene were concatenated for each sample, resulting in a total sequence length of 303-bp. Sequences were pair-wise aligned to one another in addition to nine previously published ancient sequences and 10 modern sequences of Capra spp. representing haplogroups A, B, C, D, F, and G (S4 Table). In cases where reference sequences of the control region did not accompany sequences of cytochrome b, we used cytochrome b sequences from sample KG170121. This approach is justified because of very low genetic variation in cytochrome b sequences for each Capra sp., which mainly served to improve phylogenetic clustering of wild Siberian ibex relative to domesticated goats. Sequences were trimmed to exclude the primer binding sites. A phylogenetic tree was constructed using MrBayes v. 3.2.6 [72] in Geneious v. 10.0.7. The simulation was set using the JC69 substitution model, invgamma rate variation, four gamma categories, and a sequence from Ovis aries as an outgroup; MCMC chain length was set to 500,000 with a sampling frequency of 200 and burn-in of 100,000. Additionally, a minimum-spanning haplotype network of sample and reference sequences was generated using PopART v. 1.7 [73]. Nucleotide diversity was calculated for each analytical region from the number of base substitutions per nucleotide position from averaging over all sequence pairs using the Kimura 2-parameter model and gamma distribution shape parameter of 1 [74], with standard error estimates by a bootstrap procedure with 1000 replicates in MEGA X v.10.1.1 [75]. Only sequences represented by 303-bp without gaps were included in the genetic distance analysis.

MT-CYB sequence recovery and taxonomic identification
Out of 198 ancient skeletal specimens identified as probable goats or ambiguously identified as sheep/goat (caprine), we recovered 50 MT-CYB sequences of Capra hircus, 94 MT-CYB sequences of Ovis aries, and two MT-CYB sequences of Capra sibirica ( Table 2). The relatively high overall recovery rate of 73.7% of MT-CYB sequences for Capra and Ovis spp. indicates good preservation of ancient mtDNA in the samples, which likely reflects the generally cool continental climate of sampled sites and high-elevation for sites in the IAMC, although there are appreciable differences in the recovery of MT-CYB sequences between sites ( Table 2). The Table 2. Recovery rates of MT-CYB sequences by site per taxon. MT-CYB sequences from Begash, Dali, and Tasbas were previously reported [47].

Frequencies Northern Eurasian Steppe
Chronology Samples tested Two specimens from Begash previously identified as Siberian ibex based on morphological criteria [86] were determined to be domesticated sheep, based on MT-CYB sequences that are specific to Ovis aries [47]. Two additional Siberian ibex specimens did not yield amplicons, likely due to DNA damage resulting from burning that was visible on bone surfaces. One Siberian ibex specimen from Dali dating to the Middle Bronze Age (ca. 1800-1500 cal BC), which was identified based on more reliable cranial morphology, yielded a Capra sibirica MT-CYB sequence, while one specimen previously identified as a domesticated caprine from Begash [86] yielded a Capra sibirica MT-CYB sequence [47].

MT-CYB
Misidentifications of Capra sibirica due to shared morphological characteristics with sheep and use of metrical determinations establishing large body-size, a trait typically associated with wild goats, highlights the need to exercise caution when attempting to separate wild and domesticated caprines. Under-identification of Siberian ibex in faunal assemblages limits our understanding of the importance of hunting by ancient pastoralists in Central Asia, which may have been a less common practice than previously suggested [86], and also obscures the significance of large-bodied, male rams in pastoralist livestock herds. In general, morphological criteria used to identify Central Asian sheep and goat draw from those developed for European and Near Eastern populations. Although the skeletal morphology of Siberian ibex is poorly studied, morphological similarity between domesticated sheep and wild ibex in some skeletal elements have been previously identified for both Near Eastern and European populations [95]. While collagen peptide mass fingerprinting (Zooarchaeology by Mass Spectrometry, aka "ZooMS") may be used to differentiate sheep and goat specimens as a relatively cheap biomolecular method [49,96,97], this technique is limited to genus-level resolution and, thus, cannot differentiate between wild and domesticated lineages.

Capra hircus haplogroup diversity
Out of 50 samples that yielded Capra MT-CYB sequences, we recovered 43 control region sequences. One specimen failed to generate the 193-bp amplicon but yielded a 130-bp sequence of HVR1 using primers CAP-FII and CAP-RII. Although the primers used to amplify sequences of the control region are generally Capra-specific, one sample from Taldysai produced a sequence for Ovis aries. Control region sequences were confirmed with duplicates taken from independent DNA extractions of each specimen.
Taken together, the ancient sequences of Capra hircus represent divergent mtDNA lineages, corresponding to haplogroups A, C, and D (Table 3), which are well characterized based on modern datasets consisting of more than 3000 domesticated and wild goats [15,17,18]. In our sample of Capra hircus mtDNA sequences, the Late Bronze Age (ca. 1700-1000 BC) is best represented (n = 26), followed by the early-middle Bronze Age (ca. 2700-1700 BC; n = 11) and the Iron Age (ca. 400 BC-AD 200; n = 5).
Diversity of mtDNA haplotypes was lowest for sites located in the NES, indicated by all recovered HVR1 sequences belonging to haplogroup A (Table 3) and an average genetic distance of 0.015±0.005. Low goat mtDNA diversity is in the NES is driven by the high occurrence of one HVR1 haplotype that was shared among three individual goats from Bozshakol and two from Air-Tau. This haplotype is the most common one in the dataset and was also present in five goats from sites located in central Kazakhstan (Fig 2). We are confident that the identical haplotypes from each site represent different individuals based on dental specimens being of the same body side or, in cases of specimens representing opposing body sides, individuality is confirmed by age distinction according to tooth wear stages (S1 Table). Interestingly, the remaining sample set from Taldysai is characterized by diverse HVR1 haplotypes, belonging to haplogroups A (n = 4), C (n = 2), and D (n = 1) ( Table 3; Fig 2), which together exhibit an average genetic distance of 0.042±0.008.
Mitochondrial diversity was substantially higher for sites located in the IAMC than the NES, shown by sequences belonging to haplogroups A (n = 10), C (n = 3), and D (n = 4) ( Table 3; Fig 2), which were unique, except for two C haplotypes from Begash and Uch-Kurbu (Fig 2). Notably, the IAMC is characterized by a high average genetic distance of 0.053±0.010, compared to the NES and central Kazakhstan.
Sequences of domesticated goats recovered from Dali, representing the earliest pastoralist occupations in the study [47], belong to haplogroup A (n = 4), while middle Bronze Age to Iron Age layers of IAMC sites contained an unexpectedly high 3:2 ratio of C and D haplotypes to A haplotypes (Table 3). This haplogroup distribution could indicate an influx of goats with diverse maternal origins into the IAMC after the middle Bronze Age, which strongly contrasts against the predominance of A lineages in post-Neolithic goats across the Near East. Although the sample size for the Iron Age is small, two out of three sequences are of lineage D, suggesting mtDNA diversity may have increased up to the end of the first century BC. This corresponds with the identification of one D lineage goat from Iron Age site of Bancheng in the Inner Mongolia Province of China [61] (Fig 1). Moreover, a recent study documenting mtDNA diversity in modern goats in Kazakhstan found that D lineages only occurred in low numbers in the southeastern region of the country compared to other regions that were dominated by A lineages [98]. Taken together, the diversity in mtDNA lineages from these studies and that in our sample reinforce the IAMC as an important locus connecting ancient Central Asian herders to those in east Asia via western China. Sequences representing haplogroups B, G, or F were not recovered in our study. We did not expect to identify F haplotypes, since goats belonging to the F lineage are the rarest in modern populations, having only been isolated in three Sicilian domesticated goats [99], in addition to wild bezoars from a wide geographic extent spanning the Caucasus and South Asia [15,16]. Managed goats belonging to haplogroup F (n = 5) were only present in the southern Levant during the Pre-Pottery Neolithic B, while ancient wild bezoars from southern Turkey and Armenia (n = 2) were also identified as belonging to lineage F [14]. We expected to recover B haplotypes from IAMC sites, since they are the second most abundant in modern distributions (ca. 6.5%) [15,17], and one B haplotype was found among Bronze Age goats in southern Uzbekistan [14] and two B haplotypes were found among Iron Age goats from northern China [61]. Goats of lineage B must have passed through Central Asia as pastoralists extended the social networks through which animals and crops were initially spread across Asia. The absence of G haplotypes in our sample was difficult to predict, since goats belonging to this lineage were recovered from Neolithic sites in Iran and Turkmenistan [14], which might have been transmitted into Central Asia at similar rates as goats belonging to lineage D. Today, goats of the G lineage are largely restricted to northern and central Africa and have a very low global frequency of 0.9% [15]. The absence of B and G haplotypes in our sample demands further investigation of ancient goat genetic diversity in Iran, the IAMC, and China in order to resolve processes leading to the persistence and extinction of goat lineages in particular times and places.

Network and phylogenetic analysis
Ancient goat sequences recovered in this study are scattered throughout the A clade, which is likely an outcome of A haplotypes being the most dominant in domesticated goats after large population expansions of A groups since domestication [14,16,17,19]. Haplotypes of the C lineage from Begash, Uch-Kurbu, and Taldysai cluster together and are distinct from modern reference haplotypes (Fig 2). Unexpectedly, two out of the four D haplotypes from goats in the IAMC during the Bronze and Iron Ages are newly discovered D haplotypes (Fig 2; Beg_LBA_3 and Beg_IA_1), which cluster phylogenetically as a distinct sub-clade within this haplogroup (Fig 3). This finding suggests the possibility of genetic continuity of rare domesticated goat haplotypes at Begash over a long period of time from 1600 to 20 cal BC, which spans prominent cultural turnovers during the Bronze to Iron Ages. Interestingly, no modern descendant or closely related sequences of this D lineage are known to exist. Haplotypes of lineage D are rare in contemporary world-wide distributions (0.6%) [15,17], but are more common in bezoar populations on the Iranian plateau (ca. 13%) [15]. Haplotypes of the D lineage have also been recovered from Neolithic and Bronze Age sites in Israel, Iran, Turkmenistan, and Uzbekistan [14], which cluster with the other two Bronze and Iron Age haplotypes from the IAMC on the main D clade (Fig 2). Given our recovery of a relatively high frequency of C and D haplotypes in the IAMC, ancient frequencies of these lineages were likely far higher in ancient Central Asia than in modern distributions. The presence of these lineages further indicates high mtDNA diversity of goats in the mountain regions of Central Asia compared to the NES that is suggestive of different intensities and geographic extents of social interaction between and among these regions.
Due to random genetic drift ultimately leading to monomorphism, a single lineage may become dominant in a diverse population after many generations in a population that has little inward gene flow of assorted lineages [100,101]. While it is impossible to estimate the precise population size of any particular livestock species based on zooarchaeological data, goat skeletal remains are about 10-15 times less frequent than those of sheep and about five times as infrequent as cattle remains from sites in the IAMC during the Bronze and Iron Ages [46,86,93], indicating that goat population sizes were relatively small compared to other domesticated species. Without a large population to buffer the effects of genetic drift, the persistence of C and D lineages in goats is likely due to high-rates of gene flow along the IAMC, likely connecting goat herds in the IAMC with diverse goat populations in southern Central Asia or further southwest on the Iranian Plateau at consistent rates. Frequent exchange events over relatively short distances could have linked goat herds of pastoralist communities as a large, global effective population, thus maintaining the prevalence of less common mtDNA lineages compared to the NES. A high degree of community interaction, likely driven by both seasonal mobility of pastoralists [20,45,102] and sprawling economic networks for the mining and exchange of copper and tin to produce bronzes during the Bronze Age and Iron Age [36, [103][104][105], could have been a main factor for promoting inter-regional interactions that went hand in hand with gene flows among livestock herds at inter-regional scales.
In contrast, the NES was characterized by goats exhibiting A haplotypes with low genetic diversity, as indicated by identical sequences among NES sites and a higher connectivity within the haplotype network (Fig 2). This pattern could have been the result of domesticated goats spreading eastward across the NES from a monomorphic founding population of goats, which later during the Bronze and Iron Ages was not characterized by inward gene flows of goats exhibiting C or D haplotypes that were present in the IAMC. While the NES has long been conceived as an important ecozone of pastoralist mobility and interaction, our observation of low goat mtDNA variation in this region suggests that, at least in part, there were limited communication and exchange with southern regions in Central Asia, especially the IAMC, where rare goat haplotypes and higher mtDNA diversity persisted in the Bronze and Iron Ages.
Interestingly, Taldysai located in central Kazakhstan was characterized by an intermediate level of goat mtDNA diversity relative to the NES and IAMC (Fig 2). This finding provides an important reconfiguration in how traditional cultural interaction spheres of the Eurasian steppe zone are defined, suggesting that central Kazakhstan served as an arena of interaction between both north and south regions but did not function as a corridor for goat gene flow connecting the north and south together. This finding corresponds to recent research on regional variance in domesticated sheep astragali morphology using geometric morphometrics, showing that localized populations likely were genetically isolated from one another [106]. However, the presence of C and D goat haplotypes at Taldysai suggests repeated exchange events with southern Central Asia and the IAMC, but these lineages likely did not pass along to NES sites. Human paleogenetic studies have described a substantial amount of human gene flow from pastoralist communities in the western NES to southern Central Asia during the second millennium BC [107,108], which would have most certainly been accompanied by goats predominantly characterized by A haplotypes. Similarly, the recent discovery by Narasimhan et al.
[31] of limited human gene flows from Turan and the IAMC into central Kazakhstan during this period would explain the presence of some goat lineages at Taldysai that were recovered in high frequencies at IAMC sites.
The one ancient Capra sibirica sequence recovered from Dali phylogenetically clustered as a sister group to previously published modern Capra sibirica sequences, which is supported by a 99% Mr. Bayes posterior probability (Fig 3). However, the recovered MT-CYB and HVR1 sequences from this Capra sibirica specimen did not cluster with the Capra sibirica reference sequences in the haplotype network (Fig 2), which is likely due to highly divergent mitochondrial variation between ancient and modern populations; the former probably had far greater genetic variation than the latter. To date, our Capra sibirica sequences represent the only ancient genetic sequences recovered from Capra sibirica and demonstrates a newly discovered degree of diversity in this species, while so far being poorly characterized genetically and biogeographically [109][110][111][112][113].

Conclusions
Despite small sample sizes, the results of this study suggest higher levels of gene flow among domesticated goats along the IAMC, which is likely a reflection of increased community interaction in the foothill zones, compared to the open steppe landscape. The northern "Northern Eurasian Steppe" is typically viewed as a metaphorical highway of cultural interaction between eastern Europe and Siberia, since the initial spread of pastoralism to the Altai Mountains via purported migrations of Afanasievo communities [26] and through the formation of expansive cultural horizons of the Iron Age [39][40][41][42][43][44]. Based on the low mtDNA diversity present in Bronze Age and later Iron Age domesticated goats in the NES, there may be an unexpectedly low degree of admixture between northern steppe pastoralist communities and those in adjacent regions. The goats analyzed from sites in the NES appear to be consistent with a founder population that was likely present in eastern Europe after an earlier arrival of domesticated goats characterized by A haplotypes from the Near East via Eastern Europe or the Caucasus.
In the IAMC, we observed a possible shift in the frequencies of mtDNA lineages from the early-middle Bronze Age (2700-1600 cal BC) to the late Bronze Age and Iron Age (1600-20 cal BC), which may be a result of reconfigured socio-economic networks that facilitated transfers of goats belonging to C and D lineages to the region from diverse genetic pools from agropastoralist communities located in the Iranian plateau. In the absence of individual animals moving long distances along the IAMC, genetic isolation between discrete goat populations may have been greatly diminished due to directed breeding patterns of goats by pastoralists that effectively created a macro-population of goats, thus maintaining diverse maternal lineages. These findings provide additional support that the IAMC was a dynamic arena of cultural activity spanning distant regions of Central and Inner Asia, which ultimately linked these regions to contemporaneous cultural arenas in east Asia [20]. From Taldysai, the recovery of C and D haplotypes and several A haplotypes that were identical to those from NES sites suggests the region of central Kazakhstan was discretely connected to interaction spheres of communities in the NES and IAMC but may have marked a boundary between them.
More research is needed to resolve the precise nature of livestock translocations, interregional genetic networks, and selective pressures on livestock. This study offers a muchneeded expansion of datasets for characterizing ancient goat mitochondrial diversity in Eurasia, and larger sample sizes spanning more sites along proposed transmission routes will improve the quality of inferences concerning inter-regional human interactions that can be drawn from maternal lineages. The addition of other domesticated taxa, such as sheep and cattle, to future paleogenetic analysis will help establish differences in genetic sourcing from various geographical regions. By also examining genetic variability in wild taxa, ancient biogeographic patterns may be identified to reconstruct a more complete picture of human engagements with the landscape that may reveal critical patterns in the relationship between hunting and herding strategies. The alignment of new evidence for these subsistence pursuits with recent findings of early and intensive cultivation of domesticated crops by Central Asian communities during the Bronze Age will further clarify the role that economic diversification played in shaping regional and trans-regional social interaction spheres.
Supporting information S1 Text. Descriptions of archaeological sites and references. (DOCX) S1 Table. List of aDNA samples, taxonomic identifications, haplogroup designations, skeletal element identifications, tooth wear stages, archaeological context information, and archival location. Mandibular wear stages follow Payne [65]. Maxillary wear stages are based on a relative system developed by zooarchaeologist Sarah Pleuger: scale 0-5, denoted with symbol � . (XLSX) Fig 1), in addition to previously published ancient and contemporary reference sequences from Eurasia; Ovis aries was used as the outgroup. Posterior probabilities � 60% are shown on each branch. Haplogroup designations of Capra hircus are listed to the right of the sample ID numbers.