Ancient genomes reveal tropical bovid species in the Tibetan Plateau contributed to the prevalence of hunting game until the late Neolithic.

Local wild bovids have been determined to be important prey on the northeastern Tibetan Plateau (NETP), where hunting game was a major subsistence strategy until the late Neolithic, when farming lifestyles dominated in the neighboring Loess Plateau. However, the species affiliation and population ecology of these prehistoric wild bovids in the prehistoric NETP remain unknown. Ancient DNA (aDNA) analysis is highly informative in decoding this puzzle. Here, we applied aDNA analysis to fragmented bovid and rhinoceros specimens dating ∼5,200 y B.P. from the Neolithic site of Shannashuzha located in the marginal area of the NETP. Utilizing both whole genomes and mitochondrial DNA, our results demonstrate that the range of the present-day tropical gaur (Bos gaurus) extended as far north as the margins of the NETP during the late Neolithic from ∼29°N to ∼34°N. Furthermore, comparative analysis with zooarchaeological and paleoclimatic evidence indicated that a high summer temperature in the late Neolithic might have facilitated the northward expansion of tropical animals (at least gaur and Sumatran-like rhinoceros) to the NETP. This enriched the diversity of wildlife, thus providing abundant hunting resources for humans and facilitating the exploration of the Tibetan Plateau as one of the last habitats for hunting game in East Asia.

Interactions among climate change, biogeography, and human activity have remarkably strengthened since the late Pleistocene. Abrupt climate change is attributed as a significant factor in worldwide human dispersal (1, 2) and the transformation of megafauna in the Arctic region in the late Pleistocene (3) and contributed to the agricultural revolution in the early Holocene (4). Farmers diffused across the Old World with crops and livestock during the Neolithic period (∼10,000 to 4,000 y B.P.) (5), eventually affecting the biogeography in Eurasia (6). However, the impact of the geographical distribution of wildlife on human behavior during that period remains enigmatic.

In East Asia, millet crops were domesticated ∼10,000 y B.P. (7), with a millet-based agricultural lifestyle established in the Loess Plateau ∼6,000 y B.P. (8). The millet farmers of the late Yangshao (5,500 to 5,000 y B.P.) and Majiayao (5,300 to 4,000 y B.P.) cultures expanded westward to the northeastern Tibetan Plateau (NETP) and settled extensively below 2,500 m above sea level (m a.s.l.) (9–11). Even after the introduction of agriculture to the NETP over 1,200 y B.P., hunting game remained a major subsistence strategy for the indigenous Zongri societies (∼5,300 to 4,000 y B.P.) (12) that inhabited the plateau between 2,800 and 3,000 m a.s.l. Therefore, the NETP served as one of the last habitats for the prevalence of hunting strategies in East Asia.

In contrast with previous zooarchaeological studies on early Majiayao groups, abundant wild faunal remains were unearthed from the early Majiayao site of Shannashuzha (SNSZ) (34°29′34.06″N, 104°4′24.54″E Min County, Gansu Province), which is dated to ∼5,280 to 5,050 y B.P. and located in the marginal area of the NETP at 2,323 m a.s.l. (13). The proportion of the identified wildlife remains (∼76%) far exceeded domestic animals (pig and dog) (∼24%) (14), while the plant remains were dominated by carbonized foxtail and broomcorn grains (over 80%) (13, 15). This suggests that both millet cultivation and game hunting occurred at SNSZ, despite previous zooarchaeological studies suggesting that early Majiayao groups mainly utilized livestock rather than wildlife for meat protein consumption (16). The wild animal remains at SNSZ include bones of deer (Cervidae), gorals (Nemorhaedus griseus), snub-nosed monkeys (Rhinopithecus roxellana), and rhinoceroses (Dicerorhinus). Bovid remains have also been identified at SNSZ and they were tentatively recognized as Tibetan yak or domestic cattle. However, they are too degraded for species identification using morphological approaches (14). Ancient DNA (aDNA) is a powerful tool for identifying bovid remains from archaeological sites in northern China, including domestic cattle (Bos taurus), wild water buffalo (Bubalus mephistopheles), and aurochs (Bos primigenius) (17–20). Although bovid remains have been frequently identified in many prehistoric sites on the Tibetan Plateau (7, 9, 12, 21), the species-level identification of bovid species on the Tibetan Plateau has not been reported.

Using aDNA profiles, we provide compelling evidence that ∼5,200 y B.P. the range of the present-day tropical gaurs (Bos gaurus) and Sumatran-like rhinoceroses (Dicerorhinus sumatrensis) extended as far north as the margins of the NETP at the northernmost latitude of 34.49°N. During this period the mean annual and summer temperatures in the Northern Hemisphere, especially in the middle-high latitude regions, were much higher than those during the preindustrial period (e.g., the 19th century) (22–24). These two tropical species currently inhabit low-latitude tropical or subtropical regions; however, to our knowledge (25–27), no other evidence has shown that the modern gaur’s range previously expanded into northern China. Therefore, our results provide evidence for revisiting the connection between climate change, biogeographic distributions, and the geographical patterns of human subsistence strategies in East Asia during the late Neolithic period.

Results

aDNA Authenticity, Preservation, and Radiocarbon Dating of Bovid Faunal Remains from the SNSZ Site.

The SNSZ archeological site is located in Min County, Gansu Province, the margin of the NETP. The site is adjacent to the Tao River, which is a tributary of the upper Yellow River, with surrounding high mountainous terrain to the northeast and southwest (Fig. 1). Between 1981 and 2010, the meteorological station records of Min County, which lies at an elevation of 2,315 m a.s.l. and a distance ∼26 km south of the site (www.nmic.cn), documented a mean annual temperature of ∼6.5 °C, a mean summer temperature of 15.5 °C, and a mean annual precipitation of ∼556 mm.

Fig. 1.

The distribution of gaur and Sumatran rhinos in Asia and the location of the SNSZ site. Differently colored regions represent the historic and modern distribution of the two species. The historical distribution of Sumatran rhinos was derived from range maps in Foose et al. (27), and the current distribution was derived from IUCN (2008) (26). The current distribution of gaur was derived from IUCN (25). SNSZ represents the Shannashuzha site, GH represents Gonghai Lake, XY represents Xingyun Lake, and QH represents Qinghai Lake.

According to the excavated artifacts and cultural relics, the community of the SNSZ site is regarded as the early phase of Majiayao culture, which was a late Neolithic culture that prevailed 5,300 to 4,000 y B.P. in the NETP around the valley of the upper Yellow River (16). Systematic archaeological surveys of the SNSZ reported typical mountain vegetation including conifer forests, mixed broadleaved/conifer forests, riverbank broadleaved forests, and, in particular, an abundant Bambusoideae botanical composition (14, 28). In addition, zooarchaeological classification identified a significantly higher proportion (>76%) of the number of identified specimens (NISP) of wild animal remains than of domestic animals (13) (), which include tropical inhabiting species such as bamboo rats (Rhizomys pruinosus), rhinoceroses (Dicerorhinus), and considerable unknown bovid remains (14) (). However, the species affiliation of these bovid bones is unknown due to the morphological overlap among bovid species’ faunal remains. We performed accelerator mass spectrometry (AMS) radiocarbon dating for five selected animal bones from the SNSZ site, including one bovid bone for species identification, and the 14C dating ranged from ∼5,280 to 5,050 y B.P. (Table 1).

Table 1.

Calibrated radiocarbon dates of animal bones from the excavation of the SNSZ site

Laboratory no.	Provenience	Dating method	Dating material	Radiocarbon age, y B.P.	Calibrated age, cal y B.P.
					1σ (68.2%)	2σ (95.4%)
LZU149	2013MCSIT0309⑤	AMS	Collagen from Cervidae	4,510 ± 25	5,176 ± 114	5,174 ± 124
LZU154	T0208⑤	AMS	Collagen from Bovidae	4,490 ± 25	5,165 ± 114	5,166 ± 123
Beta-539812	SNSZT0110H44	AMS	Collagen from Rhinoceros	4,450 ± 30	5,121 ± 147	5,089 ± 196
Beta-539813	SNSZTO309H113①*	AMS	Collagen from Bovidae	4,520 ± 30	5,181 ± 116	5,177 ± 127
Beta-539814	2013MCSIT0308④	AMS	Collagen from Bovidae	4,480 ± 30	5,162 ± 116	5,134 ± 155

Sample included in aDNA analysis.

We performed the aDNA laboratorial work in the aDNA dedicated laboratory of the Comparative Genomic Group (CGG) at Kunming Institute of Zoology, Chinese Academy of Sciences. Using the initial shotgun-sequencing screen data, we first assessed the aDNA authenticity and preservation situation by evaluating several key indexes, including the terminal damage rate (TDR) (), the average fragmental length, and the mapping rate (). The TDRs for all 10 samples were ≥10%, and the average fragmental lengths and the mapping rates range from 56 bp to 80 bp and 0.001 to 3.75%, respectively (), all of which indicate that all samples have molecular characteristics typical of aDNA. In addition, for the captured sequence data, we plotted the sequence fragmental size distribution for each sample, indicating a short fragment enrichment pattern, which serves as the typical pattern of aDNA ().

Molecular Identification of Species Affiliation Reveals the Present-Day Tropical Rainforest–Inhabiting Gaur and Sumatran Rhinoceros Once Roamed over the NETP ∼5,200 y B.P.

We collected 10 bone samples to perform molecular identification of these bovid faunal remains from the SNSZ site (). Identification of the bovid specimens was first achieved through mitochondrial DNA (mtDNA) analysis, with all 10 samples aligning the most reads to the modern-day gaur (NC_024818.1) mitogenome, as opposed to other bovid species such as the yak or domestic taurine cattle. By combining the shotgun and capture sequenced data, we obtained mtDNA genomes from 5.8× to 1,703.5× (median depth of 327.65×) (). To further determine the species affiliation, we analyzed these bovid mtDNA sequences together with the currently known matrilineal variability of nine modern Bovini species (taurine cattle [B. taurus], zebu [Bos indicus], gaur [B. gaurus], gayal [Bos frontalis], banteng [Bos javanicus], yak [Bos grunniens], American bison [Bison bison], European bison [Bison bonasus], and water buffalo [Bubalus bubalis]) (). Both the constructed maximum likelihood (ML) and the Bayesian phylogenetic trees confirmed the placement of the 10 SNSZ bovids as sister lineages to the modern gaurs and gayals. (Fig. 2 ). The 20 gaur and gayal mtDNA sequences clustered into three divergent haplogroups, which we named G1, G2, and G3. Moreover, we found that 9 of the 10 ancient samples were tightly clustered together (G3) and diverged from the living gaur and gayal lineages (G2) ∼18,000 y B.P. (Fig. 2), which is suggestive of a possible extinct sublineage with deep divergence from the modern lineages. However, one individual, H63, was closer to the modern gaur and gayal lineages (G2a) compared to other individuals and showed an ancestry state to the modern lineages, which is an indication of a possible transitional population between the “northern” and “southern” gaurs (Fig. 2). The relatively young ages of G2a and G3 and their star-like features could be connected to an expansion that might have occurred after the Last Glacial Maximum (LGM, 25,000 to 18,000 y B.P.).

Fig. 2.

Identification of ancient tropical animals (gaurs and rhinoceroses) in the margin of the Qinghai–Tibet Plateau. (A) Phylogeny of ancient and modern gaur mitogenomes. (B) Bayesian skyline plot showing the gaur and gayal effective population size trends. The y axis indicates the effective number of females, as inferred from our mitogenome dataset considering a generation time of 6 y. (C) An ML time scale of gaur and gayal mitogenomes. (D) The ML tree confirmed the phylogenetic affinities of present-day gaurs and ancient SNSZ gaurs using 8,961,232 autosomal sites. (E) Phylogeny of the rhinoceros estimated from partial mitochondrial genome sequences, showing strong support for the grouping of ancient SNSZ rhinoceroses with Sumatran rhinoceros.

To corroborate the above species identification and further explore the population relationships, we mapped the genome-capture data from the 10 samples to the taurine cattle reference genome (ARS-UCD1.2). We obtained genomic nuclear DNA data with sequence depths ranging from 0.12× to 3.73× (median 1.16×). We found that these 10 samples included five male and five female samples, respectively (). For the three relatively higher-genomic-coverage samples, H113, H57, and H48, we found that they possessed fourth- to third-generation affinity (), which indicated a relatively small population size of the ancient bovid species. The genomic phylogenetic neighbor-joining (NJ) phylogeny (using 990,246 high-quality autosomal sites) clustered all 10 samples into a single clade with a clear divergence from the modern gaur and gayal populations (). We also constructed the genomic ML phylogenies of Bos species, including two of the relatively higher-coverage aDNA genomes (H113 and H57) (using 8,961,232 high-quality autosomal sites and 72,720 coding region sites). Both the NJ and ML trees confirmed the above mtDNA results and were consistent with the topology of the Bos species (29) (Fig. 2 and ). We further examined the genetic relatedness between ancient SNSZ gaur and modern bovids by performing an outgroup f3 analysis on whole-genome sequence data to measure the shared drift f3(Buffalo; SNSZ, X). Our results corroborated the above ML and NJ tree results and showed that the ancient samples are genetically most similar to modern Asian wild gaur ().

For the two rhinoceros faunal remain samples, we performed shotgun sequencing and retrieved mtDNA data with a depth of 0.42× and 0.94×, covering 5,895 bp and 9,940 bp, respectively. To infer their species affiliation, we analyzed the two samples with the mitogenomes of six rhinoceros species (Sumatran [D. sumatrensis], white [Ceratotherium simum], black [Diceros bicornis], woolly [Coelodonta antiquitatis], Javan [Rhinoceros sondaicus], and Indian [Rhinoceros unicornis] rhinoceroses) and included one mitogenome of a Perissodactyla species (tapir, Tapirus terrestris) as an outgroup. The phylogenetic ML tree clearly proved that the two ancient samples were a sister clade of the Sumatran rhinoceros (Fig. 2), which currently only inhabits the islands of Sumatra and Borneo in Indonesia.

Up to now, there has been no record that gaurs once inhabited relatively high latitudes in northern China, especially in the higher-altitude NETP. Our systematic aDNA dissection revealed that the present-day low-latitude tropical rainforest gaurs and Sumatran-like rhinoceroses once rambled over northern Asia as far north as the margin of the NETP at elevation up to ∼2,300 m a.s.l. ∼5,200 y B.P.

Reconstructing the SNSZ Gaur Population Dynamics and Testing Possible Gene Flow with the Other Bovid Species.

To investigate the population dynamics of the SNSZ gaur population in the NETP, we constructed a Bayesian skyline plot (BSP) reflecting changes in the female effective population size through time using 21 mitogenomes including ancient gaur, modern gaur, and gayal mtDNAs plus one banteng (JN632605) as an outgroup. The BSP reconstruction showed a decreasing tendency in the effective female population size of gaurs and gayals beginning ∼20,000 y B.P., coinciding with the LGM. A more rapid population decline of the gaur population size occurred ∼5,000 y B.P., which persisted for a relatively long time afterward (Fig. 2). The latter decline is congruent with the paleoclimatic record (Qinghai Lake) that the millennial-scale abrupt climate changes with the lowest summer temperature occurred in the NETP during 5,000 to 3,600 y B.P. (30) (Fig. 1), which might have resulted in the shrinking of many local wild animal populations (including the gaurs and Sumatran-like rhinoceros) and indicated distinct or massive escape from the habitat region and southward migration to lower-latitude warm regions.

We applied Patterson’s D-statistic to further investigate the relatedness of the SNSZ gaur with extant bovid species and to test for past gene flow between the SNSZ gaur and other modern bovid species (31). The D-statistics demonstrates possible gene flow from SNSZ gaur into modern gaur (). Furthermore, we did not detect gene flow between SNSZ gaur and yak, indicating possible niche isolation between ancient SNSZ gaur and yak ∼5,200 y B.P. In contrast, we observed more complex introgression from modern gaur, gayal, and SNSZ gaur into Chinese indicine cattle and banteng introgression into Chinese indicine cattle (20) (). Overall, our results are consistent with the absence of yak and domesticated cattle faunal remains at the SNSZ archaeological site (14), suggesting an isolated spatial-temporal ecological niche among gaur, yak, and domesticated cattle in the NETP during the late Neolithic period.

Discussion

The gaur, which also known as the “Indian bison,” is the largest bovid species worldwide with an existing population size of ∼20,000 and is designated as a vulnerable species on the International Union for Conservation of Nature (IUCN) Red List (25, 32). The global distribution of B. gaurus both historically and presently is restricted to the tropical and subtropical regions throughout south and southeastern Asia, including in India, Nepal, Bhutan, Bangladesh, Myanmar, Thailand, China (Yunnan), Laos, Cambodia, Vietnam, and Malaysia (33) (Fig. 1). The northernmost gaur populations are found in the southern foothills of the Himalayas, with a latitude reaching northern latitudes of ∼29° N in the hilly terrain of northeastern India (32, 33) (Fig. 1). In India, the B. gaurus-like faunal remains record includes an undated partial skull and horns, which were found in the order alluvium of the Narbada River in central India (34). Bibos gaurus remains, including skulls, horns, mandibles, and teeth found in south and southwest China, date to the Pleistocene (35). The Sumatran rhinoceros, which is also known as the hairy rhinoceros or Asian two-horned rhinoceros (D. sumatrensis), is the smallest rhinoceros of the Rhinocerotidae family member and is the only extant species of the genus Dicerorhinus. The Sumatran rhinoceros was once distributed throughout a continuous range of rainforests, swamps, and cloud forests in India, Bhutan, Bangladesh, Myanmar, Laos, Thailand, Malaysia, Indonesia, and southwest China; however, currently, there are only five substantial populations (four in Sumatra and one in Borneo), with fewer than 100 individuals left in the wild (36). Currently, yak and cattle are central to the economies of the Tibetan Plateau; therefore, most archaeological bovid remains in the NETP have been traditionally identified as Tibetan yak or cattle. In this study, we provide aDNA evidence for the presence of the gaur and Sumatran-like rhinoceroses in the margin region of the NETP during the late Neolithic (∼5,200 y B.P.). These results suggest that the species geographic range likely extended up to an elevation of ∼2,315 m a.s.l. and a northernmost latitude of ∼34.49°N, updating our knowledge of the prehistoric geographic range of gaurs.

The Tibetan highlands host one of the world’s largest pastoral ecosystems (37) and is home to a wide variety of animals that may have been exploited by early hunters (38). We compiled the archaeological and paleoclimatic data from the surroundings of the Tibetan Plateau and found that rhinoceroses and wild water buffalo appeared frequently throughout northern China during 8,000 to 6,000 y B.P., revealing a relatively warm and moist environment that could have sustained the survival of cold-sensitive Sumatran rhinoceroses. These zooarchaeological data are consistent with the paleoclimatic records of monsoonal China, in which the summer temperature and precipitation reached maximum values during the Holocene (39, 40) (Fig. 3 ). Previous studies of SNSZ report the presence of badgers, snub-nosed monkeys, and Bambusoideae vegetation, all of which indicate a warmer and moister climate ∼5,200 y B.P. (28). Gaur and Sumatran rhinoceroses currently inhabit tropical or subtropical low-latitude regions that have been significantly warmer than the SNSZ in the recent past (∼6.5 °C and ∼556 mm). Additionally, subtropical species, such as badger and wild water buffalo, were among zooarchaeological remains from the Guanzhong area in the central Loess Plateau, ∼500 km to the east of the NETP, supporting higher temperature during the middle Yangshao period (∼6,000 to 5,500 y B.P.) compared to the present temperature (41, 42). A comparison with high-resolution paleoclimate records from neighboring southwest China (Xingyun Lake) and the Loess Plateau (Gonghai Lake) further demonstrated that the summer temperature (mean July temperature) and annual precipitation were likely over ∼2 °C and ∼120 mm higher ∼5,200 y B.P. than in the most recent past (Figs. 1 and 3 ). These climatic conditions could have facilitated the seasonal presence of gaur and Sumatran-like rhinoceroses in the SNSZ region. The optimum hydrothermal conditions during the Holocene occurred between ∼7,500 and 5,500 y B.P., when hunting wildlife was the primary subsistence strategy in the NETP and adjacent regions (Fig. 3). Additionally, rhinoceroses were also hunted in the western Loess Plateau during this period (Fig. 3). This highlights the intimate connection among climate change, the geographical distribution of wildlife, and human hunting behavior in northwest China during that period.

Fig. 3.

The spatiotemporal distribution of unearthed remains of typical bovids and rhinoceroses compared with paleoclimate records, zooarchaeological data, and the number of dating sites. (A) Pollen-based mean July temperature reconstructed from Xingyun Lake (red line) (39). (B) Pollen-based annual precipitation reconstructed from Gonghai Lake (blue line) (40). (C) Variations of the percentage of NISP of wild mammals along the Yellow River basin, including the NETP, Loess Plateau, and North China Plain, during 8,000 to 3,000 y B.P. (D) Number of dated sites in the NETP and adjacent regions from 8,000 to 3,000 y B.P. (E) The historical distribution of typical bovids and rhinoceroses from 8,000 to 3,000 y B.P. through China (below 38° N) (67).

The intensification and expansion of rain-fed agriculture ∼6,000 y B.P. altered the subsistence strategy in the Loess Plateau during the late Yangshao period, including the remarkable enhancement of significant animal husbandry that has dominated animal exploitation in this area since then (43). However, in the marginal areas of the mountainous environment, low human populations, combined with favorable water heating to ensure that the area became an ideal habitat for various wildlife including gaur and Sumatran-like rhinoceroses. The rich biodiversity in SNSZ and other regions of the NETP during the late sixth millennium facilitated the exploration of hunting, and therefore the emergence of a unique subsistence strategy dominated by millet cultivation and game hunting. This lifestyle was evidently different from that in the middle-lower valleys of the Yellow River, which was based on millet cultivation and animal husbandry (44).

The diversity of wildlife in northern China declined during 5,000 to 4,000 y B.P. in comparison to the former millennium (Fig. 3); however, hunting game continued to be an important strategy in the NETP above 2,500 m a.s.l. (9, 12) and probably contributed to the human influence on the local environment, such as the exploitation of woodlands during this period (45). The temperature and precipitation rapidly declined between 5,000 and 4,000 y B.P. (30, 39, 40). This climatic warm–humid/cold–dry transition was likely responsible for the rapid gaur population decline with an estimated onset of ∼5,000 y B.P. (Fig. 2) and for the change in human behavior with increased reliance on domesticated bovids. No remains of tropical wildlife including gaur and Sumatran-like rhinoceroses were identified from sites during this period in northern China, while gaur remains were reported from later southern sites (3,000 y B.P.) in the Yunnan province of southwestern China (46), suggesting that climate deterioration might have triggered the southward migration of cold-sensitive wild animals, which vacated niches for domesticated animals in the NETP (Fig. 3).

We also demonstrate that the wild animal distribution along the Yellow River basin, including the NETP, Loess Plateau, and North China Plain, was associated with human activity (Fig. 3). During 5,000 to 4,000 y B.P., although hunting game remained in northern China, the increasing human activity intensified the shrinking of wildlife populations, extinctions, and southward migrations and urged the transition of human subsistence strategies from hunter-gather to agricultural production. Based on the radiocarbon dates from archaeological sites in the NETP and the adjacent region, we found that the intensity of human colonization in that area showed an upward trend between 6,000 and 3,000 y B.P. (Fig. 3). During 4,000 to 3,000 y B.P., the intensity of human settlements significantly increased (Fig. 3). Livestock almost dominated animal exploitation in all regions of human settlements including the Tibetan Plateau, where there are vast pastures for domestic herbivorous animals, such as yaks, cattle, sheep/goats, and horses (9), which significantly reduced the habitats of wild animals. Thereafter, gaur and Sumatran-like rhinoceroses nearly vanished from areas north of ∼34°N.

In conclusion, our genome-wide aDNA analysis reveals that the wild bovids in the late Neolithic in the NETP included B. gaurus, which has since been obscured by the distribution of modern yaks and cattle. Our findings push the historical gaur distribution from ∼29°N to ∼34°N. Additionally, by integrating aDNA with paleoclimatic and zooarchaeological evidence, our interdisciplinary investigation showed that the regional mid-Holocene warm climate period facilitated the low-latitude tropical animals undertaking northward expansion, enriching the biome diversity in the NETP, providing abundant hunting resources for local early human societies, and facilitating the exploration of hunting uniquely combined with millet cultivation. However, the subsequent climatic abruption and intensified human activities during 5,000∼4,000 y B.P. narrowed the spaces for wildlife. Gradually, hunting game was replaced by herding activities after 4,000 y B.P. Indeed, after 3,600 y B.P. humans extensively and permanently settled on the high-altitude areas with cold-tolerant agropastoral economies (9). Both climatic change and human activities contributed to the disappearance of gaur and Sumatran-like rhinoceros from the middle and high latitudes, eventually shaping their modern tropical biogeographic distributions. Additionally, our aDNA analysis supports the rapid turnover of different bovid species in the NETP during the mid-Holocene. Further aDNA studies, combined with multidisciplinary approaches from more archaeological sites in northern Asia are necessary to further our understanding of the past interactions among climate change, biogeography, and human activity.

Materials and Methods

Bovidae Bone Remains Sampling from the SNSZ Archaeological Site.

The archaeological site of SNSZ dates to ∼5,280 to 5,050 cal y B.P. and is located at an altitude of 2,323 m a.s.l. in the margin region between the NETP and the Loess Plateau of northwestern China. The composition of faunal remains at this archaeological site is diverse (). Species identified via morphology include rabbits (Lepus capensis), bamboo rats (R. pruinosus), Siberian roe deer (Capreolus pygargus), barking deer (Muntiacus aurus), red deer (Cervus elaphus), sika deer (Cervus nippon), Bovidae, badgers (Meles meles), pigs (Sus scrofa domestica), Chinese goral (Naemorhedus goral), and Rhinocerotidae (14). The specific species identification of Bovidae and Rhinocerotidae specimens via morphological approaches has not been possible. Here, we randomly selected 10 Bovidae and 2 Rhinocerotidae samples () to perform molecular identification.

Radiocarbon Dating.

We collected five different animal bone fragments in total for the AMS radiocarbon dating, and one Bovidae sample of them was also used for aDNA extraction and species identification (Table 1). Two animal collagen samples from one deer and one gaur were extracted and pretreated at the Chronology Laboratory of the Ministry of Education Key Laboratory of Western China’s Environmental Systems, Lanzhou University. Then, these samples were sent to Peking University in Beijing for the dating measurement. The other three samples from two gaurs and one rhinoceros were dated at Beta Analytic, Miami, CA. All chronological data were calibrated by the Calib radiocarbon calibration program (v.7.0.2) (47) and the IntCal13 calibration curve. All reported ages are relative to AD 1950 (referred to as “cal y B.P.”)

aDNA Extraction, Library Preparation, Genome Capture, and Sequencing.

We collected 12 bone samples to perform molecular identification of these faunal remains from the SNSZ site (). aDNA extraction, the quantification of DNA preservation, and double-stranded DNA library (DSL) construction were performed following previous reports (48, 49) in the aDNA dedicated laboratory of the CGG at Kunming Institute of Zoology, Chinese Academy of Sciences. In brief, prior to DNA extraction we conducted ultraviolet (UV) irradiation on the samples for ∼2 h and removed a surface layer using a sterile dentistry drill. For each sample, we drilled ∼200 mg of bone powder for aDNA extraction. aDNA extraction was performed using the silica adsorption protocol as in previous reports (48, 49), which resulted in 200 µL of extraction for each sample.

For the screening of DNA preservation we used 10 µL of extract per sample to build DSLs, using a commercial kit from New England Biolabs (50). We performed the PCRs with a final volume of 50 µL using KAPA uracil+ polymerase (Kapa Biosystems), which can incorporate nucleotides across uracils and preserves the deamination-induced damage patterns that are an indication of authentic aDNA (51). We sequenced the constructed DSL aDNA libraries (including all later captured genomic libraries) using the Illumina Hiseq X Ten (PE-150) platform (WuXi NextCode Genomics Co., Ltd.). For the 10 Bovidae samples, we treated the remaining extract with uracil-DNA-glycosylase to remove uracil residues from the DNA to leave abasic sites and applied enzyme endonuclease VIII (Endo VIII), which cuts the DNA at the 5′ and 3′ sides of abasic sites. Following the treatment, DSLs were then built as described above. We then performed genomic capturing to these libraries using RNA baits transcribed from cattle genomic DNA as in previously reported methods (52, 53). The PCRs of both precapture and postcapture samples was conducted using the KAPA HiFi Hot Start Polymerase protocol (Kapa Biosystems).

Process of Ancient Data Mapping, Alignment, Filtering, and mtDNA Assembly.

Adapter sequences of the paired-end reads were identified and removed using AdapterRemoval version 2.2.0 (54). Reads were separately aligned to mitogenome references of the taurine cattle (V00654.1) and gaur (NC_024818.1) using Burrows–Wheeler Alignment (BWA) (55). To improve the alignment to the circular mtDNA genome, a 30-bp sequence from the end was attached to the beginning. Reads were subsequently filtered for PCR duplicates and a mapping quality <30 using SAMtools (v. 0.1.19) (56). Four resequenced modern gaurs and five resequenced modern gayals were aligned to the same B. gaurus mitogenome reference (). For the two rhino samples, we mapped the reads to the R. unicornis mitogenome reference (NC_012684.1). mtDNA coverages were calculated by Qualimap (v. 2.2) ().

We retrieved the mtDNA data of ancient gaur with the sequence depths range from 9.05× to 1,703.55× () and the coverage range from 5,895-bp to 9,940-bp mtDNA for the two Rhinocerotidae samples. The average length of the Bovidae and Rhinocerotidae molecules were 51 to 80 bp and 75 to 77 bp, respectively, the TDR of all of the shotgun sequenced samples were >10%, and the maximum mapping rate was 3.75% (T0109), all of which indicated typical aDNA characteristics (). The final mtDNA BAM files were used with Mapping Iterative Assembler v. 1.0 to assemble mtDNA consensus sequences and manually and visually inspected for all mutations against the sequence.

Alignment of Ancient Bovid NGS Genomic Reads to the B. taurus Genome.

All raw reads were processed through a standardized pipeline adapted for the ancient samples. The adapter sequences of the paired-end reads were identified and removed using AdapterRemoval version 2.2.0 (54). We merged read pairs for which the expected index was observed, required an overlap of at least 11 bp. Because there is no B. gaurus reference genome, the bovine reference genome (ARS-UCD1.2) with the addition of the Y chromosome from BosTau7 was used as a reference. Alignments were performed using the Burrows–Wheeler Alignment Tool (55) (v. BWA 0.7.5a-r405) with the parameter aln -l 16500 -n 0.01, and the parameter samse with the option -r to define the read groups and to produce unfiltered SAM files, which were subsequently converted to BAM files using the command samtools view -Sb. The BAM files were sorted using SAMtools (v. 0.1.19) (56). The BAM files were merged using picard-tools (v. 1.129), and indel realignment was performed using the Genome Analysis Tool Kit (GATK; v. 3.8) (57). Duplicates were removed using Picard-tools (v. 1.128, broadinstitute.github.io/picard/) (MarkDuplicates REMOVE_DUPLICATES = true) prior to indel realignment. The samples were filtered using SAMtools to remove unaligned reads (view -F4) and for mapping the quality of Q25 (view -q25). Quality control parameters of the alignment sequencing data and coverage statistics were calculated using Qualimap (v. 2.2) (58).

Alignment of Modern Bovid NGS Genomic Reads to the B. taurus Genome.

We generated four whole genomes of B. frontalis and combined these data with previously published genomes from modern bovids for phylogenetic analysis and species identification (), including four B. taurus, four B. indicus, three B. javanicus, five B. frontalis, four B. gaurus, four B. grunniens, one B. bison, two B. bonasus, and two B. bubalis. For the modern gaur and gayal data, all animals were assigned to five geographical groups following their sources: American zoo gaur, European zoo gaur, Asian wild gaur, Southeast Asian wild gayal, and Chinese gayal (). These data were used for “high confidence” variant datasets.

We further combined other 31 indicine cattle from South China and South Asia to investigate the potential pattern and consequences of gene flow between ancient species and modern species (). All passed quality-filtered reads were aligned to the B. taurus reference genome ARS-UCD1.2 with the addition of the Y chromosome from BosTau7 using bwa mem. Duplicated reads were filtered using Picard-tools (v. 1.128). Reads around indels were realigned by GATK (v 3.8) (57).

Variant Discovery and Genotyping.

Bcftools mpileup and call were used to call variants from the “high confidence” sample set (-q 30 -C 50 -Q 20 -B -a FORMAT/AD, FORMAT/SP, FORMAT/ADF, FORMAT/ADR) and generate vcf files (-O z -m -o). Indels and any variants within 3 bp near indels were removed using the bcftools filter. Repeat regions as defined by the UCSC Browser and Repeat-Masker files were removed (59) (genome.ucsc.edu/). Tri- and quad-allelic sites were removed. Variants with coverages of twice the mean coverage that were marked as missing (“./.”). Heterozygous variants present in a single individual or more than 75% of individuals were removed. Variant positions with missing data in any individual were then removed, resulting in a dataset of 45,964,584 sites.

For samples H113 and H57 with >2× mean coverage, the above defined “high confidence” sites were called. The sites were called using bcftools mpileup with the same options but without recalibration (bcftools mpileup -B) and without filtering for variant sites with bcftools. After indels and sites within 3 bp of the indels were removed, tri- and quad-allelic sites were removed. Sites were set to missing as defined above (greater than twice the mean coverage). Individuals were then pseudodiploidized by randomly sampling a read at each site and setting that individual as homozygous for that allele. Variant positions with missing data in any individual were then removed. All transitions were removed. This call set was then merged with the Moderns set to create an “Ancients and Moderns” dataset, resulting in a dataset of 8,961,232 sites. This dataset was used for species identification.

For samples with less than 2× depth, no minimum coverage filter was imposed and a maximum of 4 read coverage per site was permitted. Individuals were then pseudodiploidized by randomly sampling a read at each site and setting that individual as homozygous for that allele. Variant positions with 50% missing data in 10 ancient samples were then removed. All transitions were removed. This call set was then merged with the “Ancients and Modern” set to create a “Low Coverage Ancients and Modern” dataset, resulting in 990,246 sites. This dataset was used for species identification for all SNSZ samples.

We combined “Ancients and Moderns” and 31 other previously published indicine cattle genomes to investigate the gene flow between ancient genomes and modern genomes. We used the above criteria to genotype and filter single-nucleotide polymorphisms (SNPs), resulting in 4,893,675 sites.

Phylogeny Construction and Demographic Inferences with Gaur Mitogenome Sequences.

Gaur mitogenome sequences were aligned using ClustalW v2 (60). The control region, as well as any positions in the alignment that contained missing data, were removed. The resulting alignment was manually annotated using gaur (GenBank: NC_024818.1) as the reference. The consensus mitochondrial genomes of the 10 ancient gaur sequences, four B. gaurus and five B. frontalis genomes (newly assembled) were aligned with the available mitogenomes of six Bovini species, including the mitogenomes of one B. taurus, one B. indicus, one B. javanicus, one B. grunniens, one B. bison, one B. bonasus, and one B. bubalis (.

BEAST analysis was performed with BEAST (2.6.1) (61). To calibrate the molecular clock, we built a tree by including ancient gaur, modern guar, and gayal mtDNAs plus one banteng (JN632605) as an outgroup. The coalescence ages of the ancient gaurs (5.2 ± 0.15 ky; 95% CI: 4.94 to 5.46 ky) were used as an internal (recent) calibration point. The major haplogroups were considered as monophyletic to calculate their age estimates. We then obtained a BSP from the gaur phylogeny by running 50 million iterations with samples drawn every 1,000 steps.

Rhinoceros mitogenome sequences were aligned using ClustalW v2 (60). The resulting alignment (3,452 bp in length) was manually annotated in Geneious using the modern Sumatran rhinoceros (GenBank: NC_012684.1) as a reference. For interspecies comparison between ancient SNSZ and modern rhinoceros, mitogenomes for Indian, Javan, Sumatran, woolly, and white rhinoceros were downloaded and aligned with the two SNSZ specimens using ClustalW v2. Alignment columns containing missing data were not considered to enable direct comparison of genetic distances within extant species with those estimated from SNSZ assemblies. The alignment contained five rhinoceros subspecies, including one R. unicornis, one R. sondaicus, one C. simum, one C. antiquitatis, and one D. sumatrensis (). The GTR+G model of nucleotide substitution was selected by comparing the Bayesian information criterion scores in ModelGenerator v0.85 (62). An ML tree of 3,452 bp length alignment was computed using MEGA7.0 with 1,000 bootstrapping replicates () (63).

Phylogenetic Inference Using a Supermatrix of Nuclear Coding Sequences.

An NJ phylogenetic tree of 10 SNSZ samples and other species from the Bos genus was also constructed using the “Low Coverage Ancients and Modern” dataset with 990,246 sites.

The phylogenetic relationships of two ancient samples (H113 and H57) and nine modern species were inferred by the ML method with 100 bootstraps and the NJ method using the “Ancients and Modern” dataset. The low-coverage genomes of other SNSZ gaurs were not included in this analysis. We used all 8,961,232 sites and variants within coding sequences (CDSs) separately (). The protein-coding gene annotation information was performed according to the annotation of the current cattle reference ARS-UCD1.2. We obtained 72,720 sites within the CDSs. ML phylogenetic inference was performed using RaxML v8.1.3 (64) using default parameters. The resulting tree was rooted using water buffalo (B. bubalis).

A total of 8,961,232 autosomal SNPs was selected to construct an NJ tree including the two high-coverage samples with PLINK using the matrix of pairwise genetic distances ().

We further examined the genetic relatedness between SNSZ species and modern gaur and gayal by performing the outgroup f3 statistic to measure the shared drift using the ADMIXTOOLS (version 6.0) package and qp3Pop program. The test (Buffalo; SNSZ gaur; X) was performed utilizing the “Ancients and Modern” dataset, rotating X through 5 modern gaur and gayal populations.

Sex Determination and Individual Genetic Relatedness for Bovid Faunal Remains from the SNSZ Site.

The ratio of mapped reads to autosomes and X and Y chromosomes could be used to infer sex (65). We used this method to infer the sexes of the 10 samples with 990,246 sites, and we found that the sets included five male and five female samples, respectively (). For the three relatively higher genomic coverage samples H113, H57, and H48, we analyzed the kinship using KING (66) to assess their relatedness, and we found that they possessed fourth- to third-generation affinity (), which indicated a relatively small population size of the ancient bovid species.

Introgression Events Test among SNSZ Gaur with Other Bovid Species.

ADMIXTOOLS (version 6.0) software was used to determine relationship between the ancient SNSZ gaur and modern bovid groups using 4,893,675 autosomal sites without transitions. We used Patterson’s D-statistic to test whether ancient SNSZ gaur shared more derived alleles with modern bovids using the qp-Dstat program ().

Compilation of the Archaeologically Excavated Large Animals, Paleoclimate, and Human Activity Records from Northern China.

We created a database of the rhino, gaur, water buffalo, yak, and cattle presence through northern China (below 38° N). The whole database and references to the original studies are available in . The pollen-based mean July temperature reconstructed from Xingyun Lake and annual precipitation reconstructed from Gonghai Lake were also collected. The NISP of wild mammals and radiocarbon dates along the Yellow River basin from 8,000 to 3,000 y B.P. were collected (67) ().