Bovine bone remains were collected from the Taosi and Guchengzhai sites. Both sites are located at the lower and middle reaches of the Yellow River drainage basin (Figure 1) and date back to the Longshan Culture period (approximately 4,000 YBP) (Supplementary Table S1). Specifically, the Taosi site is a large urban site located between the Ta’er Mountain and the Fen River, approximately 7.5 km northeast of Xiangfen County, Shanxi Province. The Taosi site dates back 4,350–3,900 years (Brunson et al., 2016a). The paleoenvironmental conditions around the Taosi site were warmer and wetter than now (Wang et al., 2011). Many domestic animals remain, including cattle bones, were found at the Taosi site; these animals were kept primarily as animals of the wealthy inhabitants (Brunson et al., 2016a). The Guchengzhai site is located in Dafanzhuang Village, 35 km southeast of Xinmi City, Henan Province, 1.5 km northeast of the intersection of the Zhenshui River and the Yan River. The Guchengzhai site dates back 4,150–3,950 years. Similarly, the paleoenvironment of the Guchengzhai site was slightly warmer and more humid than now (Xu et al., 2013). Domestic animals, including dogs, pigs, cattle, and sheep, as well as wild animals such as bears, sika deer, and elks were extensively discovered at the Guchengzhai site (Guo, 2018). The detailed information on all samples were shown in Supplementary Table S1. Before aDNA analysis, all samples were morphologically determined to be cattle remains.

All aDNA experiments were conducted in the Ancient DNA Laboratory at China Agricultural University and repeated independently in the Animal Genetic and Evolutionary Laboratory at Foshan University, following the experimental guidelines to eliminate contamination with modern samples. First, the equipment involved in the experiment was soaked with 5% (v/v) sodium hypochlorite solution and washed with 75% (vol/vol) alcohol, followed by UV irradiation for 1 h. Thereafter, the samples were prepared by cautiously removing the adhering soils on the surface of samples with a drill, followed by washing with 5% (vol/vol) sodium hypochlorite solution and then with double-distilled water and subsequent UV irradiation. The bones were then ground to fine-grained powder using an automatic sample quick grinding machine (Jingxin Industrial Development Co., Ltd., China), and 100–200 mg of powder was subjected to DNA extraction using a silica adsorption column system. At China Agricultural University, 50 μL of extracted DNA was obtained using the QIAamp DNA Investigator Kit (QIAGEN, Germany). At Foshan University, aDNA extraction was performed following the method described by Dabney (Dabney et al., 2013) with some adaptive modifications, including setting the lytic temperature to 45°C and extending the time of incubation with double distilled water (ddH₂O) to 10 min.

The mitogenomes of the genus Bos, including Holstein (Bos taurus), zebu (Bos indicus), gayal (Bos frontalis), and yak (Bos grunniens), were obtained using long-range PCR amplification. The PCR products were then used to prepare Bos hybrid baits. A total of 50 ng DNA per sample was used to construct a double-stranded DNA library using the KAPA Hyper Prep Kit for Illumina® platforms (KAPA Biosystems, San Francisco, USA) with minor modifications. The libraries were quantified using a Qubit 4.0 Fluorometer (Life Technologies Corporation, USA) and submitted for subsequent mitogenome enrichment and capture. Four libraries of Taosi samples were equally enriched to a total of 1 ng of DNA with unique adapter index sequences. Thereafter, mitogenome enrichment and hybrid capture were performed using the xGen Hybridization and Wash Kit (Integrated DNA Technologies, Inc., USA). The enriched library was then sequenced with 150 paired end reads (PE150) using the Illumina HiSeq X Ten platform. For the resequencing samples, up to 100 ng DNA per sample was used to prepare the library with the VAHTS Universal Plus DNA Library Prep Kit for MGI (Vazyme Biotech Co., Ltd, China). The qualified libraries were sequenced with PE100 using the MGI platform.

The quality of the raw data was assessed using the FastQC v0.11.9 (Andrews, 2010). Adapter and low-quality reads were filtered using clip&merge v1.7.8 from EAGER pipeline (Peltzer et al., 2016) with the following parameters: -minlength 25, --minquality 30. Before mapping, the first 250 and 500 bp of the taurine cattle reference mitogenome (GenBank accession No.: V00654) were concatenated to the end using CircularGenerate from the EAGER pipeline (Peltzer et al., 2016). The filtered reads were mapped to the circularized reference V00654 using CircularMapper v1.0 (Peltzer et al., 2016), with parameters as -n 0.01, -l 16,500. The resulting file was converted to a bam file and sorted using Samtools v1.9 (Li, 2011). Unmapped reads and reads with mapping quality <25 were removed using Samtools v1.9. Duplicate reads were removed using MarkDuplicates of Picard-tools v2.22.9 (http://broadinstitute.github.io/picard/). The endogenous DNA, coverage rate, and mapping quality distribution were calculated using Qualimap v2.2.2-dev (Okonechnikov et al., 2016). The authenticity of the retrieved sequences was evaluated according to the characteristics of the aDNA (Hofreiter et al., 2001; Willerslev and Cooper, 2005) using the mapDamage software (Jonsson et al., 2013).

Three methods, including Mapping Iterative Assembler (MIAv.1.0) (Green et al., 2008), BCFtools (Li, 2011), and UnifiedGenotyper in the Genome Analysis Tool Kit (GATK, v3.5) (DePristo et al., 2011), were comparatively deployed to generate the mitogenome sequences. For MIA, the filtered bam files were converted to Fastq format using BEDtools (Quinlan and Hall, 2010). The program ma in MIA was used to generate complete sequences. The latter two software packages generate the complete sequence based on the mutation sites relative to the reference sequence using BCFtools consensus and manual generation, respectively. Finally, the complete mitogenome sequence of the ancient samples was obtained by comparing the differences in the complete sequences from different assembly methods and the filtered bam files.

Variant calling was carried out using BCFtools v1.8 mpileup (-d 5000 -Q 20 -a DP, AD, ADF, ADR, INFO/AD, INFO/ADF, INFO/ADR, SP, DV, DP4, DPR, INFO/DPR). Single nucleotide polymorphism (SNP) sites with <3× coverage depth were filtered out. UnifiedGenotyper in GATK3.5 was also used for the variation calling, whereas SNP sites with <3× coverage were discarded. Subsequently, the filtered SNP sites were manually checked and corrected according to the bam files visualized on the Integrative Genomics Viewer (IGV) (Thorvaldsdottir et al., 2013). Eventually, common SNP loci were retained as the true mutation loci.

All the captured and re-sequenced reads for the ancient samples were aligned to the reference mitogenome (GenBank accession No.: V00654) using the online MAFFT software with default parameters (Katoh et al., 2019). Positions embracing missing data or uncertain sites were removed using MEGA 7 (Kumar et al., 2016). Then the ratio of nonsynonymous and synonymous substitutions (dN/dS) was also calculated using MEGA 7. The 3 ancient consensus sequences were aligned to 36 others extant Bovinae mitogenomes from GenBank and called Dataset A (Supplementary Table S2), which included 1 Bison, 1 Bison bonasus, 1 Bison schoetensacki, 1 Bison priscus, 1 Bubalus bubalis, 1 Bubalus depressicornis, 1 Bubalus carabanensis, 1 Bos frontalis, 1 Bos gaurus, 1 Bos grunniens, 1 Bos mutus, 1 Bos javanicus, 3 Bos primigenius, 2 Bos indicus, and 19 Bos taurus. To further identify the haplogroup of the ancient samples, 510 cattle including Bos taurus and Bos indicus were aligned with ancient samples and called Dataset B (Supplementary Table S3).

To construct phylogenetic trees and infer population structure, the nucleotide substitution model was selected using jModelTest 2.1.1 (Darriba et al., 2012), based on the model selection strategy of Bayesian information criteria. Bayesian phylogenetic trees were constructed using MrBayes 3.2.7 (Ronquist et al., 2012), with the following parameters: 20,000,000 generations of Markov Chain Monte Carlo (MCMC) chains, checking every 10,000, sampling every 2,000, and discarding the first 10% as burn in. The consensus tree was constructed using FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).

Principal component analysis (PCA) was used to further verify the genetic affinities of the ancient samples among the genus Bos using the R Package “adegenet” (Jombart, 2008). Part of dataset A was interrogated to infer demographics using BEAST 2.6.3 (Bouckaert et al., 2014). The ages of the ancient cattle were set at 4,000 YBP, which was the time of the culture in the respective sites of origin. We used Strict Clock as the clock model and Coalescent Bayesian Skyline as a distribution of the tree priors. The MCMC chains were set to 1000,000,000 generations. The log file was checked using Tracer v1.7.1 (Rambaut et al., 2018), and the Bayesian Skyline Plot (BSP) was obtained using Tracer v1.7.1. The optimal nucleotide substitution model was derived using MEGA 7 estimation for the maximum likelihood (ML) tree. Thereafter, an ML tree was built using MEGA 7 with 1,000 bootstrapping replicates and a timetree inferred using the Reltime method (Tamura et al., 2012). The timetree was computed using two calibration constraints from TimeTree (Kumar et al., 2017), including the time of divergence between Bos taurus, Bos primigenius, Bos taurus, and Bos indicus.

T haplogroups of Dataset B was used to calculate population pairwise F _ST values, the gene flow (Nm) and the net genetic distance using DNAsp 5.10 (Librado and Rozas, 2009). Meanwhile, numbers of haplotypes and variable sites, haplotype diversity (Hd), nucleotide diversity (Pi), and the average number of nucleotide differences (K) were estimated using the same data and software. Besides, T haplogroups of Dataset B was used to conduct PCoA (Principal Coordinate Analysis) among populations. The neighbor net (Bryant and Moulton, 2004) was constructed using splitstree5 (Huson and Bryant, 2006) based on Fst value. PCoA analysis was performed using the OmicStudio tools at https://www.omicstudio.cn/tool. The final result was graphed using R Package “ggplot2” (Hadley, 2016).

A total of eight ancient samples, including four bovines (4,350–3,850 YBP) from the Taosi site and four bovines (4,150–3,100 YBP) from the Guchengzhai site, were subjected to DNA extraction. All qualified aDNA extractions were selected for shotgun sequencing. However, owing to the typically low levels of endogenous DNA in these samples, only one sample from the Guchengzhai site (GCZ4C) obtained sufficient endogenous reads (0.0001%) to yield a near-complete mitogenome. Therefore, the remaining samples were subjected to mtDNA enrichment and hybridization sequencing using Bos hybrid baits. Total sequences of the near-complete mitogenome for the two samples from the Taosi site (TS1C and TS2C) were retrieved after enrichment using baits. After capture sequencing, the endogenous mtDNA read ratio increased to 0.3599% and 0.0474% for TS1C and TS2C, respectively (Table 1), resulting in the recovery of three complete ancient mitogenomes for late Neolithic archeological sites in the Middle Yellow River region. The average length of hybrid capture reads was approximately 140 bp, and the average length of resequencing reads was approximately 57 bp (Supplementary Figure S1). Damage pattern analyses showed a high C–T mutation at the 5′ end and a high G–A mutation at the 3′ end (Supplementary Figure S2), showing typical damage and fragmentation patterns characteristic of endogenous aDNA.

By considering the effect on the integrity and accuracy of ancient mitogenome assembly, we combined three methods including GATK, BCFtools and MIA to generate mitogenome sequences. The number of duplicated sequences varied from 5.50% to 21.6% (Table 1). As a result, we obtained a 16,299 bp (coverage 99.43%) mitogenome sequence for TS1C, a 15,054 bp (coverage 92.13%) sequence for TS2C, and a 14,764 bp (coverage 90.21%) for GCZ4C (Table 1). SNP sites with ≥3× coverage were used in subsequent analyses (Supplementary Table S4).

To determine the species of the ancient bovine samples, the near-complete mitogenome sequences were combined with 36 extant Bovinae mitogenomes, including the genus Bison, Bubalus, and Bos (Supplementary Table S2). The phylogenetic trees illustrated all three ancient samples belonging to the genus Bos rather than Bison or Bubalus, specifically close to the species Bos taurus (Figure 2 and Supplementary Figure S3).

The mutation sites for each sample were discovered by alignment against the taurine cattle reference mitogenome (GenBank accession no.: V00654) (Supplementary Table S4). All ancient mitogenomes contained 13 coding genes, of which COX3 had the most variable sites and the most non-synonymous SNPs. The values of dN/dS of all genes except COX3 and ND5 were lower than 1 (Supplementary Table S5), suggesting potential purifying selection on the mitochondrial functions in the ancient bovine in North China. Based on the key SNPs in the mitogenome, all three ancient samples were classified as haplogroup T3. Specifically, TS1C from the Taosi site belonged to haplotype T3k, while another sample, TS2C, from the Taosi site and GCZ4C from the Guchengzhai site belonged to haplotype T3n (Supplementary Table S4).

To further investigate the phylogenetic relationships among the approximately 4,000-year-old cattle in North China, a Bayesian phylogenetic tree was inferred with 510 extant mitogenome sequences of Bos taurus and Bos indicus around the world (Supplementary Table S3). The Bayesian tree showed that all ancient cattle were on the T3 branch (Supplementary Figure S4). PCA (Supplementary Figure S5A and S5B) and PCoA (Supplementary Figure S6A) supported the phylogeny of the T3 haplogroup of the ancient samples. The population pairwise F _ST values, the gene flow, and the net genetic distance, illustrated that the ancient Chinese cattle had the closest genetic relationship to Chinese taurine cattle, followed by populations from East Asia, Europe, West Asia, America, and Africa (Table 2). Further subdivision of Chinese taurine cattle according to region showed that population in T haplogroups had high haplotype and nucleotide diversity (Table 3), and the ancient Chinese cattle were more closely related to the south Chinese taurine cattle than the north taurine Chinese cattle (Table 2 and Supplementary Table S6). Moreover, PCA (Figure 3) and PCoA (Supplementary Figure S6B) analysis of ancient samples and T3/T4 cattle suggested that TS1C was genetically close to modern Bos taurus from Southern China, while TS2C and GCZ4C were closer to cattle from Europe and East Asia.

To assess the population dynamics of taurine cattle in China, a BSP was produced using 25 mitogenomes. The BSP points to a rapid decrease in the female effective population size approximately 4.65 kya and a steep increase in population approximately 1.99 kya (Figure 4). To infer the divergence time between ancient and modern cattle, a ML tree was constructed with 3 Bos primigenius and 25 Bos taurus cattle, using 1 Bison bonasus as the outgroup (Figure 5). Concerning the differentiation of the major branches in the phylogenetic tree, the T3 and T4 sub-lineages possibly emerged approximately 7.44 kya, and at this time, the haplotype of sample TS1C appeared among the early domesticated cattle. TS2C, another sample from the Taosi site, was divergent to modern populations approximately 4.10 kya, in accordance with the period of rapid decrease in the female effective population size. GCZ4C was divergent to the modern T3n taurine population approximately 0.86 kya, suggesting possible long-term dilution of the ancient northern Chinese cattle.