Abstract
Cattle have been a valuable economic resource and cultural icon since prehistory. From the initial expansion of domestic cattle into Europe during the Neolithic period, taurine cattle (Bos taurus) and their wild ancestor, the aurochs (B. primigenius), had overlapping ranges leading to ample opportunities for intentional and unintentional hybridization. We performed a bioarchaeological analysis of 24 Bos remains from Iberia dating from the Mesolithic to the Roman period. The archaeogenomic dataset allows us to investigate the extent of domestic-wild hybridization over time, providing insight into the species’ behavior and human management by aligning changes with cultural and genomic transitions in the archaeological record. Our results show frequent hybridization during the Neolithic and Chalcolithic, likely reflecting a mix of hunting and herding or relatively unmanaged herds, with mostly male aurochs and female domestic cattle involved in hybridization. This is supported by isotopic evidence of ecological niche sharing, with only a few domestic cattle possibly being managed. The proportion of aurochs ancestry in domestic cattle remains relatively constant from about 4000 years ago, probably due to herd management and selection against hybrids, coinciding with other cultural transitions. The constant level of wild ancestry (~20%) continues into modern western European breeds including the Spanish Lidia breed which is bred for its aggressiveness and fighting ability, but does not display elevated levels of aurochs ancestry. This study takes a genomic glance at the impact of human actions and wild introgression in the establishment of cattle as one of the most important domestic species today.
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.
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 f3 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 f3 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 f4 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 f4 statistics are highly correlated with the first axis of the NMDS (R2=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 f4 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 f4 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 f4 statistics, although some individuals falling to the extremes of the f4 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).
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 f4 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 δ13C) 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 δ15N 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 (δ13C, δ15N) 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 δ13C, and air N2 (AIR) for δ15N. The average analytical error was <0.2‰ (1σ) as determined from the duplicate analyses of δ13C and δ15N. 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 Rx 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 f3 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-f3. 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 f2 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. f4 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.
References
- 1.Ancient pigs reveal a near-complete genomic turnover following their introduction to EuropePNAS 116:17231–17238
- 2.Animal domestication in the era of ancient genomicsNat Rev Genet :1–12
- 3.Ancient European dog genomes reveal continuity since the Early NeolithicNat Commun 8
- 4.Paleogenomics of Animal Domesticationin Paleogenomics Springer https://doi.org/10.1007/13836_2018_55
- 5.Ancient cattle genomics, origins, and rapid turnover in the Fertile CrescentScience 365:173–176
- 6.Ancient goat genomes reveal mosaic domestication in the Fertile CrescentScience 361:85–88
- 7.Herded and hunted goat genomes from the dawn of domestication in the Zagros MountainsProceedings of the National Academy of Sciences 118
- 8.Tracking Five Millennia of Horse Management with Extensive Ancient Genome Time SeriesCell 177:1419–1435
- 9.The origins and spread of domestic horses from the Western Eurasian steppesNature 598:634–640
- 10.Archaeogenetic analysis of Neolithic sheep from Anatolia suggests a complex demographic history since domesticationCommun Biol 4:1–11
- 11.Origins and genetic legacy of prehistoric dogsScience 370:557–564
- 12.Ancient Sheep Genomes reveal four Millennia of North European Short-Tailed Sheep in the Baltic Sea regionbioRxiv
- 13.Iron age genomic data from Althiburos – Tunisia renew the debate on the origins of African taurine cattleiScience 26
- 14.Early Animal Husbandry in the Northern LevantPaléorient 25:27–48
- 15.Identifying early domestic cattle from Pre-Pottery Neolithic sites on the Middle Euphrates using sexual dimorphismin The First Steps of Animal Domestication: New Archaeozoological Approaches :86–96
- 16.Oxbow BooksEvolution and domestication of the Bovini speciesAnimal Genetics 51:637–657
- 17.Genome data on the extinct Bison schoetensacki establish it as a sister species of the extant European bison (Bison bonasus)BMC Evol Biol 17
- 18.Complex Admixture Preceded and Followed the Extinction of Wisent in the WildMolecular Biology and Evolution 34:598–612
- 19.Genome-wide local ancestry and direct evidence for mitonuclear coadaptation in African hybrid cattle populations (Bos taurus/indicus)
- 20.The origin of European cattle: Evidence from modern and ancient DNAProceedings of the National Academy of Sciences 103:8113–8118
- 21.Mitochondrial genomes of extinct aurochs survive in domestic cattleCurrent Biology 18:R157–R158
- 22.Oxbow BooksIncorporation of aurochs into a cattle herd in Neolithic Europe: single event or breeding?Sci Rep 4
- 23.Large-scale mitogenome sequencing reveals consecutive expansions of domestic taurine cattle and supports sporadic aurochs introgressionEvolutionary Applications 15:663–678
- 24.Ancient DNA analysis of Scandinavian medieval drinking horns and the horn of the last aurochs bullJournal of Archaeological Science 99:47–54
- 25.Genome sequencing of the extinct Eurasian wild aurochs, Bos primigenius, illuminates the phylogeography and evolution of cattleGenome biology 16
- 26.Genetic origin, admixture and population history of aurochs (Bos primigenius) and primitive European cattleHeredity
- 27.Prehistoric contacts over the Straits of Gibraltar indicated by genetic analysis of Iberian Bronze Age cattleProceedings of the National Academy of Sciences of the United States of America 102:8431–8435
- 28.El problema de la domesticación de bovinos en el País Vasco y resto de la Región CantábricaCongresos de Estudios Vascos :123–128
- 29.Re-evaluating the Neolithic: The Impact and the Consolidation of Farming Practices in the Cantabrian Region (Northern Spain)J World Prehist 29:79–116
- 30.Genetic variation within the Lidia bovine breedAnimal Genetics 39:439–445
- 31.Intensive human contact correlates with smaller brains: differential brain size reduction in cattle typesProceedings of the Royal Society B: Biological Sciences 288
- 32.The Draft Genome of Extinct European Aurochs and its Implications for De-ExtinctionOpen Quaternary 2
- 33.Detecting the T1 cattle haplogroup in the Iberian Peninsula from Neolithic to medieval times: new clues to continuous cattle migration through timeJournal of Archaeological Science 59:110–117
- 34.Massive migration from the steppe was a source for Indo-European languages in EuropeNature https://doi.org/10.1038/nature14317
- 35.Assessing the performance of qpAdm: a statistical tool for studying population admixtureGenetics 217
- 36.Accelerated genetic drift on chromosome X during the human dispersal out of AfricaNat Genet 41:66–70
- 37.Consequences of breed formation on patterns of genomic diversity and differentiation: the case of highly diverse peripheral Iberian cattleBMC Genomics 20
- 38.Bias in estimators of archaic admixtureTheoretical population biology 100:63–78
- 39.A novel missense variant in endothelin-2 (EDN2) causes a growth and respiratory lethal syndrome in bovineAnimal Genetics 53:583–591
- 40.Where the wild things are: aurochs and cattle in EnglandAntiquity 82:1025–1039
- 41.Diet of aurochs and early cattle in southern Scandinavia: evidence from 15N and 13C stable isotopesJournal of Archaeological Science 32:855–871
- 42.Living different lives: Early social differentiation identified through linking mortuary and isotopic variability in Late Neolithic/ Early Chalcolithic north-central SpainPLOS ONE 12
- 43.Investigating palaeodietary and social differences between two differentiated sectors of a Neolithic community, La Bòbila Madurell-Can Gambús (north-east Iberian Peninsula)Journal of Archaeological Science: Reports 3:160–170
- 44.Palaeodiets of Humans and Fauna at the Spanish Mesolithic Site of El ColladoCurrent Anthropology 47:549–556
- 45.Investigating prehistoric diet and lifeways of early farmers in central northern Spain (3000–1500 CAL BC) using stable isotope techniquesArchaeol Anthropol Sci 11:3979–3994
- 46.Feeding Management Strategies among the Early Neolithic Pigs in the NE of the Iberian PeninsulaInternational Journal of Osteoarchaeology 27:839–852
- 47.Estudio de la dieta en la población neolítica de Costamar: Resultados preiliminares de análisis de isótopos estables de Carbono y NitrógenoTorre la Sal (Ribera de Cabanes, Castellón): evolución del paisaje antrópico desde la Prehistoria hasta el Medioevo :411–418
- 48.Aproximación a la dieta de la población calcolítica de La Vital a través del análisis de isótopos estables del carbono y del nitrógeno sobre restos óseosLa Vital (Gandia, Valencia). Vida y Muerte En La Desembocadura Del Serpis Durante El III y El I Milenio a.C :139–143
- 49.Isotope evidence for the use of marine resources in the Eastern Iberian MesolithicJournal of Archaeological Science 42:231–240
- 50.Late Neolithic-Chalcolithic socio-economical dynamics in Northern Iberia. A multi-isotope study on diet and provenance from Santimamiñe and Pico Ramos archaeological sites (Basque Country, Spain)Quaternary International 481:14–27
- 51.Reconstruction of human subsistence and husbandry strategies from the Iberian Early Neolithic: A stable isotope approachAmerican Journal of Physical Anthropology 167:257–271
- 52.The Beaker phenomenon and the genomic transformation of northwest EuropeNature 555:190–196
- 53.Four millennia of Iberian biomolecular prehistory illustrate the impact of prehistoric migrations at the far end of EurasiaPNAS
- 54.The genomic history of the Iberian Peninsula over the past 8000 yearsScience 363:1230–1234
- 55.Neolithic to Bronze Age economy and animal management revealed using analyses lipid residues of pottery vessels and faunal remains at El Portalón de Cueva Mayor (Sierra de Atapuerca, Spain)Journal of Archaeological Science 131
- 56.Dairying, diseases and the evolution of lactase persistence in EuropeNature 608:336–345
- 57.Metrical analysis of bovine bone remains from the Neolithic to the Bronze Age at the El Portalón site (Atapuerca, Burgos) in the Iberian contextQuaternary International 566:211–223
- 58.Upper Pleistocene and Holocene palaeoenvironmental records in Cueva Mayor karst (Atapuerca, Spain) from different proxies: speleothem crystal fabrics, palynology and archaeologyInternational Journal of Speleology 43
- 59.The contribution of Neanderthal introgression to modern human traitsCurrent Biology 32:R970–R983
- 60.Complete mitochondrial genome sequence of a Middle Pleistocene cave bear reconstructed from ultrashort DNA fragmentsProc. Natl. Acad. Sci. U.S.A 110:15758–15763
- 61.Illumina Sequencing Library Preparation for Highly Multiplexed Target Capture and SequencingCold Spring Harbor Protocols 2010
- 62.Bayesian analysis of radiocarbon datesRadiocarbon 51:337–360
- 63.The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0–55 cal kBP)Radiocarbon 62:725–757
- 64.Stable Isotope Evidence for Similarities in the Types of Marine Foods Used by Late Mesolithic Humans at Sites Along the Atlantic Coast of EuropeJournal of Archaeological Science 26:717–722
- 65.AdapterRemoval v2: rapid adapter trimming, identification, and read mergingBMC Research Notes 9
- 66.Cutadapt removes adapter sequences from high-throughput sequencing readsEMBnet.journal 17:10–12
- 67.FLASH: fast length adjustment of short reads to improve genome assembliesBioinformatics 27:2957–2963
- 68.A whole-genome assembly of the domestic cow, Bos taurusGenome Biology 10
- 69., Btau_5.0.1 -Genome -Assembly -NCBI
- 70.Fast and accurate short read alignment with Burrows-Wheeler transformBioinformatics 25:1754–1760
- 71.The Sequence Alignment/Map format and SAMtoolsBioinformatics 25:2078–2079
- 72.Analysis of high-throughput ancient DNA sequencing datain Methods in Molecular Biology (Clifton, N.J Springer :197–228
- 73.A Molecular Approach to the Sexing of the Triple Burial at the Upper Paleolithic Site of Dolní VěstonicePLOS ONE 11
- 74.Whole genome analysis of water buffalo and global cattle breeds highlights convergent signatures of domesticationNat Commun 11
- 75.Whole-genome resequencing reveals world-wide ancestry and adaptive introgression events of domesticated cattle in East AsiaNat Commun 9
- 76.The Yak genome database: an integrative database for studying yak biology and high-altitude adaptionBMC Genomics 13
- 77.ANGSD: Analysis of Next Generation Sequencing DataBMC Bioinformatics 15
- 78.DomeTree: a canonical toolkit for mitochondrial DNA analyses in domesticated animalsMolecular Ecology Resources 15:1238–1242
- 79.R: A Language and Environment for Statistical Computing(R Foundation for Statistical Computing
- 80.On the limits of fitting complex models of population history to f-statisticseLife 12
- 81.vegan: Community Ecology Package
Article and author information
Author information
Version history
- Preprint posted:
- Sent for peer review:
- Reviewed Preprint version 1:
Copyright
© 2024, Günther et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics
- views
- 332
- downloads
- 22
- citations
- 0
Views, downloads and citations are aggregated across all versions of this paper published by eLife.