Email updates

Keep up to date with the latest news and content from BMC Evolutionary Biology and BioMed Central.

Open Access Research article

Population dynamic of the extinct European aurochs: genetic evidence of a north-south differentiation pattern and no evidence of post-glacial expansion

Stefano Mona1, Giulio Catalano2, Martina Lari2, Greger Larson3, Paolo Boscato4, Antonella Casoli5, Luca Sineo6, Carolina Di Patti7, Elena Pecchioli8, David Caramelli2 and Giorgio Bertorelle1*

  • * Corresponding author: Giorgio Bertorelle ggb@unife.it

  • † Equal contributors

Author Affiliations

1 Dipartimento di Biologia ed Evoluzione, Università di Ferrara, Ferrara, Italy

2 Dipartimento di Biologia Evoluzionistica, Università di Firenze, Firenze, Italy

3 Department of Archaeology, Durham University, Durham, UK

4 Dipartimento di Scienze Ambientali, Università di Siena, Siena, Italy

5 Dipartimento di Chimica Generale e Inorganica, Chimica Analitica, Chimica Fisica, Università di Parma, Parma, Italy

6 Dipartimento di Biologia Animale, Università di Palermo, Palermo, Italy

7 Museo di Geologia e Paleontologia G Gemellaro, Università di Palermo, Palermo, Italy

8 IASMA Research and Innovation Centre, Fondazione Edmund Mach, Trento, Italy

For all author emails, please log on.

BMC Evolutionary Biology 2010, 10:83  doi:10.1186/1471-2148-10-83

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1471-2148/10/83


Received:15 October 2009
Accepted:26 March 2010
Published:26 March 2010

© 2010 Mona et al; licensee BioMed Central Ltd.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Background

The aurochs (Bos primigenius) was a large bovine that ranged over almost the entirety of the Eurasian continent and North Africa. It is the wild ancestor of the modern cattle (Bos taurus), and went extinct in 1627 probably as a consequence of human hunting and the progressive reduction of its habitat. To investigate in detail the genetic history of this species and to compare the population dynamics in different European areas, we analysed Bos primigenius remains from various sites across Italy.

Results

Fourteen samples provided ancient DNA fragments from the mitochondrial hypervariable region. Our data, jointly analysed with previously published sequences, support the view that Italian aurochsen were genetically similar to modern bovine breeds, but very different from northern/central European aurochsen. Bayesian analyses and coalescent simulations indicate that the genetic variation pattern in both Italian and northern/central European aurochsen is compatible with demographic stability after the last glaciation. We provide evidence that signatures of population expansion can erroneously arise in stable aurochsen populations when the different ages of the samples are not taken into account.

Conclusions

Distinct groups of aurochsen probably inhabited Italy and northern/central Europe after the last glaciation, respectively. On the contrary, Italian and Fertile Crescent aurochsen likely shared several mtDNA sequences, now common in modern breeds. We argue that a certain level of genetic homogeneity characterized aurochs populations in Southern Europe and the Middle East, and also that post-glacial recolonization of northern and central Europe advanced, without major demographic expansions, from eastern, and not southern, refugia.

Background

Ancient DNA offers the unique opportunity to broaden the time depth of population genetic analyses. Moreover, it is the only tool capable of deciphering the demographic histories of extinct species [1-4]. The aurochs, Bos primigenius, was one of two wild bovine species in Europe (Bison bonasus was the other). Probably as a consequence of human hunting and habitat reduction [5], the aurochs went extinct in 1627, when the last individual died in Poland. This species most likely originated in India between 1.5 and 2 MYA ago and later occupied many parts of Asia, North Africa and Europe. The first documented remains in Europe dates back to 275,000 years ago [5]. Its demographic history is largely unknown and, in particular, it is unclear to what extent it was shaped by Pleistocene climate changes. Indeed, the genetic structure and population history of many European species were largely influenced by the alternating of glacial and inter-glacial periods [3,6-8].

When modern samples are analysed, refugial areas generally harbour greater genetic diversity and multiple phylogenetic clusters with deeply rooted branches that indicate their antiquity beyond the most recent Ice Age; conversely, if we exclude the contact zones where distinct migratory routes mix [9], non-refugial areas display reduced diversity and show only a subset of the phylogenetic clusters present in the refugia [7]. Based mainly on genetic evidence of this type, four models have been proposed to describe the post-glacial recolonization routes in Europe [7,10,11]. The models differ in the relative contribution of the refugia populations to the expansion process, and the genetic data analysed for several species seemed always to fit the predictions of at least one of these models. Recently, however, a growing number of studies reports unexpected patterns of diversity [e.g., [12-14], and even within refugia, the genetic structure might have been much more complex than previously acknowledged [15]. Recent theoretical studies (reviewed in [16]) also suggest that a range expansion from a single refugium can lead to a geographically structured population, which can be easily misinterpreted as the result of an expansion from multiple distinct refugia. Therefore, the geographic distribution of genetic clusters might not always be informative for identifying recolonization routes.

The analysis of DNA from ancient samples represents a key step to overcome these problems. It opens a genetic-variation window directly onto the past, useful to investigate ancient population dynamics and to compare the demographic history of refugial and non-refugial populations [3,17]. In particular, ancient DNA is a fundamental tool available to investigate aurochs dynamics. In fact, the aurochs surviving descendents, cattle, have a complex history of selection and displacement which heavily interferes with any attempt to reconstruct past events that occurred in the wild.

DNA data from ancient aurochs specimens can also be informative for the reconstruction of the domestication process (e.g., [18,19]). Recent archaeological evidence suggests that taurine cattle (Bos taurus) were initially domesticated in the upper Euphrates Valley around10,000-11,000 years B.P. [20] This single origin hypothesis implies that present day European breeds, which belong to the taurine type, all descended from Fertile Crescent ancestors. The phylogenetic dichotomy between northern and central European Bos primigenius (P and E) and Bos taurus (T) mtDNA clades was interpreted as a clear evidence supporting the single origin hypothesis [19,21,22]. But five Italian pre-domesticate aurochsen with T haplotypes [18], and a P haplotype recently found in one Korean cow generically classified as "beef-cattle" [23] raise doubts about the primigenius-taurus mtDNA distinctiveness, and, as a consequence, about the single origin hypothesis. As suggested also by the large genetic diversity in the Italian, compared to other European breeds [18,23-26], the domestication process in Europe, or at least in Italy, requires a careful reconsideration.

For this study we generated partial mtDNA control region sequences from 14 ancient aurochsen samples collected from a geographically diverse range of locations across Italy. Samples from Southern refugia are particularly valuable since DNA preservation is uncommon in warmer climates [27]. Indeed, several aurochs sequences have been analyzed in northern and central Europe [19,21,22] but only few sequences have been described so far from Southern areas [18,22,28]. The main goals of our study were: i) to investigate the genetic variability and the demographic history of Italian aurochsen, especially in comparison with the northern and central European groups (i.e., comparing a refugia vs. a non-refugia population); ii) to study whether and how the analysis of aurochs data collected at different time intervals (heterochronous sampling) may produce erroneous conclusions when simple statistics which assume simultaneous (isochronous) sampling are used; and iii) to test previous insights suggesting a role for the Palaeolithic aurochs in the domestication process in Italy [18].

Results

aDNA laboratory results

The values of amino acid racemization obtained in all our samples suggested good biochemical preservation (< 0.10 for aspartic acid). This result appears compatible with DNA survival [29,30], although some concern have been raised on the utility of this approach [31,32]. The alanine D/L and glutamine D/L are lower compared to the aspartic acid D/L, which is expected since the racemization of alanine and glutamine is slower than that of aspartic acid. This difference suggests that modern (contaminant) aminoacids are absent [30]. Sporadic contamination is unlikely when the number of molecules that PCR will use as template (or target DNA) is greater than 1,000 [30]. By using real-time PCR, we quantified the number of DNA molecules in the extracts from 20 samples, selected to be representative of different Palaeolithic and Late Mesolithic sites (Table S1, see Additional file 1); we found well preserved DNA in 18 of the 20 samples, with numbers of DNA molecules varying from 1,000 to more than 5,000 (Table S1, see Additional file 1). After positive extraction and amplification of 14 samples, and cloning of PCR products for all them, we obtained the entire control region sequence for one sample; 13 other samples yielded partial sequences from at least one of the three fragments in which HVR-I had been subdivided (Table S2, see Additional file 2). Independent replication of extraction, amplification of the highly variable first fragment (134 bp), and cloning was performed at CEA (Trento) ancient DNA laboratory for seven bone samples (Table S2, see Additional file 2). The replicated sequences matched those obtained in the Florence laboratory (see Table S2, see Additional file 2). Only sporadic G-to-A change and C-to-T change were observed in some clones, likely due to post-mortem cytosine deamination [33]. The Taq misincorporation value for all 202 cloned sequences was estimated to be 1.21 substitutions, and 87% of the clones possessed the consensus sequence. These results suggest that the DNA templates were not significantly damaged.

Additional file 1. Table S1. Summary of the major characteristics of the samples, including radiocarbon dates.

Format: XLS Size: 24KB Download file

This file can be viewed with: Microsoft Excel ViewerOpen Data

Additional file 2. Table S2. The clone sequences.

Format: DOC Size: 104KB Download file

This file can be viewed with: Microsoft Word ViewerOpen Data

Analysis of haplotype diversity

On the basis of the diagnostic sites identified in previous studies [19,22,25], 12 sequences were assigned to the T haplogroups and two sequences belong to the P haplogroup (Table S3, see Additional file 3). In other words, the vast majority of the Italian aurochsen analysed in this study showed the same type of mtDNA sequence observed previously in 5 Italian aurochsen [18] and in all the extant European cattle breeds analysed so far, though two individuals share the same haplogroup with the aurochsen from northern and central Europe (and with a single cattle sequenced in Korea generically classified as "beef cattle" [23]). All the DNA sequences are available from GenBank under the accession numbers from GU434117 to GU434130.

Additional file 3. Table S3. The consensus sequence for each sample.

Format: DOC Size: 34KB Download file

This file can be viewed with: Microsoft Word ViewerOpen Data

The dataset was then reduced to ten sequences using an overlapping fragment of 120 bp from position 16042 to 16161. These sequences were jointly analysed with 55 previously published aurochs sequences from Italy and northern/central Europe [18,22], and 800 modern cattle representing 53 breeds from Europe, Turkey, Middle East and North Africa (the list of breeds and the sample sizes are available in Table S4, see Additional file 4). In this data set, the group of Italian aurochsen (hereafter, Au-It) and the group of northern and central European aurochsen (hereafter, Au-NCE) were composed of 15 and 50 individuals respectively. Even though they are derived from different geographical areas, the Au-NCE group included specimens with highly similar sequences. In fact, when Au-NCE was split into three geographic groups (Britain, east central Europe, and west central Europe), no significant genetic divergence was found (AMOVA: Φst = 0.02, P = 0.16). The Φst genetic distances between the two groups of aurochsen and the modern breeds were graphically represented by a MDS analysis (Figure 1).

Additional file 4. Table S4. The geographic origin, name, samples size, and reference for all the taurine breeds considered in this study.

Format: XLS Size: 28KB Download file

This file can be viewed with: Microsoft Excel ViewerOpen Data

thumbnailFigure 1. Two-dimensional MDS plot based on Φst distances between aurochs population and 53 taurine cattle breeds; stress value = 0.053. Au-NCE = Northern and central European aurochs; Au-It = Italian aurochs; MOL (Morocco Olmez), LYB (Lybian Local), EGB (Egyptian Baladi) and EGY (Egypt) are all cattle breeds from Africa. The complete list of the cattle breeds used for this figure is reported in Table S4 (see Additional file 4).

As observed in a previous study based only on five Italian aurochsen [18], the Italian Bos primigenius fell into the cluster that included all modern breeds in Europe. North African breeds, due to the prevalence of T1 haplotypes in this region, appear related, but clearly distinct from European breeds. Conversely, the Bos primigenius specimens excavated outside Italy, joined here in the Au-NCE group, and all characterized by P sequences, are distantly related to all the other modern or ancient samples. The Φst value between the two ancient groups (Au-It vs. Au-NCE) was very high (0.727, P < 0.0001), which is consistent with their haplogroup frequencies. However, this value needs to be viewed with caution since the two groups both included samples with different ages; nevertheless, these results clearly suggest the existence of a strong north-south geographic structure in this species before domestication, with very limited gene flow across the Alps.

The phylogeographic structure of the aurochs haplotypes (see Figure 2, where Au-NCE is subdivided as before into geographic areas) was characterized by three major clades: 1) T sequences only found in Italy, 2) P sequences separated among them by no more than two mutations and distributed throughout the European regions (including Italy), and 3) a highly divergent E group found only in one German Neolithic sample [22]. These results, which are in agreement with previous analyses based on longer mtDNA fragments [22], indicate that the 120 bp fragment used in our data-set contains enough information to discriminate between these clades.

thumbnailFigure 2. Median-joining network of the 120 bp aurochs haplotypes. The branch length is proportional to the number of substitutions; the node diameter is proportional to the haplotype frequency. The names of the major haplogroups are shown.

Considering the P group of sequences as representative of the genetic variation in northern and central Europe (thus excluding the divergent E sequence), a star-like pattern suggesting a demographic increase in these area was identified in previous studies [21,22]. The phylogenetic structure in Italy is clearly different: the highly divergent T and P clades, typical of modern breeds and northern and central European aurochs populations, respectively, are both observed in pre-Neolithic samples.

Choice of the best fitting demographic model

The demographic history of Bos primigenius in Italy and northern/central Europe, and the possible differences in these two areas, was further analysed using BEAST under two models: constant population size and exponential change. This analysis was performed separately for the Au-It and the Au-NCE groups.

The Au-NCE data set was reduced to the 37 sequences published and analysed by [22] so that the sequences retained their full 360 bp fragment size, and we could also compare our results with those in [22]. Using the natural logarithmic scale of Jeffreys to interpret the Bayes Factor (as reported in [34]), we found "positive" evidence in favour of the exponential change model (BF = 3.7). The median value of the posterior distribution of the growth rate, g, was small and negative (around -2*10-5), i.e. compatible with a limited population contraction. However, the 95% high posterior density, though clearly left-shifted, included zero (Table 1) suggesting that demographic stability could not be excluded. In other words, in contrast with the analysis reported in [22], the data for the northern and central European aurochs appeared inconsistent with a population expansion: a demographic model with either population size stability or a very limited decline seems more compatible with the genetic variation pattern.

Table 1. Demographic parameters estimated by BEAST under the a priori model of exponential growth.

When the same analysis was applied to the Italian aurochs (Au-It data set), the Bayes Factor was much smaller (1.6 in favour of the exponential change model, "not worth more than a bare mention" following the Jeffreys natural logarithmic scale), as well as the estimated growth rate (see Table 1). A demographic model with constant population size can therefore explain the Italian aurochs data set, and Au-It showed a long-term effective population size twice as large as Au-NCE (Table 1; population sizes are estimated assuming a generation time of 7 years [35]).

Simulating the effects of heterochronous sampling

The BEAST results contrast with both the very similar analysis of the Au-NCE data set by [22], and with the common assumption that the star-like haplotype network in Au-NCE indicates a population expansion. We will discuss the first point in the next section, after we first report the results of our investigation of the second point using a simulation approach. Specifically, we wanted to understand how readily a star-like haplotype network could be interpreted as evidence for a population expansion when the samples are heterochronous. The empirical distributions of Tajima's D [36] and Fu's Fs [37] were used as indicators for the network shape. The expectation of both these statistics, assuming an isochronous sampling, is zero under a constant population size model, but becomes negative under a population expansion model which also produces a star-like genealogy. The effect of pooling DNA data with different ages (heterochronous sampling) on these expectations seems to depend strictly on details of the temporal sampling scheme[38]. The results we obtained when the samples ages mirror the aurochs data set reported in [22] clearly indicate that the distribution of Tajima's D and Fu's Fs are skewed towards negative values not only when a demographic expansion is simulated, but also when the population size is constant (see Table 2 and Figure 3).

Table 2. Simulation results, based on 1,000 replicates for each scenario (sample size is fixed to 37 individuals, as in the real Au-NCE sample).

thumbnailFigure 3. Box-plot of Tajima's D (a) and Fu's Fs (b) obtained from 1,000 simulations under a constant population size model with various Ne and assuming either an heterochronous (HET) or an isochronous (IS) sampling scheme (see text for details).

Under constant population size, therefore, an artificial signal of population expansion is likely to arise when the individuals are not sampled from the same time period. The relationship between the fraction of negative and significant D and Fs values and the different parameters of the simulations is complex and requires more systematic investigations (as in a recent publication [38])). It seems, however, that in all cases where isochronous sampling produces small fractions of significant statistics (i.e., for g = 0 or g = 0.001), the fraction of significant statistics produced by a heterochronous sampling is several times larger.

We have so far demonstrated that an expansion signal can be generated by heterochronous sampling when the temporal information is not taken into account, and that coalescent analyses of both Au-It and Au-NCE were suggestive of a constant, rather than an expanding, population. To further explore this issue, we performed a Bayesian skyline analysis for both aurochs populations using the ages of the sequences as temporal information. We then compared the reconstructed skylines with those estimated assuming that the ancient sequences were sampled simultaneously (isochronous sampling). For clarity of comparison, the isochronous skylines were depicted assuming that the sequences were all sampled at 2,000 and 7,000 years B.P. for Au-NCE and Au-It respectively (i.e., the ages of the most recent sample). The temporal dynamic of the median values is shown in Figure 4. Neither Au-NCE nor Au-It showed evidence of population size changes. However, a strong signal of population expansion emerged when the same groups were considered to be composed by isochronous samples (Figure 4). The 95% highest posterior density (not shown in Figure 4) for all the four analyses was wide, indicating that the results were not very informative regarding the absolute estimate of Ne. However, the 95% limits followed the same trend of the median, suggesting that the overall temporal pattern of Ne was robust.

thumbnailFigure 4. Bayesian skyline plot of the effective population sizes through time for Au-It (Au-It HET) and Au-NCE (Au-NCE HET). The same groups were analyzed assuming an isochronous sampling at 7,000 and 2,000 years B.P. respectively (Au-It IS and Au-NCE IS; see text for details); only the median values are shown.

Discussion

Over the last decade, the study of ancient DNA has proven to be particularly helpful in deciphering the evolutionary history of populations [39-42], often revealing unexpected demographic dynamics [2-4,14,43-45]. In this study we investigated the mitochondrial genetic variability of the extinct aurochs.

The distributions of mtDNA haplotypes in Italy and northern/central European areas are very different. Sequences belonging to the P clade, typical of northern and central Europe and not observed in five Italian specimens analysed by [18] were found only in two of the additional 14 sequences here produced. Overall, the vast majority of the Italian aurochsen analysed so far (17 out of 19) have mtDNA sequences belonging to the T clade, which is the phylogenetic group of sequences observed in virtually all European cattle breeds.

The genetic similarity between the Italian aurochs and the European cattle breeds may have important implication for the study and the understanding of the domestication process in this species. Interpretations of the genetic data over the last 15 years, including two recent publications [26,28], have supported both single and multiple origin narratives. It seems, however, that these two hypotheses are based on extreme, oversimplified models developed when wild and domestic status could be determined by mtDNA sequences. Both the process and the data are more complex. European wild and domestic forms, which coexisted at least since the beginning of the Linearbandkeramik culture, around 7,500 years ago [46], shared many genotypic as well as phenotypic traits, and secondary domestication or introgression events possibly occurred in some areas but not in others.

The remarkable difference between the patterns of aurochs genetic variation in Italy and in northern/central Europe is one of the two main results of our study, and confirms the conclusions of a previous study based on only five Italian aurochsen [18]. The other is the lack of demographic size fluctuations in this species in both geographical regions. This latter result conflicts with the conclusions of two previous studies [22,47] which supported a demographic, possibly postglacial, expansion in northern and central Europe (Italy was excluded from these analyses). We discuss in the next three paragraphs the likely reasons for this discrepancy.

Star-like shapes in sequence networks are commonly taken as evidence for population expansion, and this was also inferred for the aurochs considering the structure of the P clade [21,22]. The star-like expectation [48], however, assumes that all samples have the same age (isochronous sampling), which is not the case for the aurochs data. As recently shown by [38], the analysis of samples with different ages (heterochronous sampling) can introduce a systematic bias in some classical population genetics analyses, but the direction of the bias cannot be easily predicted. Using a simulation approach, we found that a signal of population expansion, as estimated by neutrality tests strictly dependent on the network shape, can be generated in a population with a constant size when individuals are sampled at the same time points as in the aurochs data set. In other words, the shape of the network in mtDNA sequences observed in northern and central Europe cannot be taken per se as evidence for population expansion.

More specifically, [22] concluded for the northern and central European aurochs that an exponential growth model gave a better fit than a constant size model, as indicated by marginal posterior probabilities. Using the same data set, and comparing the models not on the basis of the marginal posterior probabilities but on the basis of the marginal likelihoods using the more statistically sound Bayes factor, we reached different conclusions. There was no support for an exponential growth model over a constant size model. A similar result was obtained analysing the Italian aurochsen, where the Bayes factor in favour of a demographic change model, as compared to a constant population size model, was even smaller. The genetic data, therefore, do not support a demographic expansion and were more compatible with a model of constant population size. These results were confirmed when the Bayesian skyline plot was properly generated taking into account the temporal structure of the data.

The demographic expansion model for the aurochs was also suggested in a recent study on several bovine species [47]. However, in this analysis the estimated effective population size for the aurochs had a large, uninformative 95% credible interval of 3*104 - 9*109 and a median value 2*109 (after two independent MCMC with 108 steps using the input file provided in Ho et al. 2008). In addition, MCMC showed poor mixing and high autocorrelation for the population size and the growth rate. This multi-species approach is possibly uninformative regarding the aurochs demography, and requires further investigations (Simon Ho, personal communication).

The Pleistocene glacial and inter-glacial periods heavily impacted genetic variation in many species [6,8,10], and the steppe-tundra environment in central and northern Europe during the last glaciation was probably not suitable for the aurochs [5]. Does the genetic pattern we observed for the aurochs tell us something about the population dynamics in this species during the last severe climatic change in Europe?

First of all, genetic variation is higher in Italy compared to northern/central Europe (see also the effective population size estimates in Table 1), suggesting that a typical Southern refugium, Italy, might have acted as a variation reservoir during the glaciations. We can predict a certain degree of genetic homogeneity between Italy and other South-Eastern Mediterranean areas, considering that i) the Italian aurochs, and the European and Middle-Eastern cattle breeds, which likely descend mainly from aurochsen domesticated in the Fertile Crescent, share very similar mtDNA haplotypes; and ii) three early Neolithic specimens from Syria and Turkey, with a possible attribution to Bos primigenius [22,28], had T sequences, i.e. sequences belonging to the mtDNA clade typical in the Italian aurochs. It is therefore possible that aurochsen living in several Mediterranean regions, sharing similar ecological conditions and belonging to a species with a large migratory potential, exchanged enough genetic material to maintain a certain degree of homogeneity. Second, the low frequency of P sequences in Italy, which are the only sequences found in northern areas, indicates that these areas were not re-colonized by Italian or South-Eastern migrants. The high genetic divergence between Italian and northern/central European aurochsen is compatible with a long-term separation and small levels, if any, of gene flow. Third, the Bayesian analysis suggests that post-glacial warming and likely range expansions were not accompanied in the aurochs by a demographic expansion, at least an expansion intense enough to leave a clear signature on the pattern of genetic variation.

If these conclusions are correct, we should finally address the following question: where is the geographical origin of the aurochsen living after the last glaciation in the northern and central Europe, a region which was probably re-colonized in warmer times by a group of migrants without a significant demographic increase?

Armed with the current data we cannot rule out the hypothesis that the Iberian peninsula was a genetically divergent refugium in Southern Europe which contributed to the post-glacial colonization. We note however that the occurrence of continuous gene flow across southern Europe (including Iberia) was already proposed for the brown bear [14,49], supporting the idea that the model of southern refugia in Europe inhabited by genetically differentiated populations can be violated in some species. At the same time, the role of eastern, and not south-eastern, refugia, i.e. areas from the Carpathian to the Caspian regions, and further East in Asia, might have been important for the aurochs, as well as for other species [7,10,50,51]. We should not forget that some regions in Eastern Asia hosted temperate forests with large areas of grassland during the last glaciation [52], where relict populations of aurochsen could have survived before a post-glacial westwards migration. In this regard, it is intriguing to consider that a domestic individual with uncertain origin (possibly belonging to the Holstein breed) but with a P sequence was recently reported from a laboratory in Korea [23]. More recently, a survey of cytochrome b sequences in modern cattle found a P haplotype in another Holstein cattle as well as in an individual of the indigenous Korean Youngju Yellow breed [53]. Additional ancient DNA data from different areas will possibly clarify this issue.

A final warning

Unfortunately, our conclusions regarding the northern and central European aurochsen, as well as those reported in all the other studies on the same topic (e.g.,[18,19,28]), rely on the assumption that T sequences were not present before domestication in that area. This is a plausible view, but the evidence supporting it is limited, since i) the number of pre-Neolithic specimens excavated and successfully typed from that area is small; ii) several individuals mainly from Germany, morphologically classified as "Bos sp", or "Bos taurus/Bos primigenius?", or "Bos primigenius?" or "Bos taurus?" had T haplotypes ([28], Supplementary Materials). It seems therefore that the idea that T haplotypes were absent in northern and central European aurochsen is not really based on the accumulation of data in favour of this hypothesis. Rather, this idea may have more to do with the results of early studies based on few aurochs sequences, and with some level of uncertainty in the morphological attribution of ancient skeletal remains. In the words of the authors, a more conservative conclusion should be "...there is no conclusive evidence of aurochs with a T haplotype in Europe" ([22], Supplementary Materials).

Of course, this conclusion should not be applied to Europe in general, since we have seen that a large fraction of Italian pre-Neolithic aurochsen (where the morphological attribution does not represent a problem) have T sequences. In conclusion, it is possible that the geographical distribution of T haplotypes before the arrival in Europe of domestic breeds was wider than commonly assumed. If this is true, the pattern of north-south population divergence in the aurochs could be less extreme than is described here, and the hypothesis that modern cattle breeds have ancestors not only in the Middle East and possibly in Italy[18], but also in other European regions [54], should be subjected to a robust investigation.

Conclusion

The remarkable genetic difference between Italian aurochsen and the group of northern/central European aurochsen, and the similarity between the former and modern taurine breeds, is fully supported by the ancient mtDNA sequences typed in this study. Using extensive simulation and Bayesian based coalescent approaches we found no evidence of a demographic expansion in either Italian or northern/central European aurochsen groups. In fact, our observations lead us to conclude that the previously reported signatures may be artefactual. We interpret our results as evidence that the aurochs population sizes in southern Europe were not affected by climatic change in the last 30,000 years, and that central and northern areas in Europe were recolonized after the last glaciation by eastern immigrants without a major demographic expansion.

Methods

Samples collection and laboratory analyses

A total of 56 Bos primigenius samples from different well-studied sites across Italy were analyzed. To avoid the possible attributions of fragmented remains from domestic animals to Bos primigenius, we considered only samples dating before the Neolithic spread of domestication. Radiocarbon dating of associated remains (47 samples) or stratigraphic contests (for nine samples) were used to select the samples to analyse [55-57].

Details of each sample are provided in Table S1 (see Additional file 1). DNA was extracted for all samples in the Florence laboratory, and, for a subset of samples (see below), also in the Trento laboratory. Both laboratories are dedicated exclusively to ancient DNA work.

DNA extractions and PCR set up were carried out in rooms physically separated from those in which PCR cyclings and post-PCR analysis were conducted. Multiple measures were undertaken to exclude contamination and potential artifactual DNA changes (including post-mortem damage).

Disposable masks, gloves, and laboratory coats were worn during the experiments and were changed frequently. Pipettors were UV-irradiated in between use. DNA extractions and PCR reactions included negative controls. In each set of extractions or purifications, we included a negative control, represented by all the reagents except the bone powder, and these negative controls, together with blanks (all amplification reagents minus DNA), were regularly analysed in every PCR experiment to control for presence of exogenous DNA. All steps of the analysis were replicated at least twice. Seven bone specimens were independently extracted, amplified, cloned and sequenced in the ancient DNA laboratory of the Centre of Alpine Ecology (Edmund Mach Fundation, Trento). To test for preservation of other macromolecules as an indirect evidence for DNA survival [30,58] we estimated the degree of amino acid racemization in each sample. Approximately 5 mg of bone powder were used in this analysis [30]. The amount of target DNA was quantified in a sub-sample of specimens by Real Time (RT) PCR.

Ancient DNA analysis in Florence

DNA extraction

All DNA-preparation and extraction methods followed strictly specific ancient DNA requirements [27,29,59]. To prevent contamination from prior handling, the outer layer of bones was removed with a rotary tool, and the fragments were briefly soaked in 10% bleach. All samples were then irradiated (1 hour under UV light) and powdered. DNA was extracted by means of a silica-based protocol (modified from [60]). Two independent extracts were obtained from each remain. Multiple negative controls were included in each extraction.

Quantification of DNA Molecules

Real-time PCR amplification was performed using Brilliant® SYBR® Green QPCR Master Mix (Stratagene) in MX3000P (Stratagene), using 0.5 μM of appropriate primers (reverse primer located at L16030 and forward primer located at H16165). Thermal cycling conditions were 95°C for 10 min, 40 cycles at 95°C for 30 s, 53°C for 1 min and 72°C for 30 s, followed by SYBR® Green dissociation curve steep. Ten-fold serial dilutions of the purified and quantified standard were included in the experiment to create the standard curve in order to know the number of initial DNA molecules in the samples.

Amplification of mtDNA

Two μl of DNA extracted from the bone were amplified with this profile: 94°C for 10 min (Taq polymerase activation), followed by 50 cycles of PCR (denaturation, 94°C for 45 sec, annealing, 52°C for 1 min and extension, 72°C for 1 min) and final step at 72°C for 10 min. The 50 μl reaction mix contained 2 U of AmpliTaq Gold (Applied Biosystems), 200 μM of each dNTP and 1 μM of each primer. The HVR-I (Hyper variable region one of D-loop Mitochondrial DNA) was subdivided in three overlapping fragments of 134, 116, and 104 bais pairs, respectively. Primers used were: L16030 5'-ACATTAAATTATATGCCCCATGC-3'; H16165 5'-TTCACGCGGCATGGTA-3'; L16159 5'-'TTCCTTACCATTAGATCACGAGC-3; H16276 5' GATGAGATGGCCCTGAAGAA-3'; L16267 5'-CAATGAATTTTACCAGGCAT-3'; H00034 5'-CCAAATGTGACAGCACAG-3'. Each extract was amplified at least twice. Since overlapping primers were used throughout the PCR amplifications, it is highly unlikely that we amplified a nuclear insertion rather than the organellar mtDNA.

Cloning and Sequencing

PCR products were cloned using the TOPO TA Cloning Kit (Invitrogen) according to the manufacturer's instructions. Screening of white recombinant colonies was accomplished by PCR, transferring the colonies into a 30 μl reaction mix (67 mM Tris HCl [pH 8.8], 2 mM MgCl2, 1 μM of each primer, 0.125 mM of each dNTP, 0.75 units of Taq Polymerase) containing M13 forward and reverse universal primers. After 5 min at 92°C, 30 cycles of PCR (30 sec at 90°C, 1 min at 50°C, 1 min at 72°C) were carried out and clones with insert of the expected size were identified by agarose gel electrophoresis. After purification of these PCR products with Microcon PCR devices (Amicon), a volume of 1,5 μl was cycle-sequenced following the BigDye Terminator kit (Applied Biosystems) using the supplier's instructions. The sequence was determined using an Applied BioSystems 3100 DNA sequencer. The clones sequences for each sample are reported in Table S2 (see Additional file 2).

Independent replications in Trento

The first DNA fragment, which includes most of the variable sites, was replicated for seven bone samples in the CEA laboratory in Trento. The same ancient DNA protocols followed in Florence were adopted in Trento, with the exception of the Strata Clone TM PCR Cloning Kit.

Statistical analyses

A median-joining network [61] among aurochs haplotypes, defined by the segregating sites of the partial control region stretch (120 bp), was constructed using the software NETWORK 4.5.0.0 (available at http://fluxus-engineering.com webcite) with default settings. The level of divergence between populations/breeds was estimated using the analysis of molecular variance [62], as implemented in ARLEQUIN 3.11 [63]. Pairwise Φst values were simply computed using the number of differences as a molecular distance between alleles, with standard errors based on 1,000 bootstrap replicates. The distance matrix was plotted in two dimensions by means of multidimensional scaling (MDS) using the ALSCAL algorithm implemented in the SPSS 12.0 software package.

The demography of the Italian aurochs population was investigated with the software BEAST 1.4.7 [64] using two different coalescence prior models (constant size and exponential change) and the HKY model of nucleotide substitution, as selected by MODELTEST 3.7 according to the hierarchical likelihood ratio test method [65]. The transition/transversion ratio could not be estimated from the data (because of the presence of a single transversion in one sample only) and was set as in [22]; final results are based on Ti/Tv = 50. The demographic models were run three times for 10,000,000 MCMC iterations with a 20% burnin and a thinning interval of 1,000, and were compared by means of the Bayes Factor. Bayes Factor was computed as twice the difference between the log of the marginal likelihoods, which were approximated using the harmonic mean as suggested by [34]. Convergence was checked by examining the generation plot visualized with TRACER 1.3 [66] and the final estimates were based on the pooled runs (after removing the burnin). TMRCA of the whole sample as well as of the T and P haplogroups were recorded within both demographic models, after checking that T and P were monophyletic with posterior probabilities of 1.00 (hence confirming that the partial control region contained enough information to distinguish between the two haplogroups). The mutation rate was estimated from the data, using the ages of the sequences to calibrate the clock. This procedure is considered more appropriate than using an external calibration point [67].

To check if the demography of the Italian aurochs differed from that estimated for the northern and central European aurochs, we analysed the data in [22] with the same settings as above. We also explored the population history of the European aurochs to determine whether heterochronous sampling could skew the distribution of summary statistics used to detect changes in the effective population size. In this case, we performed a grid of simulation with various combinations of Ne and growth rates under a temporal sampling scheme analogous to that in [22]. Generation time was set to seven years [35] and mutation rate to 7.6 × 10-7 per site per year (mean of the posterior distribution of the mutation rate estimated by [22]). The software SERIALSIMCOAL [68] was used for the simulations. For each scenario we performed 1,000 replicates, and Tajima's D and Fu' Fs were computed with ARLEQUIN 3.11 [63] as indices of population size changes. The empirical distribution of these statistics was reconstructed under various models, and compared to the null distributions expected assuming an isochronous sampling scheme under the same coalescent models. For each simulated scenario we computed the percentage of significant tests using the procedure implemented in ARLEQUIN (i.e., assuming an isochronous sampling scheme).

Finally, a Bayesian skyline plot [69] was constructed for both the Italian and the northern/central European aurochs population. We ran the analysis twice for each group. Length of the MCMC was set to 20,000,000 iterations with a 10% burnin and a thinning interval of 1,000. The gene genealogy was divided into three internodes groups and effective population size function was fitted with a piecewise constant function of population size change. The reconstructed skyline was also compared for both groups with the skyline estimated assuming an isochronous sampling scheme.

Abbreviations

MYA: Millions of years; MCMC: Monte Carlo Markov Chain; TMRCA: Time to the Most Recent Common Ancestor; BF: Bayes Factor.

Authors' contributions

GC, ML, EP and GL performed aDNA laboratory analyses. AC performed the racemization analysis. PB, LS and CDP provided samples and radiocarbon/stratigraphic information. SM and GB performed the statistical analyses. GB and DC conceived the project. GB, DC, SM and GC wrote the paper. All authors read and approved the manuscript.

Acknowledgements

We are grateful to Roscoe Stanyon and Guido Barbujani for critical comments on the manuscript, Matteo Girardi for technical assistance, and Simon Ho for discussions regarding aurochs demography. GL was supported by an EMBO fellowship. BioTethed and the University of Florence supported this work (ex 60%) and the Ph.D program of GC. We thank the University of Ferrara (FAR funds) and the Italian Ministry for Research and Education (PRIN funds) for supporting this work and a fellowship to SM.

References

  1. Gilbert MT, Drautz DI, Lesk AM, Ho SY, Qi J, Ratan A, Hsu CH, Sher A, Dalen L, Gotherstrom A, et al.: Intraspecific phylogenetic analysis of Siberian woolly mammoths using complete mitochondrial genomes.

    Proc Natl Acad Sci USA 2008, 105(24):8327-8332. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  2. Hofreiter M, Munzel S, Conard NJ, Pollack J, Slatkin M, Weiss G, Paabo S: Sudden replacement of cave bear mitochondrial DNA in the late Pleistocene.

    Curr Biol 2007, 17(4):R122-123. PubMed Abstract | Publisher Full Text OpenURL

  3. Hofreiter M, Serre D, Rohland N, Rabeder G, Nagel D, Conard N, Munzel S, Paabo S: Lack of phylogeography in European mammals before the last glaciation.

    Proc Natl Acad Sci USA 2004, 101(35):12963-12968. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  4. Larson G, Albarella U, Dobney K, Rowley-Conwy P, Schibler J, Tresset A, Vigne JD, Edwards CJ, Schlumbaum A, Dinu A, et al.: Ancient DNA, pig domestication, and the spread of the Neolithic into Europe.

    Proc Natl Acad Sci USA 2007, 104(39):15276-15281. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  5. Van Vuure C: Retracing the Aurochs: History, Morphology and Ecology of an extinct Wild Ox.

    Pensoft 2005. OpenURL

  6. Hewitt G: The genetic legacy of the Quaternary ice ages.

    Nature 2000, 405(6789):907-913. PubMed Abstract | Publisher Full Text OpenURL

  7. Hewitt GM: Genetic consequences of climatic oscillations in the Quaternary.

    Philos Trans R Soc Lond B Biol Sci 2004, 359(1442):183-195.

    discussion 195.

    PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  8. Taberlet P, Fumagalli L, Wust-Saucy AG, Cosson JF: Comparative phylogeography and postglacial colonization routes in Europe.

    Mol Ecol 1998, 7(4):453-464. PubMed Abstract | Publisher Full Text OpenURL

  9. Petit RJ, Aguinagalde I, de Beaulieu JL, Bittkau C, Brewer S, Cheddadi R, Ennos R, Fineschi S, Grivet D, Lascoux M, et al.: Glacial refugia: hotspots but not melting pots of genetic diversity.

    Science (New York, NY) 2003, 300(5625):1563-1565. OpenURL

  10. Hewitt G: Post-glacia re-colonization of European biota.

    Biological Journal of the Linnean Society 1999, 68:87-112. Publisher Full Text OpenURL

  11. Schmitt T: Molecular biogeography of Europe: Pleistocene cycles and postglacial trends.

    Front Zool 2007, 4:11. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  12. Hansson B, Hasselquist D, Tarka M, Zehtindjiev P, Bensch S: Postglacial colonisation patterns and the role of isolation and expansion in driving diversification in a passerine bird.

    PLoS ONE 2008, 3(7):e2794. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  13. Tougard C, Renvoise E, Petitjean A, Quere JP: New insight into the colonization processes of common voles: inferences from molecular and fossil evidence.

    PLoS ONE 2008, 3(10):e3532. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  14. Valdiosera CE, Garcia-Garitagoitia JL, Garcia N, Doadrio I, Thomas MG, Hanni C, Arsuaga JL, Barnes I, Hofreiter M, Orlando L, et al.: Surprising migration and population size dynamics in ancient Iberian brown bears (Ursus arctos).

    Proc Natl Acad Sci USA 2008, 105(13):5123-5128. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  15. Canestrelli D, Cimmaruta R, Nascetti G: Population genetic structure and diversity of the Apennine endemic stream frog, Rana italica--insights on the Pleistocene evolutionary history of the Italian peninsular biota.

    Mol Ecol 2008, 17(17):3856-3872. PubMed Abstract | Publisher Full Text OpenURL

  16. Excoffier L, Ray N: Surfing during population expansions promotes genetic revolutions and structuration.

    Trends Ecol Evol 2008, 23(7):347-351. PubMed Abstract | Publisher Full Text OpenURL

  17. Rowe KC, Heske EJ, Brown PW, Paige KN: Surviving the ice: Northern refugia and postglacial colonization.

    Proc Natl Acad Sci USA 2004, 101(28):10355-10359. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  18. Beja-Pereira A, Caramelli D, Lalueza-Fox C, Vernesi C, Ferrand N, Casoli A, Goyache F, Royo LJ, Conti S, Lari M, et al.: The origin of European cattle: evidence from modern and ancient DNA.

    Proc Natl Acad Sci USA 2006, 103(21):8113-8118. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  19. Troy CS, MacHugh DE, Bailey JF, Magee DA, Loftus RT, Cunningham P, Chamberlain AT, Sykes BC, Bradley DG: Genetic evidence for Near-Eastern origins of European cattle.

    Nature 2001, 410(6832):1088-1091. PubMed Abstract | Publisher Full Text OpenURL

  20. Helmer D, Gourichon L, Monchot H, Peters J, Saña Segui M: Identifying early domestic cattle from prepottery Neolithic sites on the middle Euphrates using sexual dimorphism. In The First Steps of Animal Domestication. Edited by Vigne JD, Peters J, D H. Oxford: Oxbow Books; 2005:86-95. OpenURL

  21. Bailey JF, Richards MB, Macaulay VA, Colson IB, James IT, Bradley DG, Hedges RE, Sykes BC: Ancient DNA suggests a recent expansion of European cattle from a diverse wild progenitor species.

    Proc Biol Sci 1996, 263(1376):1467-1473. PubMed Abstract | Publisher Full Text OpenURL

  22. Edwards CJ, Bollongino R, Scheu A, Chamberlain A, Tresset A, Vigne JD, Baird JF, Larson G, Ho SY, Heupink TH, et al.: Mitochondrial DNA analysis shows a Near Eastern Neolithic origin for domestic cattle and no indication of domestication of European aurochs.

    Proc Biol Sci 2007, 274(1616):1377-1385. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  23. Achilli A, Olivieri A, Pellecchia M, Uboldi C, Colli L, Al-Zahery N, Accetturo M, Pala M, Kashani BH, Perego UA, et al.: Mitochondrial genomes of extinct aurochs survive in domestic cattle.

    Curr Biol 2008, 18(4):R157-158. PubMed Abstract | Publisher Full Text OpenURL

  24. Negrini R, Nijman IJ, Milanesi E, Moazami-Goudarzi K, Williams JL, Erhardt G, Dunner S, Rodellar C, Valentini A, Bradley DG, et al.: Differentiation of European cattle by AFLP fingerprinting.

    Anim Genet 2007, 38(1):60-66. PubMed Abstract | Publisher Full Text OpenURL

  25. Pellecchia M, Negrini R, Colli L, Patrini M, Milanesi E, Achilli A, Bertorelle G, Cavalli-Sforza LL, Piazza A, Torroni A, et al.: The mystery of Etruscan origins: novel clues from Bos taurus mitochondrial DNA.

    Proc Biol Sci 2007, 274(1614):1175-1179. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  26. Achilli A, Bonfiglio S, Olivieri A, Malusa A, Pala M, Kashani BH, Perego UA, Ajmone-Marsan P, Liotta L, Semino O, et al.: The multifaceted origin of taurine cattle reflected by the mitochondrial genome.

    PLoS One 2009, 4(6):e5753. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  27. Willerslev E, Cooper A: Ancient DNA.

    Proc Biol Sci 2005, 272(1558):3-16. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  28. Bollongino R, Elsner J, Vigne JD, Burger J: Y-SNPs do not indicate hybridisation between European aurochs and domestic cattle.

    PLoS ONE 2008, 3(10):e3418. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  29. Paabo S, Poinar H, Serre D, Jaenicke-Despres V, Hebler J, Rohland N, Kuch M, Krause J, Vigilant L, Hofreiter M: Genetic analyses from ancient DNA.

    Annu Rev Genet 2004, 38:645-679. PubMed Abstract | Publisher Full Text OpenURL

  30. Poinar HN, Hoss M, Bada JL, Paabo S: Amino acid racemization and the preservation of ancient DNA.

    Science 1996, 272(5263):864-866. PubMed Abstract | Publisher Full Text OpenURL

  31. Collins MJ, Penkman KE, Rohland N, Shapiro B, Dobberstein RC, Ritz-Timme S, Hofreiter M: Is amino acid racemization a useful tool for screening for ancient DNA in bone?

    Proc Biol Sci 2009, 276:2971-2977. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  32. Schwarz C, Debruyne R, Kuch M, McNally E, Schwarcz H, Aubrey AD, Bada J, Poinar H: New insights from old bones: DNA preservation and degradation in permafrost preserved mammoth remains.

    Nucleic acids research 2009, 37(10):3215-3229. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  33. Hofreiter M, Jaenicke V, Serre D, Haeseler Av A, Paabo S: DNA sequences from multiple amplifications reveal artifacts induced by cytosine deamination in ancient DNA.

    Nucleic Acids Res 2001, 29(23):4793-4799. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  34. Kass R, Raftery A: Bayes factors.

    Journal of the American Statistical Association 1995, 90:773-795. Publisher Full Text OpenURL

  35. Gautier M, Faraut T, Moazami-Goudarzi K, Navratil V, Foglio M, Grohs C, Boland A, Garnier JG, Boichard D, Lathrop GM, et al.: Genetic and haplotypic structure in 14 European and African cattle breeds.

    Genetics 2007, 177(2):1059-1070. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  36. Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism.

    Genetics 1989, 123(3):585-595. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  37. Fu YX: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection.

    Genetics 1997, 147(2):915-925. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  38. Depaulis F, Orlando L, Hanni C: Using classical population genetics tools with heterochroneous data: time matters!

    PLoS One 2009, 4(5):e5541. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  39. Barnes I, Matheus P, Shapiro B, Jensen D, Cooper A: Dynamics of Pleistocene population extinctions in Beringian brown bears.

    Science 2002, 295(5563):2267-2270. PubMed Abstract | Publisher Full Text OpenURL

  40. Belle EM, Ramakrishnan U, Mountain JL, Barbujani G: Serial coalescent simulations suggest a weak genealogical relationship between Etruscans and modern Tuscans.

    Proc Natl Acad Sci USA 2006, 103(21):8012-8017. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  41. Dalen L, Nystrom V, Valdiosera C, Germonpre M, Sablin M, Turner E, Angerbjorn A, Arsuaga JL, Gotherstrom A: Ancient DNA reveals lack of postglacial habitat tracking in the arctic fox.

    Proc Natl Acad Sci USA 2007, 104(16):6726-6729. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  42. Shapiro B, Drummond AJ, Rambaut A, Wilson MC, Matheus PE, Sher AV, Pybus OG, Gilbert MT, Barnes I, Binladen J, et al.: Rise and fall of the Beringian steppe bison.

    Science 2004, 306(5701):1561-1565. PubMed Abstract | Publisher Full Text OpenURL

  43. Haak W, Forster P, Bramanti B, Matsumura S, Brandt G, Tanzer M, Villems R, Renfrew C, Gronenborn D, Alt KW, et al.: Ancient DNA from the first European farmers in 7500-year-old Neolithic sites.

    Science 2005, 310(5750):1016-1018. PubMed Abstract | Publisher Full Text OpenURL

  44. Leonard JA, Vila C, Fox-Dobbs K, Koch PL, Wayne RK, Van Valkenburgh B: Megafaunal extinctions and the disappearance of a specialized wolf ecomorph.

    Curr Biol 2007, 17(13):1146-1150. PubMed Abstract | Publisher Full Text OpenURL

  45. Leonard JA, Wayne RK, Cooper A: Population genetics of ice age brown bears.

    Proc Natl Acad Sci USA 2000, 97(4):1651-1654. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  46. Bökönyi S: History of domestic mammals in central and eastern Europe. Budapest: Akadémiai Kiadó; 1974.

  47. Ho SY, Larson G, Edwards CJ, Heupink TH, Lakin KE, Holland PW, Shapiro B: Correlating Bayesian date estimates with climatic events and domestication using a bovine case study.

    Biol Lett 2008, 4(4):370-374. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  48. Slatkin M, Hudson RR: Pairwise comparisons of mitochondrial DNA sequences in stable and exponentially growing populations.

    Genetics 1991, 129(2):555-562. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  49. Valdiosera CE, Garcia N, Anderung C, Dalen L, Cregut-Bonnoure E, Kahlke RD, Stiller M, Brandstrom M, Thomas MG, Arsuaga JL, et al.: Staying out in the cold: glacial refugia and mitochondrial DNA phylogeography in ancient European brown bears.

    Mol Ecol 2007, 16(24):5140-5148. PubMed Abstract | Publisher Full Text OpenURL

  50. Provan J, Bennett KD: Phylogeographic insights into cryptic glacial refugia.

    Trends in ecology & evolution (Personal edition) 2008, 23(10):564-571. PubMed Abstract | Publisher Full Text OpenURL

  51. Sommer RS, Fahlke JM, Schmolcke U, Beneke N, Zachos FE: Quaternary history of the European roe deer Capreolus capreolus .

    Mammal Rev 2009, 39:1-16. Publisher Full Text OpenURL

  52. Liew P-M, Kuo C-M, Huang S, Tseng M-H: Vegetation change and terrestrial carbon storage in eastern Asia during the Last Glacial Maximum as indicated by a new pollen record from central Taiwan.

    Global and Planetary Change 1998, 16-17:85-94. Publisher Full Text OpenURL

  53. Stock F, Edwards CJ, Bollongino R, Finlay EK, Burger J, Bradley DG: Cytochrome b sequences of ancient cattle and wild ox support phylogenetic complexity in the ancient and modern bovine populations.

    Anim Genet 2009, 40(5):694-700. PubMed Abstract | Publisher Full Text OpenURL

  54. Gotherstrom A, Anderung C, Hellborg L, Elburg R, Smith C, Bradley DG, Ellegren H: Cattle domestication in the Near East was followed by hybridization with aurochs bulls in Europe.

    Proc Biol Sci 2005, 272(1579):2345-2350. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  55. Boscato P: Vado all'Arancio (Massa Marittima, GR), studio delle faune.

    Rassegna di Archeologia 1996, 13:159-176. OpenURL

  56. Calattini M: Il livello finale di Grotta delle Mura- Area A.

    Rivista di Scienze Preistoriche, Supplemento 1: Askategi, miscellanea in memoria di Georges Laplace 2005, 305-320. OpenURL

  57. Palma di Cesnola A: Il Paleolitico Superiore in Italia.

    Garlatti e Razzai Editori, Firenze, Italy 1993, 38:645-679. OpenURL

  58. Poinar HN, Stankiewicz BA: Protein preservation and DNA retrieval from ancient tissues.

    Proc Natl Acad Sci USA 1999, 96(15):8426-8431. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  59. Cooper A, Poinar HN: Ancient DNA: do it right or not at all.

    Science 2000, 289(5482):1139. PubMed Abstract | Publisher Full Text OpenURL

  60. Caramelli D, Milani L, Vai S, Modi A, Pecchioli E, Girardi M, Pilli E, Lari M, Lippi B, Ronchitelli A, et al.: A 28,000 years old Cro-Magnon mtDNA sequence differs from all potentially contaminating modern sequences.

    PLoS ONE 2008, 3(7):e2700. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  61. Bandelt HJ, Forster P, Rohl A: Median-joining networks for inferring intraspecific phylogenies.

    Mol Biol Evol 1999, 16(1):37-48. PubMed Abstract | Publisher Full Text OpenURL

  62. Excoffier L, Smouse PE, Quattro JM: Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data.

    Genetics 1992, 131(2):479-491. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  63. Excoffier L, Laval G, Schneider S: Arlequin (version 3.0): An integrated software package for population genetics data analysis.

    Evolutionary Bioinformatics 2005, 1:47-50. OpenURL

  64. Drummond AJ, Rambaut A: BEAST: Bayesian evolutionary analysis by sampling trees.

    BMC Evol Biol 2007, 7:214. PubMed Abstract | BioMed Central Full Text | PubMed Central Full Text OpenURL

  65. Posada D, Crandall KA: MODELTEST: testing the model of DNA substitution.

    Bioinformatics 1998, 14(9):817-818. PubMed Abstract | Publisher Full Text OpenURL

  66. Rambaut A, Drummond AJ: TRACER. Version 1.3. [http://beast.bio.ed.ac.uk/Tracer] webcite

    2004.

  67. Ho SY, Phillips MJ, Cooper A, Drummond AJ: Time dependency of molecular rate estimates and systematic overestimation of recent divergence times.

    Mol Biol Evol 2005, 22(7):1561-1568. PubMed Abstract | Publisher Full Text OpenURL

  68. Anderson CN, Ramakrishnan U, Chan YL, Hadly EA: Serial SimCoal: a population genetics model for data from multiple populations and points in time.

    Bioinformatics 2005, 21(8):1733-1734. PubMed Abstract | Publisher Full Text OpenURL

  69. Drummond AJ, Rambaut A, Shapiro B, Pybus OG: Bayesian coalescent inference of past population dynamics from molecular sequences.

    Mol Biol Evol 2005, 22(5):1185-1192. PubMed Abstract | Publisher Full Text OpenURL