Mitochondrial DNA analysis shows a Near Eastern Neolithic origin for domestic cattle and no indication of domestication of European aurochs

(a) Samples

In this study, 106 Bos bones from many locations across Europe were assessed for survival of ancient DNA. Eighty-three bones had been previously differentiated as wild aurochs, rather than domestic cattle, on the basis of size or date, by the researchers who carried out the archaeozoological studies on the various sites. However, 9 bones were only tentatively labelled as aurochs, 11 either did not have determinations associated with them or were labelled as Bos (i.e. either B. p. primigenius or B. taurus), 2 were labelled as possibly bison or buffalo and 1 was differentiated as a domestic cattle bone. A further five samples, D740 (Bailey et al. 1996) and TP65, CHWF, NORF and CPC98 (Troy et al. 2001), were re-amplified as part of this study, and the aurochs sequence D812 (Bailey et al. 1996) was also included in the analyses. This made a total of 112 Bos specimens under consideration. Detailed information on sample provenance and analytical results is given in table S1 in the electronic supplementary material.

(b) Extraction and polymerase chain reaction amplification

The samples were analysed in three different laboratories: the Smurfit Institute of Genetics at Trinity College Dublin; the Institute of Anthropology at the University of Mainz; and the Henry Wellcome Ancient Biomolecules Centre at the University of Oxford. The analytical location for each sample is indicated in table S1 in the electronic supplementary material. Bone samples were prepared using previously described protocols (Yang et al. 1998; MacHugh et al. 2000; Burger et al. 2004; Shapiro et al. 2004). All primers were designed to be genus specific, if not species specific, and amplified fragments of the hyper-variable control region of the mitochondrial genome (see figure 4 in electronic supplementary material for the primer strategy employed by each of the three amplification laboratories). Polymerase chain reaction (PCR) set-up was conducted in laboratories dedicated solely to pre-amplification ancient work. PCR conditions and primer details were previously described: Bollongino et al. (2006) for Dublin and Mainz, and Shapiro et al. (2004) for Oxford. Cloning was as described in Bollongino et al. (2006). Second-round PCR was not undertaken on any samples that did not amplify in the first round. In Dublin, all non-amplifiable samples were tested for presence of inhibitors (Edwards et al. 2004).

(c) Statistical and phylogenetic analysis

The criteria for authenticating mitochondrial haplotypes were as previously described (Edwards et al. 2004; Bollongino et al. 2006). mtDNA sequences were aligned by eye. A reduced median network was constructed (figure 2) from the control region data using a median algorithm (Bandelt et al. 1995). Sequences were analysed using an HKY model of nucleotide substitution (Hasegawa et al. 1985), which was selected using the hierarchical likelihood ratio tests implemented by ModelTest v. 3.7 (Posada & Crandall 1998). The neighbour-joining method (Saitou & Nei 1987) was used to construct a dendrogram (table S4 in the electronic supplementary material) from genetic distances generated by the DNADIST program (Felsenstein 1989). Representative sequences of the five T (T, T1, T2, T3 and T4) and two Z (Z1 and Z2) cattle haplogroups, which predominate in taurine and zebu cattle, respectively, covering the 361 bp region, were located on GenBank, as well as representative outgroups of yak (Bos grunniens, AF083355), gaur (Bos gaurus, AF083371) and European bison (Bison bonasus, AF083356). Bootstraps were calculated from 1000 pseudoreplicates of the data.

Bayesian estimates of the mutation rate and the age of the most recent common ancestor (TMRCA) of the aurochs sequences were obtained using Beast v. 1.3; available from http://evolve.zoo.ox.ac.uk/beast/ (Drummond et al. 2006), from 37 P samples, covering a comparable 360 bp region. Sequences were analysed using an HKY model of nucleotide substitution (Hasegawa et al. 1985), which was selected using the hierarchical likelihood ratio tests implemented by ModelTest v. 3.7 (Posada & Crandall 1998). As the transition: transversion ratio could not be estimated from the data due to the lack of observed transversions, three different values (5, 50 and 500) were tested and yielded almost identical results. Two demographic models were tested for the coalescent prior: exponential growth and constant population size. The two models were compared using the average marginal posterior probabilities of the data, given the model. Mutation rate estimates were calibrated using the ages of the sequences, obtained by radiocarbon dating. The incorporation of these dates serves as sufficient calibration information for the estimation of the rate and divergence times (Drummond et al. 2002), and is more appropriate than using an external calibration point (Ho & Larson 2006). Posterior estimates of the mutation rate and age of the TMRCA were obtained by Markov Chain Monte Carlo analysis, with samples drawn every 500 steps over a total of 5 000 000 steps, following a discarded burn-in of 2 000 000 steps. Two separate runs were performed and the results combined. Adequate sampling and convergence to the stationary distribution were checked using Tracer v. 1.3 (Rambaut & Drummond 2004). Posterior estimates of parameters were all found to be distinctly unimodal (although with wide 95% highest posterior densities), and all parameters appeared to be identifiable, despite the relatively low information content in the sequences and the small age range of the sequences.

(a) Population history of the aurochs

This study is the most comprehensive, to date, of aurochs mtDNA sequence variation and has increased the number of available aurochs sequences to around fourfold. Previous studies have found the P haplogroup in aurochs from Britain and Iberia (Bailey et al. 1996; Troy et al. 2001; Anderung et al. 2005), but here we document its prevalence across the European continent (Hungary, Slovenia, Slovakia, Austria, Germany and France), dating from the Mesolithic to the Bronze Age. All P haplotypes fall into a star-like phylogeny, consistent with a past population expansion. Through Bayesian phylogenetic dating analysis, the coalescence of all P haplotypes was estimated to occur between 10 050 and 30 230 cal. years BP; a finding that is reconcilable with Late Glacial expansion from a refugial population.

(b) A novel aurochs haplotype

Bayesian phylogenetic analysis of the novel E haplotype, from the Early Neolithic German site of Eilsleben, estimates the age of the split between E and P lineages to be 52 700 years ago (95% HPD: 16 020–108 000 years). The bone analysed was an unusually robust distal humerus and was determined to be definitely a Bos sp. individual from morphological analysis. As this novel E sequence is an outlier to the P lineage, but groups within the B. taurus part of the tree (figure 3), it is plausible that its presence in Europe points to a greater variability of cattle prior to the Last Glacial Maximum. The location of E, between the P/T haplogroups and the Z haplogroup (figure 3), could also indicate a glacial refugial origin further to the East.

(c) Prior demonstration of pre-Neolithic T3 aurochs in Italy

The recent reporting of five T haplotypes in southern Italian aurochs (Beja-Pereira et al. 2006) is surprising and, on first viewing, contradictory to our data from northern and central Europe and the Near East. At the last glacial maximum, Italy was one of the three Southern European refugia for the Holocene fauna that subsequently repopulated Europe. Postglacial expansion of the majority of Italian refugial species was restricted by the barrier of the Alps to the North (Hewitt 1999). Therefore, with the presence of T haplotypes in the Italian aurochs population, we would expect them to remain there until the local aurochs population died out, rather than contributing to the northwards postglacial expansion of aurochs (characterized by P haplotypes in central and northern Europe), which instead is more likely to have been derived from the Iberian and Balkan refuges (Taberlet et al. 1998). A similar case is found in the native Italian wild boar, which are genetically distinct from the rest of Europe (Larson et al. 2005), supporting the notion that many Italian wild faunas are phylogenetically distinct. It is also possible that, at some point in the past, the Near Eastern and Italian aurochs belonged to the same ancestral population, perhaps connected via the Northern Adriatic basin, which remained submerged late during the Late Glacial transgression (van Andel 1989). Although we examined a Mesolithic aurochs sample from northern Italy, this bone yielded no reproducibly amplifiable DNA. Therefore, our one sequence from the region that could have connected the Near Eastern and Italian populations is the Mesolithic LJU3 from Slovenia, which exhibited a P haplotype. However, as it is situated at the northern extremity of the Balkan Peninsula, additional sampling from Italy and the Balkans is required to fully address this possibility. Given the above, and the results presented here, the unexpected haplotypes found in the Italian aurochs reported by Beja-Pereira et al. (2006) raise many questions and warrant additional verification.

(d) Post-domestication introgression between wild and domestic Bos populations

Eight European samples, classified on morphological grounds as aurochs by the corresponding archaeologists, yielded a T haplotype. These may have arisen due to: (i) morphological misclassification of domestic samples, as already postulated (Anderung et al. 2005; Bollongino et al. 2006), (ii) sporadic introgression from domestic cattle, or (iii) low-frequency haplotypes that survived in the predominantly P haplotype northern and eastern European aurochs population. Introgression (scenario (ii)) is certainly possible and plausible. Zooarchaeological diagnoses are inherently error prone, as they are based on size differences, and bone size depends upon many different genetic and environmental factors. It is conceivable that introgression might have been more frequent in the Early Neolithic, which would explain the predominantly Early Neolithic dates for the anomalous specimens. Scenario (iii) is also possible. However, as there are no pre-domestication T haplotypes apart from Italy and we did not find T haplotypes in our morphologically robust specimens, we favour hypotheses (i) and (ii).

Critical archaeozoological and archaeological examination of each of the eight discordant European samples leads to the conclusion that a domestic origin cannot be entirely excluded for any of them (see electronic supplementary material), again favouring the morphological misclassification explanation (scenario (i)). As the P haplogroup is not observed in extant cattle, births arising from hybridization are likely to have been rare, suggesting that domestic cattle were kept isolated from wild aurochs by humans. Even if P haplotypes were present in domestic cattle, but were subsequently lost by genetic drift, it is unlikely that such haplotypes were ever at high frequency in the B. taurus population. In sum, in all the samples that were morphologically identifiable, there is a clear-cut division between P haplotypes, found in morphologically robust specimens classified as aurochs, which were already present in Europe during the Pleistocene, and T haplotypes, found in less robust specimens that arrived in Europe during the Early Holocene. This shows that introgression was unlikely to have been a common event. However, as this study is based on mitochondrial data alone, we cannot rule out any introgression of male aurochs. Nevertheless, stable isotope data support the notion that these two subpopulations did not interact extensively. ¹³C/¹²C ratios yielded by bovine remains, from several Neolithic contexts in the Paris Basin (Balasse et al. 1997; Bocherens et al. 2005) and Denmark (Noe-Nygaard et al. 2005), clearly show that aurochs and domestic cattle, at least in these places, were feeding on distinct plant sources, the former in the forest and the latter in more open environments.

(e) Near Eastern origin of taurine haplotype

Modern genetic and archaeological evidence suggests that the domestication centre of European cattle was in the Near East. Here, we present, for the first time, ancient biomolecular data that establish the presence of the T family of haplogroups in an archaic Near Eastern population. The sample Syria17 (table 1; table S1 in the electronic supplementary material) is a Bos specimen from Dja'de el Mughara in the Middle Euphrates Valley. This site was occupied mainly during the Early Pre-Pottery Neolithic (PPNB), ca 10 650–10 250 cal. BP (Coqueugniot 2000). Archaeozoological Bos material from this site has long been considered as wild, but recent morphometric studies have revealed a very small decrease of size and overall a slight and significant decrease of the sexual dimorphism that convincingly suggest an incipient domestication (Helmer et al. 2005). This result is in concordance with the observation that domestic cattle were transferred to Cyprus from the mainland as early as ca 10 250–10 150 cal. BP (Vigne et al. 2000, 2003). The Dja'de sample, which lay undisturbed 3 m beneath the surface of a Tell, yielded a T3 haplotype, a result that was replicated in both Dublin and Mainz. Although previous work has shown the occurrence of the T haplotype in early Near Eastern cattle (Edwards et al. 2004), and we report here a Bronze Age aurochs specimen from Maral Tappeh, Iran, with a T haplotype, this Syrian sample is the earliest dated DNA evidence for the occurrence of this haplogroup in the Near East. A taxonomic discrimination between B. taurus and B. p. primigenius forms that are chronologically and geographically close to the centre of domestication is arbitrary, and thus the bone sample Syria17 may be an example of a very early taurine domesticate, or alternatively an Early Holocene aurochs. In either case, this is the first confirmation of the presence of a common B. taurus haplotype close to the suggested centre of domestication, at the very beginning of the Neolithic.