Introduction

Domestication has been the dominant and most enduring innovation of the transition from a hunter-gathering lifestyle to farming societies, representing the direct exploitation of genetic diversity of wild plants and animals for human benefit. Ancient DNA (aDNA) has proved crucial to understanding the domestication process and the interaction between domesticated species and their wild relatives both within domestication centers and throughout the regions that the domestics expanded into (1–13). The origins of the European domestic taurine, Bos taurus, are located in the Fertile Crescent (14, 15) and unlike dogs, pigs and goats, where the wild forms are still extant, the wild cow (the aurochs) went extinct in 1627. Aurochs, B. primigenius, was present throughout much of Eurasia and Africa before the expansion of domestic cattle from the Levant that accompanied the first farmers during the Neolithisation of Europe. Upon arrival, these early incoming domesticates inevitably coexisted with their wild counterparts in great parts of Europe facilitating gene flow in both directions. In general, taxa within the genus Bos can hybridize and produce fertile offspring (16) which may have facilitated and contributed to domestication, local adaptation and even speciation (5, 17–19). Mitochondrial DNA studies have previously indicated gene flow between domestic cattle and aurochs outside their domestication center (20–24) and more recently, genomic studies have shown the presence of European aurochs ancestry in modern taurine cattle breeds (25, 26). Although cattle have represented a significant economic resource and a prominent cultural icon for millennia, and have been studied through aDNA for more than a decade (20, 21, 27), to date only a few aurochs’ genomes have been sequenced (5, 25). This limits our understanding of the interaction of early cattle herds and wild aurochs, the latter’s local input into modern domestic breeds and human management of these animals in the past.

Aurochs have been widely exploited by humans since the European Palaeolithic, nevertheless, archaeological evidence indicates that the species survived in Europe until historical times. In Iberia, the most recent evidence is found at a Roman site in the Basque Country (28). Domestic cattle were introduced into Iberia with the Mediterranean Neolithic expansion and reached the northern coast of the peninsula around 7000 years cal BP (29). Consequently, aurochs and domestic cattle have coexisted in Iberia for about five millennia. Since then, cattle have played an important role in Iberian societies as a source of food and labour as well as cultural events such as bullfighting. Currently, there are more than 50 bovine breeds officially recognized in the Iberian Peninsula including the Lidia breed, a primitive, isolated population selected for centuries to develop agonistic-aggressive responses with the exclusive purpose of taking part in such socio-cultural events (30). Recently, it has been reported that Lidia breed individuals have the largest brain size among a comprehensive data set of European domestic cattle breeds and are the most similar to wild Aurochs (31). The combination of aggressiveness and larger brain size in the Lidia breed may suggest a higher proportion of aurochs ancestry compared to other cattle breeds.

Here, we present the genomes and stable isotope data of Iberian Bovine specimens ranging from the Mesolithic into Roman times from four archeological sites (Fig. 1A). We explore the extent of interbreeding between wild aurochs and domestic cattle over time and the correlation of genetic ancestry with metric identification and ecological niches. Finally, we compare the results to genomic data obtained from modern Iberian cattle breeds to estimate the genetic contribution of the now-extinct aurochs to the Iberian farming economy.

Figure 1.

Results

Exploratory analysis

We successfully sequenced 24 Bovine specimens excavated at four prehistoric sites in Iberia (Fig. 1A). Nine of these individuals were inferred to or suspected to represent aurochs based on morphology or chronology. Direct radiocarbon dates and contextual dating placed the individuals between the Mesolithic (oldest sample moo001, 8641-8450 cal BP) and the Roman Age (youngest sample moo010p, 2260-2150 CalBP). Based on the number of reads mapping to the X chromosome, 13 individuals were identified as female and eleven as male for the samples with sufficient amounts of reads for this analysis. Sequencing coverage of the nuclear cattle genome was low to medium, reaching up to 4.7x with a mean of 0.38x (Dataset S1). The sequence data for non-UDG-treated libraries showed damage patterns characteristic of ancient DNA (Supplementary Information, Figure S1).

Nine individuals were assigned to the mitochondrial P1 haplogroup, one to haplogroup T1, one to haplogroup T4, and 13 to haplogroup T3 (Dataset S1). P haplogroups are dominant and thought to be endemic to European aurochs (32). The prevalence of the T3 haplogroup in our samples is expected; this haplogroup is dominant among modern European B. taurus and is the most common haplogroup in ancient Western European domestic cattle. T3 was found in directly dated Neolithic samples from different sites providing direct evidence for the arrival of domestic cattle in northern Iberia during the Neolithic (Dataset S1). The specimen assigned to T1 (moo009a) is notable since this individual was previously used to argue for Bronze Age contact between Iberia and Africa, where the T1 haplogroup is thought to have originated (27, 33). T4 is usually considered to be restricted to Asian breeds with rare finds in Europe, restricted to the Balkans (23). The presence of T4 in Chalcolithic Iberia suggests that this haplogroup must have been distributed across Western Europe at low frequencies in prehistory. Furthermore, the fact that some specimens that were morphologically identified as aurochs carry domestic T haplogroups implies some level of interbreeding between the two groups.

As mitochondrial genomes only reflect the maternal line of ancestry, they are not informative about the exact extent of interbreeding in our dataset. To avoid being constrained by the variation in modern domestic breeds as with common approaches such as projected PCA (Supplementary Information, Figure S2), we performed non-metric multi-dimensional scaling (NMDS) ordinations on a matrix of pairwise outgroup f₃ statistics to explore the genomic ancestry of the sequenced individuals. For reference, we included early cattle genomes from Anatolia and the Balkans as well as aurochs excavated from Morocco (Th7), Armenia (Gyu2), Anatolia (Ch22) and Britain (CPC98) (5, 25), and calculated pairwise outgroup f₃ statistics. The NMDS ordination outcome (Fig. 1B) seems to represent a separation between domestic autosomal ancestry (to the left) and European aurochs ancestry (to the right). In contrast, aurochs from other regions (Th7 and Gyu2) seem genetically distinct. Many early domestic samples from Iberia fall close to early cattle from the Balkans and Anatolia as well as the Anatolian aurochs (Ch22). Notably, three of the Iberian samples in this cluster (moo004, moo007, moo015) were morphologically identified as aurochs. Eight Iberian samples fall to the right or top-right in the plot, together with the aurochs from Britain (CPC98); only one of these (moo019) carried a domestic T3 mitochondrial genome, the others are all P. Out of nine samples that were presumed aurochs based on their morphological features, only six would be considered aurochs based on this analysis. Moreover, 2 of the 8 samples that fall to the right of this plot were not allocated a species. This highlights a substantial overlap between measurements or criteria that are used to distinguish wild and domestic Bos based on morphometrics.

This analysis suggests that one can use the British aurochs (CPC98) as a reference for Western European aurochs as it seems similar to our three low-coverage Mesolithic Iberian samples. Indeed, we do not see any evidence against them forming a clade to the exclusion of all other samples used for the NMDS (Supplementary Table S2). Using CPC98 and Sub1, a Neolithic Anatolian individual as references, we can perform f₄ statistics to measure which Iberian individuals share more alleles with one or the other (Fig 1C). Despite the relatively low coverage of some samples, the f₄ statistics are highly correlated with the first axis of the NMDS (R²=0.77, p=1.75×10^-8) implying that they detect the same pattern. Non-overlapping confidence intervals also confirm that the high genetic differentiation between Western European aurochs and domestic cattle allows confident assignment even with low coverage data. The three Mesolithic individuals as well as an additional six up until the Late Neolithic/Chalcolithic share most of their alleles with aurochs. Three individuals from the Neolithic and Late Neolithic/Chalcolithic share most of their alleles with domestic Anatolian cattle while two individuals (moo013a and moo039) are intermediate and their confidence intervals are not overlapping with the others, suggesting that they are the result of hybridization. More recent samples from the Bell Beaker period onwards all appear to have similar amounts of allele sharing with mostly domestic ancestry but some level of aurochs introgression.

Quantifying the extent of introgression

While f₄ statistics measure allele sharing but do not directly quantify the amount of introgression in the different specimens, we employed the qpAdm framework (34, 35) to model each Iberian individual from British aurochs (CPC98) and/or Anatolian Neolithic cattle (Sub1) as sources (Fig. 2). While the point estimates of qpAdm are noisier than the f₄ statistics, qpAdm also rejects ancestry models that would not be compatible with the data for each individual. For all 24 samples, models that fit the data were found. From these, seven individuals could be modeled from an aurochs source alone while eight individuals could be modeled with only domestic ancestry. For these individuals, a single source is the most parsimonious explanation but it should be noted that for eleven of these fifteen, a model of mixed aurochs/domestic ancestry would fit the data as well. However, it is noteworthy that nine of these eleven have very limited amounts of sequence data, making it possible that more data would lead to a rejection of the single source model for some of them. Seven samples could not be modeled from a single source and only fitted as a mixture between aurochs and domestic cattle implying that they carried substantial proportions of both ancestries. The ancestry estimates in the 2-source model vary substantially between individuals ranging from ~25% (moo009a) to ~90% (moo014) aurochs ancestry. These results are largely consistent with the f₄ statistics, although some individuals falling to the extremes of the f₄ range also show signals of admixture. For some Mesolithic individuals, the two-source model is consistent with the data but the depth of coverage for these individuals was generally low and proportions of domestic ancestry are not significantly different from zero. During the Neolithic, two individuals are exclusively of Anatolian domestic ancestry without aurochs introgression while others show ancestry from both sources with predominantly wild ancestry in most cases. All Chalcolithic individuals are consistent with the two-source model with aurochs ancestry estimates ranging from 7 to 84%. The wide range of different ancestry estimates in the Neolithic and Chalcolithic samples suggests that interbreeding between wild aurochs and domestic cattle was frequent in prehistoric Iberia. However, at some point, the herding practices must have changed since modern Iberian breeds show approximately 20-25% aurochs ancestry (Supplementary Information, Table S2). In fact, we observe that after the Bell Beaker period the ancestry estimates change to predominantly domestic ancestry and aurochs proportions similar to those found in modern western European breeds (Fig. 3A & C).

Figure 2.

Figure 3.

A limitation of this analysis is the use of the British aurochs genome as a reference for the aurochs source population. There is a risk that this analysis will miss certain ancestries that are exclusive to Iberian aurochs while including some that are exclusive to the British aurochs, for example, contributions from glacial refugia other than Iberia. Our Mesolithic Iberian aurochs contained too little endogenous DNA to be used as a proxy aurochs reference and the three Neolithic and Chalcolithic samples estimated at >90% aurochs ancestry (including the 2.7x genome of moo014) already carry low levels of domestic ancestry (Supplementary Information, Figure S8).

An important question that remains unexamined is the exact process that led to the bidirectional introgression since this would provide insight into human management practices. The fact that some individuals with predominantly aurochs ancestry carry T haplogroups (e.g. moo019) and that some individuals with predominantly domestic ancestry carry P haplogroups (e.g. moo009x) implies that females contributed in both directions. To assess whether the admixture process was sex-biased, we compared the allele-sharing patterns on the X chromosome and autosomes. Since females carry two X chromosomes and males only have one, we can assume that an excess of a certain ancestry on the X chromosome indicates more females from that particular source population. For the Mesolithic samples, as expected, we do not see any significant difference in the level of allele sharing between X chromosome and autosomes (Fig. 3B). Similarly, there is no evidence of sex bias for most individuals with predominantly aurochs ancestry (negative autosomal f₄ values), i.e. individuals that are most likely wild rather than domestic or managed. In contrast, all other individuals with feasible and non-rejected two-source models of ancestry show an excess of allele sharing with Anatolian domestic cattle than with aurochs on the X chromosome, when compared to the autosomes. This difference is significant for eight of those twelve individuals including an individual carrying a P mtDNA haplogroup (moo009x). Even one of the individuals for which only the single source domestic qpAdm model is accepted shows an excess of Anatolian domestic ancestry on the X chromosome (moo004), suggesting that sex-biased processes already took place prior to the arrival of cattle to Iberia. This general observation is consistent with the relatively clear (but not perfect) discrimination between wild and domestic based on mitochondrial haplogroups where modern breeds mostly carry Near Eastern haplogroups (5, 21). In the absence of aurochs Y chromosomal data, however, it is difficult to assess sex-biased processes from uniparental data alone. The comparison of X chromosomes and autosomes should generally have more power to detect such processes as they are less sensitive to genetic drift due to their recombining nature (36). Overall, this illustrates that the contribution of wild ancestry into domestic cattle was mostly through aurochs bulls while no clear trend can be seen for admixture into wild populations.

Aurochs ancestry in modern breeds and the Spanish Lidia cattle breed

We estimated aurochs ancestry in a set of Western European cattle breeds (26) as we performed for the prehistoric samples. Previous studies have used D statistics for pairwise comparisons between breeds (25, 26, 37). Such D statistics, however, are sensitive to biases including gene flow from populations not included in the analysis (38). Furthermore, qpAdm provides the possibility to reject scenarios not fitting the data. Our point estimates for the aurochs ancestry range between 20% and 25% across all breeds (Fig. 3C) and do not show an increase in aurochs ancestry in Iberian breeds (37). This result differs from the previous studies which suggested geographic differences in western and central Europe and we believe this could be due to ancestry from other, non-European groups in some commercial breeds (Supplementary Information). Importantly, not all tested breeds did fit the simple two-source model Anatolian Neolithic domestic + British aurochs, likely representing low levels of contributions from other groups, e.g. indicine cattle (26). For instance, the two-source model was rejected for four of the ten commercial breeds in the dataset.

Cattle have played an important role in Iberian culture during the last centuries as they have been part of numerous traditional popular events including bullfighting. The Lidia breed, a heterogeneous group of Iberian cattle that is mainly bred for aggressive behavior, has commonly been used for such popular festivities (30). Even though Lidia cattle has only been actively bred for agonistic behavior for about 200 years, some people attribute their aggressiveness and appearance as an indication of high levels of aurochs ancestry (31). We additionally use medium coverage genomes of six Lidia individuals (39) to estimate their proportion of aurochs ancestry. Lidia cattle in the (26) data set had a point estimate of 21.6% (95% confidence interval: [16.1, 27.1]) aurochs ancestry and estimates in the individual genomes ranged from 15.8% [8.3, 23.4] to 26.1% [18.8, 33.5] (Supplementary Information, Figure S5) – all overlapping with the observed range for other western European breeds. Despite some variation between individuals, which might be attributable to noise due to low coverage sequencing data in the reference populations, we do not observe a systematically elevated level of aurochs ancestry compared to other modern breeds or ancient samples since the Bronze Age. While these results reject the idea that the specifics of Lidia cattle can be attributed to a substantially increased genome-wide aurochs ancestry, it does not rule out the possibility that the roots of their aggressiveness and appearance are indeed due to aurochs variants at key loci responsible for those traits. An in-depth investigation of such questions would require a larger dataset of aurochs genomes as well as a more comprehensive Lidia sampling due to their fragmentation in highly distant genetic lineages (30).

Stable Isotope Analysis

In addition to their ancestry, we studied the ecology of the bovids through stable isotope analysis of bone collagen. Lynch et al. (40) suggested that stable isotope data could be used to infer niche separation between the species in Britain, with domestic cattle in more open settings, while aurochs (about 1‰ more depleted in δ¹³C) were habitually in more forested areas, or wet ground. Most likely facilitated through human management of the domestic cattle, separating them from their wild counterparts. In contrast, Noe-Nygard et al. (41) failed to observe such an effect in samples from Denmark and northern Germany.

Considering our dataset (using genetically defined categories; B. primigenius and predominantly B. primigenius against B. taurus and predominantly B. taurus - to distinguish the two categories of Bos) and other data published on Iberian cattle (categorised on morphology/date) (42–51) the nitrogen isotope means are statistically different (see Dataset S1 and Supplementary Information). This is mostly due to some domestic cattle with δ¹⁵N values greater than 6.5‰. This could be explained by some taurine cattle having exclusive habitual access to high nitrogen isotope ratio resources, for example, human management such as corralling on manured ground, or feeding with manured crops, would produce this effect. Nevertheless, there is generally a large amount of overlap in the isotope values for the two groups.

Discussion

We generated and analyzed biomolecular data from B. primigenius and B. taurus spanning more than 9000 years in the same region. Cattle are important livestock in the Iberian Peninsula today, and our results illustrate the interaction between domestic cattle and their wild relatives in the past. The two groups show signs of frequent hybridization starting soon after the arrival of cattle to the peninsula, as evident in our oldest directly dated Neolithic individual (moo039, 7426-7280 CalBP) where signals of carrying both ancestries are clear as qpadm only considers the admixture model as fitting the data. Throughout the Neolithic, we observed large variations in the wild versus domestic ancestry per individual, but this pattern later stabilized (to 20-25% aurochs ancestry) from the Chalcolithic/Bronze Age onwards. This could reflect a transition from hunting and herding to predominantly herding and it is possible that systematic herd management led to the nearly constant levels of aurochs ancestry over the last 4000 years. This period also coincides with several other societal changes; including the Bell Beaker complex, and the introduction of human ancestry from the Pontic steppe and domestic horses into the Iberian Peninsula (9, 52–54). Around this time, humans also started processing a significantly higher amount of dairy products connected with the “secondary product revolution” (55, 56). Aurochs were probably present in Iberia until Roman times (28) leaving possibilities for interbreeding but we cannot exclude that various factors such as hunting or changing vegetation had led to a substantial decline in the wild aurochs population around the early Bronze Age. A previous study on cattle morphology from the site of El Portalón described a decrease in size from the Neolithic to the Chalcolithic and a further significant size decrease from the Chalcolithic to the Bronze Age (57) and associated this change in size to the aridification of the area at this time (58). Indeed, this climatic change could also be related to a reduction of the aurochs population contributing to the stabilization of the levels of ancestry in domestic cattle from the Bronze Age to the present. Nonetheless, our stable isotope results suggest that wild and domesticated groups often did not occupy substantially different niches on the Iberian Peninsula. Material excavated from Denmark suggested that aurochs changed their niches over time (41) demonstrating some flexibility depending on local vegetation and the possibility of aurochs adapting to changing environments.

The reduced level of aurochs ancestry on the X chromosome (compared to the autosomes) in admixed individuals suggests that it was mostly aurochs males who contributed wild ancestry to domestic herds. Consequently, the offspring of wild bulls and domestic cows could be born into and integrated within managed herds. It is unclear how much of this process was intentional but the possibility of a wild bull inseminating a domestic cow without becoming part of the herd suggests that some level of incidental interbreeding was possible. Modern breeders are still mostly exchanging bulls or sperm to improve their stock which manifests in a lower between-breed differentiation on the X chromosome (37).

The lack of correlation between genomic, stable isotope and morphological data highlights the difficulties of identifying and defining aurochs to the exclusion of domestic cattle. All of these data measure different aspects of an individual: their ancestry, ecology or appearance, respectively. While they can give some indication, none of them are a direct measurement of how these cattle were recognised by prehistoric humans, whether they were herded or hunted. It remains unclear whether our ancestry inferences had any correlation to how prehistoric herds were managed and how much intentional breeding is behind the observed pattern of hybridization. It is even possible that all hybrids identified in this study were part of domestic herds.

Even though wild aurochs populations went extinct, European aurochs ancestry survived into modern cattle with a relatively uniform distribution across western European breeds. Isolated Iberian Lidia, bred for their aggressiveness, appears to be no exception to this pattern. This rejects the notion that an overall increased proportion of aurochs ancestry causes the distinctiveness of certain breeds, but considering the functional relevance of archaic introgression into modern humans (59), it is possible that aurochs variants at functional loci may have a substantial influence on the characteristics of modern cattle breeds. Our low coverage sequencing data did not allow us to investigate this but future bioarchaeological studies combining different types of data will have the possibility to clarify the role of the extinct aurochs ancestry in modern domestic cattle.

Conclusions

Using a bioarchaeological approach we have demonstrated that since cattle arrived in Iberia there has been bidirectional introgression with the local aurochs population, but the wild to domestic direction is sex biased with mainly aurochs bulls contributing to domestic herds. Admixture proportions vary for the first few millennia but stabilize during the Early Bronze Age at approximately 20% of wild ancestry in domestic herds, a level that is still observed in modern Iberian breeds, including the more aggressive Lidia breed. This is likely to be the result of loose management of herds, becoming more controlled over time.

The amount of hybridisation observed in the ancient cattle makes it difficult to define what a domestic or wild Bos is, bringing into doubt the validity of such categorisations. Our interpretation is made more difficult by the overlap in morphological and metric data, creating further difficulties in species determination (especially in hybrids) and niche sharing as revealed by stable isotopes. To some extent, our interpretation is moot, as the salient matter is, how did prehistoric humans interact with cattle? What was their sense of wild and domestic and hybridisation? While we have recognised individual hybrids, to what extent these were part of domestic herds or intentionally bred and managed is uncertain.

Another source of uncertainty in our determinations is the limited knowledge about the genetic diversity in European aurochs, forcing us to use a British aurochs as a proxy for unadmixed Iberian aurochs. Further regional (and temporally longitudinal) aurochs genomes would aid future genomic studies defining the genetic variation in the European aurochs population.

Materials and Methods

Data generation

We attempted DNA extractions of 50 archaeological remains from which we successfully extracted DNA from 24 individuals identified as domestic cattle and aurochs excavated from four prehistoric sites in Iberia: El Portalón de Cueva Mayor (n=18), Artusia (n=1), Els Trocs (n=2) and Mendandia (n=3). Teeth and bones were UV irradiated (6 J/cm2 at 254 nm) and the first millimeter of bone/tooth surface abraded using a Dremel™ tool. DNA was extracted in a dedicated ancient DNA facility using a silica-based DNA extraction protocol (60). For each sample, 100-200mg of bone or tooth powder were incubated for 24 h at 37ºC, using the MinElute column Zymo extender assembly replaced by the High Pure Extender Assembly (Roche High Pure Viral Nucleic Acid Large Vol. Kit) and performed twice for each sample. DNA extracts were subjected to UDG treatment for the removal of deaminated cytosines and were further converted into blunt-end double stranded Illumina multiplex sequencing libraries (61). Between seven and fifthteen qPCR cycles were performed to amplify the DNA libraries using indexed primers (61). These were subsequently pooled at equimolar concentrations and shotgun sequenced on Illumina HiSeq and Novaseq sequencing platforms.

Radiocarbon dates

Eight Bone and teeth were directly radiocarbon dated (AMS) at Waikato University in New Zealand and two teeth at Beta Analytics in the United States. Radiocarbon dates were calibrated using the OXcal 4.4 program (62) and the IntCal20 calibration curve (63). Three samples from the site of Mendandia were conventionally radiocarbon dated at Groningen (Netherlands) radiocarbon laboratory and calibrated as above.

Stable Isotopes Analysis

Many of the samples analysed here were radiocarbon dated and stable isotope data (via IRMS) were generated in this process, to augment this data we also produced stable isotope data for some additional samples in this dataset, where they were available. The additional samples underwent bone collagen or tooth dentine collagen extraction at the Laboratorio de Evolución Humana (Universidad de Burgos) following the protocol of (64). In brief, this is a cold acid demineralization, followed by Milli Q water rinsing, gelatinization at pH3 (24 hrs at 70ºC), Ezee filtering and lyophilization. Collagen yields (as % mass of starting bone) were recorded. Stable isotope values (δ¹³C, δ¹⁵N) and %C, %N were measured in duplicate at the Universitat Autònoma de Barcelona, unless only one sample was successful in the analysis. Collagen samples (approx. 0.4 mg) were analysed using a Flash IRMS elemental analyser (EA) coupled to a Delta V Advantage isotope ratio mass spectrometer (IRMS), both from Thermo Scientific (Bremen, Germany) at the Institute of Environmental Science and Technology of the Universitat Autònoma de Barcelona (ICTA-UAB). International laboratory standard IAEA-600 was used, with measurements made relative to Vienna PeeDee Belemnite (V-PDB) for δ¹³C, and air N₂ (AIR) for δ¹⁵N. The average analytical error was <0.2‰ (1σ) as determined from the duplicate analyses of δ¹³C and δ¹⁵N. In house standards used was dog hair collected and homogenized for interlaboratory comparisons.

Data processing

HiSeq X10 reads have been trimmed and merged using AdapterRemoval (65) while adapters for NovaSeq 6000 reads have been trimmed with cutadapt (66) and merging was performed with FLASH (67) requiring a minimum overlap of 11bp. Single-end reads of at least 35bp length were then mapped to the cattle reference genomes UMD3.1 (68) and Btau5 (69) using bwa (70) with the non-default parameters: -l 16500, -n 0.01, and -o 2. Different sequencing runs per sample were merged with samtools (71) and consensus sequences were called for duplicate sequences with identical start and end coordinates (72). Finally, reads with more than 10% mismatches to the reference genome were removed. Biological sex was assigned to the samples mapped to the Btau_5 reference genome (as UMD3.1 does not contain a Y chromosome assembly) using the R_x method (73) modified for 29 autosomes.

For comparative purposes, we also processed published data from (5, 17, 18, 25) using the same bioinformatic pipeline. Furthermore, we downloaded sequence data for six Spanish Lidia cattle (39), a single modern water buffalo (Bubalus bubalis, Jaffrabadi-0845) (74) and a single zebu cattle individual (Sha_3b) (75) and processed them with our ancient DNA mapping pipeline. To obtain a pseudohaploid Yak (Bos grunniens) sequence, we followed the approach by (26) splitting the Yak reference genome (76) contigs into 100bp fragments and mapping them to the UMD3.1 reference genome.

Data Analysis

Mitochondrial consensus sequences were called using ANGSD (77) and the options-doFasta 2-doCounts 1-minQ 30-minMapQ 30. Mitochondrial haplogroups were then assigned to the whole mitogenome sequences using the Python script MitoToolPy (78).

For population genomic analysis, we used the panel of modern European breeds presented by (26) which were genotyped at ~770,000 SNPs. Prior to genotype calling, all ancient BAM files were modified such that Ts in the first 5 bases of each fragment and As at the last 5 base pairs of each fragment have a base quality of 2. To generate pseudohaploid representations of each individual, we randomly draw a single read with mapping and base quality of at least 30 at each SNP position. If the allele carried by the ancient individual was not one of the two known alleles, we removed the site from the panel.

To conduct an ordination of the nuclear data, sequences of 43 ancient Eurasian cattle and two aurochs were obtained from (25) and (5). Outgroup f₃ statistics were calculated for all pairs of our Iberian Bos samples, using a Yak (Bos grunniens) genome as an outgroup, and a distance matrix for all samples was calculated as 1-f₃. All f-statistics were calculated in R version 4.1.2 (79) package ‘admixtools2’ (80). The distance matrix was used to compute scores for non-metric multi-dimensional scaling (NMDS) ordinations using the metaMDS function in the ‘vegan’ R package and 10000 random starts (81).

We used admixtools2 (80) and qpAdm (34, 35) to model the ancestry proportions in the samples. CPC98 was used as a source for European aurochs ancestry while the domestic Anatolian Neolithic Sub1 was used as a source for domesticated cattle ancestry. As “right” populations, we used Bes2, Gyu2, B. indicus, Yak and Bison bonasus bonasus PLANTA. qpAdm was run with auto_only=FALSE and allsnps=TRUE. The same setup was used for modern western European breeds from (26) excluding breeds from Italy and the Balkan as non-taurine ancestry (26) in them would lead to a rejection of the models. We also used admixtools2 to reconstruct admixture graphs with the find_graphs() function. We first generated an f₂ matrix from the genotype data with maxmiss=0.15 and auto_only=FALSE and then ran find_graphs with stop_gen = 10000, stop_gen2 = 30 and plusminus_generations = 10 as suggested in (80). For between 0 and 5 migration events, we started the graph search with 1000 different seeds to select the model with the best likelihood for each number of migrations. Graph fits for different numbers of migrations were then compared using qpgraph_resample_multi() and compare_fits() using 100 bootstraps.

To obtain unbiased genotype calls for autosomes and the X chromosome for the analysis of sex bias, we called SNPs across all ancient and modern BAM files using ANGSD (77) and the parameters -GL 1 -doHaploCall 1 -doMajorMinor 1 -SNP_pval 1e-7 -doMaf 1 –minMaf 0.05 -minQ 30 -minMapQ 30 -remove_bads 1 -uniqueOnly 1 -rmTrans 1 -doCounts 1 -doGeno -4 -doPost 2 -doPlink 2 -minInd 3. f₄ statistics were then calculated using admixtools2 (80).

Acknowledgements

We are incredibly grateful to the excavation teams at the four archaeological sites. Sequencing was performed at The National Genomics Infrastructure (NGI) Stockholm. The computations and data handling were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) at Uppmax, partially funded by the Swedish Research Council through grant agreements no. 2022-06725 and no. 2018-05973. This project was supported by grants from the Royal Physiographic Society of Lund (Nilsson-Ehle Endowments) to T.G. and C.V., Vetenskapsrådet (2017-05267) to T.G., Ramón y Cajal (RYC2018-025223-I) to C.V and a Beatriz Galindo Fellowship (BGS220-461AA-69201) to C.S.