Past climate changes, population dynamics and the origin of Bison in Europe
- Diyendo Massilani†1,
- Silvia Guimaraes†1,
- Jean-Philip Brugal2, 3,
- E. Andrew Bennett1,
- Malgorzata Tokarska4,
- Rose-Marie Arbogast5,
- Gennady Baryshnikov6,
- Gennady Boeskorov7,
- Jean-Christophe Castel8,
- Sergey Davydov9,
- Stéphane Madelaine10,
- Olivier Putelat11, 12,
- Natalia N. Spasskaya13,
- Hans-Peter Uerpmann14,
- Thierry Grange†1Email author and
- Eva-Maria Geigl†1Email authorView ORCID ID profile
© Massilani et al. 2016
Received: 20 July 2016
Accepted: 11 October 2016
Published: 21 October 2016
Climatic and environmental fluctuations as well as anthropogenic pressure have led to the extinction of much of Europe’s megafauna. The European bison or wisent (Bison bonasus), one of the last wild European large mammals, narrowly escaped extinction at the onset of the 20th century owing to hunting and habitat fragmentation. Little is known, however, about its origin, evolutionary history and population dynamics during the Pleistocene.
Through ancient DNA analysis we show that the emblematic European bison has experienced several waves of population expansion, contraction, and extinction during the last 50,000 years in Europe, culminating in a major reduction of genetic diversity during the Holocene. Fifty-seven complete and partial ancient mitogenomes from throughout Europe, the Caucasus, and Siberia reveal that three populations of wisent (Bison bonasus) and steppe bison (B. priscus) alternately occupied Western Europe, correlating with climate-induced environmental changes. The Late Pleistocene European steppe bison originated from northern Eurasia, whereas the modern wisent population emerged from a refuge in the southern Caucasus after the last glacial maximum. A population overlap during a transition period is reflected in ca. 36,000-year-old paintings in the French Chauvet cave. Bayesian analyses of these complete ancient mitogenomes yielded new dates of the various branching events during the evolution of Bison and its radiation with Bos, which lead us to propose that the genetic affiliation between the wisent and cattle mitogenomes result from incomplete lineage sorting rather than post-speciation gene flow.
The paleogenetic analysis of bison remains from the last 50,000 years reveals the influence of climate changes on the dynamics of the various bison populations in Europe, only one of which survived into the Holocene, where it experienced severe reductions in its genetic diversity. The time depth and geographical scope of this study enables us to propose temperate Western Europe as a suitable biotope for the wisent compatible with its reintroduction.
KeywordsAncient DNA Bison Population dynamics Evolution Climate Paleoenvironment Next generation sequencing Sequence capture
Drastic climatic fluctuations during the Pleistocene in the northern hemisphere led to population contractions, extinctions, re-expansions, and colonizations of fauna and flora . Bison, along with other large ungulates, thrived during the middle and late Pleistocene . Numerous cave paintings and engravings in France and Spain, such as those in the caves of Chauvet, Lascaux, and Altamira, attest to the important role this impressive animal played for the late Paleolithic hunter-gatherers. The steppe bison (Bison priscus (Bojanus, 1827)) appears in the fossil record during the early middle Pleistocene, replacing another archaic but smaller forest-adapted bison (B. schoetensacki ), which went extinct ca. 700 kiloyears ago (kya) . Since B. priscus was adapted to the cold tundra-steppe, occurrence of its remains is considered indicative of open environments . It roamed over Europe and Asia, and also crossed the Bering Strait during the middle Pleistocene to populate North America, where it evolved into the American Bison B. bison [3, 6, 7]. The numerous fossil remains display a pronounced sexual dimorphism, and a large initial body size, gradually decreasing throughout the Pleistocene . Differences in morphology related to climatic, environmental, and topographic conditions have led several authors to propose a high diversity for the Pleistocene cold-steppe bison expressed as subspecies or ecomorphotypes [5, 6].
The taxonomy, evolutionary history, and paleobiogeography of the genus Bison in Eurasia, and of the European bison or wisent B. bonasus (Linnaeus, 1758) in particular, is still patchy despite a rich fossil record and its current endangered status (e.g., [8–11]). Indeed, two opposing hypotheses on the evolution of bison in Eurasia coexist . Traditionally, it has been considered that bison developed within one single phylogenetic line (B. schoetensacki – B. priscus – B. bonasus), but it has also been proposed that at least two parallel lines of bison existed, one being the line of forest bison from B. schoetensacki (Freudenberg, 1910) to the recent B. bonasus and the other being the line of the steppe bison B. priscus (for a review see ). Thus, the phyletic relationships between B. schoetensacki, B. priscus, and B. bonasus, as well as the approximate date and geographical origin of the wisent, remain elusive, due in part to the limited power of paleontological studies to resolve species-level mammalian taxonomy issues or to detect broad-scale genetic transitions at the population level .
B. priscus disappeared from the fossil record of Western Europe at the end of the Pleistocene, around 12–10 kya, and relict populations of B. priscus seem to have survived until the beginning of the middle Holocene (7–6 kya) in Siberia (e.g., [13, 14]). In Europe, B. priscus is believed to have been replaced at the end of the Pleistocene or during the Holocene by the morphologically (eidonomically) distinguishable wisent B. bonasus [2, 10, 15, 16]. At least two sub-species are recognized: (1) B. b. bonasus Linnaeus, 1758 from the Lithuanian lowland and the Polish Białowieża ecosystem, and (2) the Caucasian highland B. b. caucasicus (Turkin and Satunin, 1904) . B. priscus was adapted to forest-steppe and steppe, and B. bonasus to forest and mountain-forest environments. B. priscus and B. bonasus are anatomically much closer to each other than to other more ancient bison, such as B. schoetensacki. B. bonasus has a relatively more massive rear quarter and shorter horns compared to B. priscus, which has longer and slightly curved horns and a smooth double-humped appearance [15, 16]. B. priscus and B. bison (Linnaeus, 1758), both of which are grazers, have a lower head position than B. bonasus, which is a mixed feeder . It is, however, very difficult to assign fossil bison bones to either species . Paleolithic paintings of bison from caves in France and Spain are often classified as belonging to either B. priscus or B. bonasus . The diversity of the cave art depictions and the large range of their occurrence is interpreted as indicating an origin of B. bonasus in the area between southern Europe and the Middle East and of its existence well before the end of the late Pleistocene at a time when B. priscus was still present .
Both extant bison species narrowly escaped extinction. The American B. bison was almost wiped out during the 19th century through commercial hunting and slaughter, but also due to introduced bovine diseases and competition with domestic livestock . The wisent also almost went extinct at the beginning of the 20th century. Indeed, similar to other large herbivores, such as the aurochs, intensification of agriculture since the Neolithic period pushed the wisent into the forests of Eastern Europe , where it was strictly protected for several centuries as royal game . During the First World War, however, a diminished population size followed by poaching led to its extinction in the wild . The entire population living today is descended from just 12 out of the 54 surviving animals at the beginning of the 1920s .
The wisent is still poorly characterized genetically. While genetic markers from the autosomes and Y chromosomes of American bison and wisent are closer to each other than to the other members of the genus Bos and they can reproduce and give rise to fertile offspring, their mitochondrial genomes are phylogenetically separated [9, 21, 22]. Indeed, mitochondrial sequences of the American bison and the yak Poephagus mutus f. grunniens (Linnaeus, 1758) form a distinct cluster, while the wisent occupies a phylogenetic position closer to Bos primigenius f. taurus (Linnaeus, 1758), a phenomenon that has been explained by incomplete lineage sorting or ancient hybridization [21, 22]. European, Siberian, and American B. priscus mitogenomes were shown to be phylogenetically closer to B. bison than to B. bonasus [7, 14, 23].
Ancient DNA studies have the potential to better resolve taxonomy than paleontological studies, in particular at the species level, and have revealed a far more dynamic picture of megafaunal communities, biogeography, and ecology, including repeated localized extinctions, migrations, and replacements (e.g.,  and citations therein). To better understand the phylogeography and evolution of the wisent and to reconstruct its origin, we performed a study of the mitochondrial genome over the past 50,000 years (50 kyr) across Europe and Asia. Using Pleistocene skeletal remains from Siberia (Yakutia), the Caucasus and Western Europe (France, Switzerland), as well as Neolithic and Medieval samples from France and pre-bottleneck B. bonasus samples from the early 20th century from the Polish lowland and the Caucasian highland lines, we constructed a phylogenetic framework based on both the hypervariable region (HVR) of the mitochondrial DNA as well as complete mitochondrial genomes of selected specimens. Our results give new insights into the climate-driven dynamics of the bison in Europe.
The evolution of the mitochondrial HVR
The two other clades, Bb1 and Bb2, are more closely related to the extant wisent B. bonasus. The Bb1 clade, significantly divergent from Bb2, contains most of the specimens from the Mezmaiskaya cave northwest of the Caucasus in Russia, dating from ca. 48–43 kya, as well as seven specimens from southern France (Arquet, Les Plumettes), dating from ca. 47–34 kya. Thus, the geographical range of the Bb1 clade stretched from at least the northern Caucasus to the south of France, from 48 to 34 kya during the overall mild Marine Isotope Stage (MIS) 3 (ranging from 57 to 29 kya , http://www.lorraine-lisiecki.com/LR04_MISboundaries.txt). None of the haplotypes of the Bb1 clade was found in the more recent specimens, suggesting that the corresponding population may have been the first to disappear. This Bb1 population in France was apparently replaced by the steppe bison of the Bp clade around 30 kya, the latter remaining there throughout the glacial period of MIS2 (ranging from 29 to 14 kya , http://www.lorraine-lisiecki.com/LR04_MISboundaries.txt).
Clade Bb2 includes both ancient specimens from the Caucasus and Western Europe as well as all recent and extant B. bonasus specimens. Nearly all ancient samples belonging to this clade are distinct from the more recent populations and include a ca. 49 kyr-old specimen from the Mezmaiskaya cave, two specimens from the Kudaro cave in the central part of the southern slope of the Greater Caucasus dated at 38 and 22 kya, and specimens from Western Europe – one from the Kesslerloch cave (Switzerland) dated at 14 kya, and, in France, two bison from the upper stratigraphic sequence of the Gral dated at 12 kya, two 5.2-kyr-old bison samples from the Neolithic site of Chalain, as well as three medieval (7th to 8th century CE) specimens from Alsace. The members of this clade represent the western European Pleistocene-Holocene lineage of B. bonasus and display a high mitochondrial diversity. This lineage appears to have replaced the B. priscus lineage, at least in France, at the end of the Upper Pleistocene between 15 and 12 kya, coinciding with the onset of a more temperate climate, and to have persisted in France up to the Middle Ages. Apart from the sequence found in the sample from Kesslerloch (12.2 kya), none of the ancient Bb2 sequences are present in the extant mitochondrial gene pool.
Within the Bb2 clade, a compact group of closely related sequences, are the Upper Pleistocene Kesslerloch specimen and the 1898–1917 pre-bottleneck wisents from both Poland and the Caucasus that almost became extinct at the end of the First World War. Five out of 24 of these pre-bottleneck bison have a HVR sequence identical to that of extant B. bonasus, whereas the rest differ from the extant sequence by only one or two single nucleotide polymorphisms (SNPs). Thus, the pre-bottleneck mitochondrial diversity appears only slightly higher than at present and much lower than that observed in older samples. The 14-kyr-old Kesslerloch specimen reveals the first occurrence of the mitogenome lineage of the extant B. bonasus population. Thus, the modern population corresponds to a minor fraction of the diversity that was present in Europe during the Late Pleistocene when B. bonasus replaced B. priscus.
A major reduction in the intrapopulation diversity is apparent from the Late Pleistocene to the early 20th century (Additional file 1: Table S2). The Pleistocene populations of Siberia, the Caucasus, and Europe are characterized by a high diversity at both the haplotype (H = 1.00) and nucleotide levels (as estimated with Pi and various Theta estimators), the nucleotide diversity of the European B. priscus being lower than that of its Siberian population. The Western European Holocene population of B. bonasus experienced a major reduction of its diversity between the Middle Ages and the beginning of the 20th century.
Node age estimates and clock rate estimates obtained through Bayesian analyses
Age estimates of the nodes (kya)
Age from HVR (95 % HPD)
Age from mitogenome (95 % HPD)
Bos Taurus/B. bonasus
B. indicus/B. taurus
B. grunniens/B. priscus
All B. bonasus
B. bonasus Bb2
Clock rate estimates (per site * year)
Rate from HVR (95 % HPD)
Rate from mitogenome (95 % HPD)
5.4 (3.9–6.9) × 10–7
4.6 (3.7–5.6) × 10–7
2.5 (2.1–3.0) × 10–8
Coding (first to second positions)
2.4 (2.0–2.8) × 10–8
Coding (third position)
9.9 (8.5–11.4) × 10–8
The phylogenetic analyses of the B. bonasus haplotypes reveal that the ancient Caucasian population had a deep root and was highly diverse. Based on the HVR, the age of the MRCA of the Bb1 and Bb2 clades is estimated at 438 (643–284) kya, which is significantly older than the MRCA of B. priscus (Fig. 2 and Table 1). This observation indicates that the Upper Pleistocene Caucasian population had a deeper root than the contemporaneous and more northern B. priscus population, even though the latter may have occupied in the past a much larger territory from Western Europe to the American continent. The Bb2 clade encompasses several branching events, the separation of a branch represented by a ca. 50-kyr-old northern Caucasus specimen from Mezmaiskaya occurring first, i.e., 105 (158–75) kya, followed, 61 (89–43) kya, by the separation of a branch represented by a ca. 37.8-kyr-old specimen from the Kudaro cave in the southern Caucasus. Then, 42 (60–26) kya, the subgroup comprising a 22.2-kyr-old specimen from the Kudaro cave as well as all French specimens between 12.4 and ca. 1.2 kya separated from the haplogroup encompassing the modern B. bonasus sequences. Finally, in this latter haplogroup, the 14-kyr-old Kesslerloch specimen and the early 20th century Central European and Caucasus specimens have a MRCA estimated at 16 (33–14) kya. This sequential order of the radiation events and the previously mentioned reduction of the population diversity indicates that, at the transition between the Pleistocene and the Holocene, Western and Central Europe were populated by a subset of a wisent population established in an area of Southwest Asia including the Caucasus, and that only a part of this population survived into the 20th century.
The evolution of the Bison clades based on entire mitogenomes
Complete mitogenomes confirm the observation made with the HVR sequences that the European B. priscus specimens from the Late Pleistocene were more homogeneous genetically than the contemporaneous Siberian population. Indeed, the two ca. 14–19 kya French mitogenomes from the Gral (Gral232), and Trois-Frères cave  are very similar, with only 19 SNPs distinguishing them, whereas the Siberian mitogenomes (all three Yaku and Rauchua ) are much more divergent, with 8-fold more SNPs distinguishing them (ca. 150 SNPs separating Rauchua from the Yaku and ca. 50 SNPs separating the Yaku samples from each other). The more ancient, ca. 38 kya French sample of La Berbie (LBN6a) is, however, more distant from the two French B. priscus (ca. 50 SNPs) and more closely related to one of the Siberian samples (Yaku118, 23 SNPs). Strikingly, the age of the MRCA of the three Yakutian and three French B. priscus samples is estimated at 59 (66–52) kya. Since the MRCA of these sequences is also the MRCA of the group of HVR sequences comprising all other French B. priscus samples (Fig. 1), this indicates that the steppe bison inhabiting France between 39 and 14 kya originated from a migration from North-East Eurasia that occurred not earlier than 59 (66–52) kya. Presumably, the low genetic diversity of the Western European B. priscus population is due to a founder effect during this migration.
The phylogenetic relationships between the B. bonasus mitogenomes indicate that modern wisent in central Europe are more closely related to the Late Pleistocene/Holocene Western European population than to the Late Pleistocene Caucasian population. Indeed, the 14-kyr-old specimen from the site of Kesslerloch at the Swiss-German border is closely related to modern wisent, with only 15 SNPs distinguishing the ancient and modern mitogenome sequences. This demonstrates that the range of the ancestral population of the modern wisent encompassed Middle Europe 14 kya. The Late Pleistocene Caucasian population was highly diverse with a deep root separating most of the specimens from the Mezmaiskaya cave in the northern Caucasus, belonging to the Bb1 clade, from the specimens from the Kudaro cave in the southern Caucasus, belonging to the Bb2 clade and including a specimen from the northern Caucasus. About 300 SNPs distinguished the Bb1 (Mez127, Mez128) and the Bb2 (Mez130, Kud133, Kud136) mitogenomes and approximately 100 SNPs distinguished the northern (Mez 130) and southern Caucasian (Kud133, Kud136) mitogenomes of the Bb2 clade. In contrast, the members of the Bb2 clade in Western Europe were more similar, with only 39 SNPs distinguishing the French (Gral125) and Swiss (KSL) specimens and three SNPs distinguishing the two French samples (Gral76, Gral125), in agreement with the reduction of diversity observed when comparing the HVR of the Western European and Caucasian B. bonasus populations. Thus, complete mitogenomes confirm the observations from the HVR analyses of a larger sample size while providing a more accurate phylogenetic analysis, in particular with respect to the dating of the MRCA of the various mitogenomes.
European bison population turnover during the late Pleistocene and the Holocene
The marked climatic variations that occurred during the Pleistocene implied drastic environmental changes that triggered the profound reorganization of floral and faunal biomes at different geographic and temporal scales. Response of biota to climatic stimuli was regionally/continentally distinct, sometimes occurring synchronously but most often diachronously (e.g., [2, 28]). The repeated cyclical climate changes progressively promoted landscape renewal, locally modifying, obliterating, and creating specific and varied ecological niches . The various expansions and contractions of bison populations in response to these climatic fluctuations led to apparent alternating occupations of Western Europe by either B. priscus or B. bonasus during the last 50 kyr, each presumably respectively originating from either a northeastern territory including Siberia, or a southern territory including the Caucasus. We hypothesize that these alternating occurrences result from climate-driven changes of two different habitats to which B. priscus, a grazer, and B. bonasus, a mixed feeder, were adapted, i.e., the open tundra-steppe for the first and open woodlands for the second. Indeed, the diet of B. priscus during the LGM included typical steppe and grassland (C3) vegetation and lichens , whereas the wisent’s diet in the Holocene was more flexible and included a higher content of shrubs . Similarly, on the American continent, B. priscus adapted to the climatic and environmental changes of the Holocene and evolved in two recently divergent forms, the plain bison thriving on the grasslands of the Great Plains and the wood bison inhabiting the boreal forest in North America. This suggests that, on the Eurasian continent, the competition with B. bonasus, a species seemingly better adapted to a more temperate environment, may have prevented a similar adaptation of B. priscus to habitat changes during the warmer periods. Finally, in Western Europe, local variations in ecological and physical barriers could have affected the speed of bison population turnover as a reaction to climatic shifts. In the future, a higher resolution genetic study involving a much higher sample number with denser time and space sampling may provide a more accurate and nuanced view of these population turnovers.
The mitochondrial lineages of B. bonasus that were present 40 kya have an older root than those of B. priscus (246 (283–212) vs. 114 (134–95) kya). This indicates that the B. bonasus did not experience as severe a bottleneck prior to 100 kya as the population reduction the steppe Bison experienced between 150 and 100 kya, before later thriving in northern Eurasia, Beringia, and in North America between 80 and 20 kya. All the various population expansions and contractions that occurred in response to the climatic and environmental changes characterizing the late Pleistocene and the transition to the Holocene gave rise to a reduction of the population diversity and to the extinction of lineages, like the Bb1 lineages, which apparently was the first to disappear, and the Bp lineage, which disappeared at the beginning of the Holocene on the Eurasian continent. Finally, the reduction of the diversity of the Bb2 lineage seems to involve both migrations and local extinctions with the major and most severe reductions occurring between the Early Holocene and historic times, most likely owing to human impact through hunting pressure and habitat fragmentation.
Radiation of the Bos and Bison lineages
Since the initial classification of Linnaeus in 1758 of Bos and Bison within a single genus, there has been debate about whether they should be placed in separate genera . Indeed, the members of these genera can still be crossed, despite reduced fertility in some combinations of crosses, and B. p. taurus genetic material can be found in a number of present-day American bison and wisent [32, 33]. Thus, gene flow could have occurred between the various lineages, in particular in the early phases of their differentiation. Such gene flow between separated populations has already been considered as an alternative possibility to incomplete lineage sorting to explain that the mitogenome of B. p. taurus was closer to that of B. bonasus than to that of B. bison [21, 22]. The phylogenetic estimates of the age of the various common ancestors of these two lineages indicate, however, that the most parsimonious interpretation of incomplete lineage sorting suffices to account for the relative affiliations of these mitogenomes (Additional file 3: Figure S2). Indeed, as mentioned above, the age estimates indicate that the MRCA of the Bos and Bison mitogenome lineages is not much older than that of the B. p. taurus and B. bonasus lineages, which are 927 (1064–790) kya and 768 (886–657) kya, respectively, with a 35–40 % overlap of the confidence interval (95 % HPD). This indicates that the two bifurcation events have occurred within a relatively short evolutionary period. Rapid speciation during this short evolutionary period appears unlikely since these species have still not yet totally lost interfertility almost a million years later. Our radiation date estimates are within the same range of the earliest fossils that were clearly attributed to the Bos genus and that are dated at 1 mya . They are also consistent with the analysis of the complete genome of a modern wisent that estimated that the wisent and bovine species diverged between 1.7 and 0.85 mya through a speciation process involving an extended period of limited gene flow with some more recent secondary contacts posterior to 150 kya . Thus, incomplete lineage sorting of mitogenomes in a metapopulation of the Bos and Bison ancestors during the period of divergence of these species could account for the affiliation patterns of these mitogenomes without the need to postulate a more recent post-speciation gene flow. It is interesting to note that this radiation event could be coincident with the onset and intensification of high-latitude glacial cycles (100 kyr-periodicity) around 1.2–0.8 mya. Incomplete lineage sorting, however, does not preclude later sporadic introgression of nuclear DNA at any point up to the present owing to the persistent interfertility between these species, as evidenced by the detection in the wisent genome of ancient gene flow from B. p. taurus [35, 36].
The analysis of DNA preserved in ancient bison remains from Eurasia covering the last 50,000 years allowed us to retrace some of the population dynamics that took place during the Late Pleistocene and the Holocene, including population migrations, extinctions, and replacements. We could trace back the origin of the wisent to Eastern Europe, including the Caucasus. We observed several alternating expansion waves of wisent and steppe bison in Western Europe that were related with climatic fluctuations. The wisent was prevalent when the climate was more temperate leading to a more forested vegetation, similar to present-day Western Europe, which is compatible with its reintroduction in this area. In contrast, the steppe bison population originating from northern Eurasia was predominant in Western Europe during the colder periods of the Late Pleistocene with their open environments. These fluctuations may have been recorded in Paleolithic cave paintings, in particular in the cave of Chauvet that had been occupied by humans over a long period and where two distinct bison types are depicted.
Details of samples are described in Additional file 1: Table S1, Additional file 5: Document S1 and Additional file 6: Figure S4. All pre-PCR procedures were carried out in the high containment facility of the Jacques Monod Institute physically separated from areas where modern samples are analyzed and dedicated exclusively to ancient DNA analysis using the strict procedures for contamination prevention previously described .
The external surface of the specimens was removed with a sterile blade to minimize environmental contamination. For each bone sample, roughly 0.2 g was ground to a fine powder in a freezer mill (Freezer Mill 6750, Spex Certiprep, Metuchen, NJ), which was then suspended in 2 mL of extraction buffer containing 0.5 M EDTA, 0.25 M di-sodium hydrogen phosphate (Na2HPO4) at pH 8.0, 0.14 M 2-mercaptoethanol, and 0.25 mg/mL of proteinase K and incubated under agitation at 37 °C for 48 hours. Blank extractions were carried out for each of eight total extraction series.
Samples were then centrifuged and the supernatant was purified with the Qiaquick PCR Purification Kit (Qiagen, Hilden, Germany), as previously described .
PCR amplification of the HVR
Primers were designed to target conserved regions with minimal primer dimer propensity from a multiple alignment of all the bison mitochondrial sequences present in GenBank in 2012 using the software Oligo 7 as previously described [38, 39], and were then tested for efficiency and dimer formation using quantitative real-time PCR (qPCR) [38, 39]. Four primer pairs amplifying short (118–152 bp) subregions of the mitochondrial HVR were selected (Additional file 1: Table S6 and Additional file 2: Figure S1). To protect against cross-contamination between samples, we used the UNG-coupled quantitative real-time PCR system  and, to avoid the production of erroneous data due to the presence of bovine DNA in reagents, reagent decontamination was performed as previously described . Amplifications were performed in a final volume of 10 μL containing 2 mM MgCl2, 1 μM primers, 0.04 mM of dA/G/CTPs, 0.08 mM of dUTP, 0.01 U/μL of UNG, 1 U/μL of FastStart Taq (Roche Applied Science, Penzberg, Germany), and 1× qPCR home-made reaction buffer . Blank amplification controls were included for each amplification. In total, 415 amplification blanks were carried out during the various amplifications of 85 samples. No products were observed in any of the amplification blanks. Amplifications were performed using a LightCycler 1.5 (Roche Applied Science) with the following cycling program: 15 min at 37 °C (carry-over contamination prevention through digestion by UNG of dUTP-containing amplicons), 10 min at 95 °C (inactivation of UNG and activation of the Fast Start DNA polymerase) followed by 60 cycles at 95 °C for 15 sec, at 60 °C for 40 sec (for primer pairs Bon1, 2 and 3), and 95 °C for 15 sec, or 56 °C for 15 sec and 67 °C for 20 sec (for primer pair BB3r4m), followed by a final melting curve analysis step. All extracts were tested for inhibition as previously described . Each sequence was determined on both strands from at least two independent amplifications using capillary electrophoresis sequencing. Sequencing data analyses were performed using the software Geneious 6.1.8 .
Capture and sequencing of whole mitogenomes
Dual barcoded libraries for Illumina sequencing were constructed in the dedicated ancient DNA facility at the Jacques Monod Institute (Paris, France) using a double-stranded procedure previously described [37, 44]. Ancient DNA extracts were quantified on a Qubit 2.0 Fluorimeter (Thermo Fisher Scientific, Waltham, MA). Libraries were prepared using between 3 and 124 ng total DNA. The USER Enzyme, End-repair and adapter ligation mixtures were decontaminated using Ethidium monoazide  to inactivate bovine DNA associated with BSA present in these reagents that could interfere with the final results. Prior to library preparation, uracil residues were removed through treatment with 1.5 units/30 μL reaction volume of USER Enzyme (New England Biolabs, Ipswich, MA) for 60 min at 37 °C in NEBNext End Repair buffer (New England Biolabs). End repair was performed using the NEBNext End Repair Module, using only half of the enzyme mix recommended by the manufacturer. End-repaired DNA was purified on MinElute columns using the Qiagen Gel extraction protocol, and DNA was eluted twice, each with 16 μL of γ-irradiated H2O . Blunt-end ligations were performed using the NEBNext Quick Ligation Module (New England Biolabs); 1 μL each of N700 and N500 barcoded adaptors (20 μM) from the Nextera XT series (Illumina Inc., San Diego, CA) and 1.5 μL of T4 ligase were added to the purified reactions and incubated 30 min at 20 °C in ligation buffer. After ligation, elongation of the adapters was performed by adding 1 volume of OneTaq DNA polymerase 2× Master Mix (New England BiolabsI) and incubating for 20 min at 60 °C in an Eppendorf MasterCycler epGradientS (Eppendorf, Hamburg, Germany). A 3-μL aliquot of a 20 μM mix of primers IS7 and IS8 (Illumina) were added to the tubes after the elongation step and libraries were directly amplified by PCR. PCR was performed as follows: 95 °C for 2 min, followed by 9 amplification cycles (denaturation at 95 °C for 20 sec, primer annealing at 46 °C for 30 sec, and extension at 60 °C for 1.5 min). This amplification protocol was modified to reduce the amplification bias of GC and AT-rich DNA fragments ( and unpublished data). The purification of the libraries following nine cycles of amplification was carried out using 1.3× SPRI magnetic beads from the NucleoMag NGS clean-up and Size Select kit (Macherey Nagel, Düren, Germany). DNA was eluted with 50 μL of γ-irradiated H2O (irradiation with 2 kGy). This first eluate was purified a second time using the same protocol with SPRI magnetic beads and eluted in 30 μL of EB buffer. Amplified and purified libraries were visualized and quantified on a Bioanalyzer 2100 (Agilent, Santa Clara, CA) using a High Sensitivity DNA assay chip.
DNA sequence capture of mitochondrial genomes
For the production of complete mitogenome sequences, we developed a sequence capture approach using biotinylated RNA baits .
Capture: bait synthesis
The complete Bos taurus mitogenome was amplified using 11 overlapping PCR products (ca. 1.5 kb) containing the T3 RNA polymerase promoter sequence on the 5′ ends of either the forward or reverse strand. The 22 PCR reactions were performed in a final volume of 50 μL containing 12 ng of total genomic B. taurus DNA, 1 U of FastStart Taq, 1× FastStart buffer, 120 μM dATP, dCTP, dGTP and dUTP, 1 mg/mL BSA, and 0.4 μM of one of the 22 primer pairs (Additional file 1: Table S3). PCRs were performed on an Eppendorf MasterCycler epGradientS as follows: polymerase activation at 95 °C for 10 min, 45 amplification cycles (95 °C for 20 sec, 58 °C for 30 sec, 72 °C for 3 min), and a final extension step at 72 °C for 5 min. The amplification of each fragment was visualized by agarose gel electrophoresis in a 1 % agarose gel (100 volts for 30 min), and the products were purified using the QIAquick PCR purification kit. The purified mitochondrial fragments were then quantified on a Nanodrop ND-1000 spectrophotometer. The products were pooled in equimolar amounts in two different lots: one containing the fragments with the T3 RNA polymerase promotor sequence at the 5′ end of the forward strands, and the other with fragments containing the T3 RNA polymerase promotor sequence at the 5′ end of the reverse strands.
A total of 500 ng of each PCR product pool was used as a template in a 40 μL MEGAscript T3 transcription reaction (Ambion, Foster City, CA) containing 0.75 mM of each ATP, UTP and GTP, 0.375 mM of CTP, and 0.375 mM of Biotin-14-CTP (Invitrogen, Thermo Fisher Scientific). After an overnight incubation at 37 °C, the DNA template was digested with 2 U/40 μL reaction volume of TURBO DNase (Ambion) for 30 min at 37 °C. The ca. 1.5 kb transcripts were purified using an RNeasy MinElute Cleanup Kit (Qiagen) and eluted twice in 35 μL of nuclease-free water. The eluate was treated a second time with 5 U/40 μL reaction volume of Recombinant DNase1 (Roche) for 15 min at 35 °C in 1× Recombinant DNase1 buffer (Roche) in order to hydrolyze potentially remaining DNA template molecules. The DNase was removed through another purification of the transcripts using the RNeasy MinElute Cleanup Kit and the RNA was eluted as above. The purified ca. 1.5 kb transcripts were fragmented to generate RNA fragments ranging from 100 to 600 nucleotides with an average size of 300 nt using the NEBNext magnesium RNA Fragmentation Module kit (New England Biolabs) for 4 min at 94 °C. The fragmented transcripts were ethanol precipitated using three volumes of absolute ethanol and 0.3 M sodium acetate pH 5.5. The transcripts were washed with 500 μL of 70 % ethanol, dried and resuspended in 40 μL of nuclease-free water. The fragmented concentrated transcripts were visualized on a 2 % agarose gel, and the products were quantified using a Nanodrop ND-1000 spectrophotometer.
Capture: production of blocking oligos
PCR products corresponding to each N500 and N700 barcoded sequencing adaptor of the Nextera XT series (Illumina) containing the T7 RNA polymerase promotor sequence at the 5′ end of the forward strands were amplified in a 100 μL final volume reaction containing 100 pM of Nextera barcoded adaptor primer, 2 U FastStart Taq, 1× FastStart buffer, 120 μM of a mix of dNTPs, 1 mg/mL of BSA, and 0.5 μM of each N500 and N700 blocking oligo primers (N500 adaptor blocking oligos: T7BO + P5-F (5′ATGTAATACGACTCACTATAGGGAATGATACGGCGACCAC3′) and BO + P5-R (5′GGAAGAGCGTCGTGTAGG3′); N700 adaptor blocking oligos: T7BO + P7-F (5′ATGTAATACGACTCACTATAGGGAGATCGGAAGAGCACACG3′) and BO + P7-R (5′CAAGCAGAAGACGGCATAC3′). PCRs were performed on an Eppendorf MasterCycler epGradientS using a polymerase activation step at 95 °C for 10 min, followed by 25 amplification cycles (95 °C for 20 sec, 60 °C for 30 sec, 72 °C for 45 sec), and a final extension step at 72 °C for 2 min. PCR products were visualized on a 2 % agarose gel (100 volts for 1 hour), and then quantified using a Qubit 2.0 Fluorimeter. All PCR products were pooled in equimolar amounts and purified using the Qiagen Gel extraction protocol. DNA was eluted with 50 μL of nuclease-free water, and each purified pool was quantified using a Qubit 2.0 Fluorimeter; 50 ng of each PCR product pool were used as a template in a 20 μL MAXIscript T7 transcription (Ambion) containing 1 mM of each ATP, UTP, GTP, and CTP. After an overnight incubation at 37 °C, the DNA template was removed by a TURBO DNase treatment for 30 min at 37 °C. Unincorporated nucleotides, buffer components and DNase were removed by phenol-chloroform purification followed by alcohol precipitation. The blocking oligo transcript RNA was visualized on a 2 % agarose gel (100 volts for 1 hour), and quantified using a Nanodrop ND-1000 spectrophotometer.
Capture: hybrid selection
Each library was individually hybridized to the biotinylated probes in the presence of RNA blocking oligos complementary to the Illumina adapters for 48 hours at 62 °C. For each 38 μL final volume reaction, 15 μL of ancient DNA library (80–270 ng) was heated for 5 min at 95 °C, incubated for 5 min at 62 °C, and mixed with ca. 4 μg of each pool of blocking oligos (N500 or N700), prewarmed (62 °C) 2× hybridization buffer (10× SSPE, 10 mM EDTA and 0.2 % Tween-20), and ca. 500 ng prewarmed (2 min at 62 °C) mix of biotinylated RNA. After incubation, the hybridization mix was added to 500 ng (50 μL) M-280 streptavidin Dynabeads (Invitrogen), that had been washed three times and resuspended in 162 μL 1 M NaCl, 10 mM Tris-HCl, pH 7.5, 1 mM EDTA, and 1× Denhardt’s solution. After 40 min at room temperature, the beads were pulled down and washed once at room temperature for 15 min with 0.2 mL 1× SSC/0.1 % Tween-20, followed by three 10 min washes at 62 °C with 0.2 mL 0.1× SSC/0.1 % Tween-20, resuspending the beads once at each washing step. Hybrid-selected DNA was eluted twice with 50 μL 0.1 N NaOH. After 10 min at 20 °C, the beads were pulled down and the supernatant transferred to a tube containing 75 μL 1 M Tris-HCl, pH 7.5. The captured single stranded DNA fragments were then purified using the Qiagen Gel extraction protocol. DNA was eluted with 30 μL of EB buffer.
Capture: PCR amplification of enriched library
Amplification of the enriched library was performed using 25 μL of each library and FastStart Taq with 1× FastStart buffer, 100 μM each of dNTPs, and 0.4 μM of each Illumina primer IS7 and IS8 in a 50 μL final volume on an Eppendorf MasterCycler epGradient S. After an 8 min activation of the FastStart Taq at 95 °C, PCR was carried out for the appropriate number of cycles (95 °C for 20 sec, 46 °C for 30 sec, and 60 °C for 2 min, followed by a 2 min extension) to avoid plateau, as determined by a previous qPCR under similar conditions using 1 μL of each captured library in 10 μL total volume for 30 cycles. Amplified enriched libraries were then purified using the QIAquick PCR purification kit and quantified on a Nanodrop spectrophotometer.
A second round of capture was then performed as above, on all amplified and purified enriched libraries. After purification and amplification, the double-enriched libraries were then purified using 1.3× SPRI magnetic beads from the NucleoMag NGS clean-up and Size Select kit. DNA was eluted with 50 μL of EB buffer, and the eluates were purified a second time using the SPRI magnetic beads protocol with a final elution in 30 μL of EB buffer. The double-enriched amplified and purified libraries were visualized and quantified on a Bioanalyzer 2100 using a High Sensitivity DNA assay chip.
Enriched purified libraries were pooled in equimolar amounts, and paired-end sequenced using 2 × 75 cycles on a V3 flowcell for MiSeq Illumina platform at the Jacques Monod Institute (Paris, France) according to the manufacturer’s protocol.
Bioinformatic processing and mapping of the reads
Paired-end reads were merged using leeHom using the --ancientdna parameter . Merged reads were then mapped to modern Bison bison and Bison bonasus mitogenomes using bwa as described previously . The first 100 nt of each linearized reference mitogenome used for mapping was duplicated at the 3′ end to allow mapping of fragments overlapping the junction. Mapped read duplicates were then removed using samtools rmdup as previously described . The resulting bam files were then imported into Geneious 6.1.8  and remapped onto the appropriate mitogenome sequence without the 100 nt duplication (B. bison for B. priscus sequences, B. bonasus for the Bb1 and Bb2 sequences). Consensus sequences were generated in Geneious and verified by visual inspection of the aligned reads. Geneious was used to measure coverage depth and the number of covered bases displayed in Additional file 1: Table S4.
Sequence alignments were performed using the Muscle algorithm and were visually inspected and adjusted using Geneious 6.1.8 . The maximum likelihood analyses presented in Fig. 1 were computed using PHYML 3.0, using an HKY substitution model with a gamma-distributed rate of variation among sites (+G) and invariant sites (+I) . Robustness of the nodes was estimated using 500 bootstraps. RaXML 8.2.3 was used to generate the maximum likelihood bootstrap support values for the complete mitogenome alignment shown in Fig. 3 .
Phylogenetic analyses conducted under the Bayesian framework were performed using the program BEAST v. 1.8.2, which allows estimation of mutation and population history parameters simultaneously from temporally spaced sequence data . Nucleotide substitution models were chosen following comparisons performed with jModelTest 2.1.7 using the Bayesian Information Criteria . The HVR analysis presented in Fig. 2 was performed considering a TN93 model for the nucleotide substitution model, a gamma-distributed rate of variation among sites (+G) with four rate categories and invariant sites (i.e., TN93 + I + G model). For the complete mitogenome analysis presented in Fig. 3, we used four partitions, the HVR, the first and second positions of the codons within the coding region, the third position, and the RNA genes. We considered the HKY + I + G model for the first two partitions, and the TN93 + G for the last two. Default priors were used for all parameters of the nucleotide substitution model. For the analysis of Fig. 2, we used a strict molecular clock with a lognormal prior for the substitution rate (mean = –15.0, stdev = 1.4) corresponding to a median of 2 × 10–7 substitutions per site and per year (95 % HPD 1.3 × 10–8 to 3.2 × 10–6) based on the estimation for Bison HVR substitution rate . For the various partitions of the mitogenome of Fig. 3, we used estimates for the human mitogenome substitution rate to set the priors : lognormal priors: HVR, mean = –16.1 stdev = 2.0, corresponding to a median of 1.0 × 10–7 (2 × 10–9 to 5 × 10–6); RNA, mean = –18.65 stdev = 2.0, corresponding to a median of 8.0 × 10–9 (1.6 × 10–10 to 4 × 10–7); first and second positions, mean = –18.5 stdev = 2.0, corresponding to a median of 9.0 × 10–9 (1.8 × 10–10 to 4.7 × 10–7); third position, mean = –17.7 stdev = 2.0 corresponding to a median of 2.0 × 10–8 (4× 10–10 to 1 × 10–6). Finally, a standard coalescent model was considered for the tree prior with a Bayesian skyline plot to model populations (5 and 10 populations with default parameters for the HVR and the mitogenomes respectively). The prior for the tree height followed a log-normal distribution, mean = 14.6 stdev = 0.6, truncate to 8.0 × 106 and 1.0 × 105, corresponding to a median of 2.2 × 106 and a 95 % HPD of (6.3 × 106 to 6.7 × 105), which integrates the various fossil finds assumed to correspond to ancestors of cattle and bison [54, 55].
To estimate the posterior distribution of each parameter of interest, we used the Markov Chain Monte Carlo algorithm implemented in the BEAST software. We ran five independent chains with initial values sampled as described above and an input UPGMA tree constructed using a Juke-Cantor distance matrix. Each of these chains was run for 50,000,000 iterations and for each parameter of interest, 18,000 samples (one every 2500 generated ones) were drawn after discarding a 10 % burn-in period. The BEAST output was analyzed with the software Tracer v. 1.6 (http://tree.bio.ed.ac.uk/software/tracer/). Visual inspection of the traces and the estimated posterior distributions suggested that each MCMC had converged on its stationary distribution. Using Logcombiner v. 1.8.2, we further combined all the results from the five independent chains. The maximum clade credibility tree with the median height of the nodes was finally calculated using TreeAnnotator v. 1.8.2 and visualized using FigTree v. 1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).
We thank G Baldacci for continuous support. We thank Lou Saïer for help with the production of HVR data for some samples. We like to thank Marie-Hélène Moncel, Camille Daujeard, and Marie-Anne Julien for providing samples that did not yield genetic results. We acknowledge the contribution to sampling of Myriam Boudadi-Maligne, Jean-Baptiste Mallye, Pierre Pétrequin, Sylvie Lourdaux, and René Rémond. We are grateful to the French ministry of culture, Centre national de Préhistoire for providing the image of the wall paintings of the Chauvet cave.
The study was supported by the French national research center CNRS. The paleogenomic facility obtained support from the University Paris Diderot within the program “Actions de recherche structurantes”. The sequencing facility of the Institut Jacques Monod, Paris, is supported by grants from the University Paris Diderot, the Fondation pour la Recherche Médicale (DGE20111123014), and the Région Ile-de-France (11015901). MT acknowledges financial support through the “Biodiversity of East-European and Siberian large mammals on the level of genetic variation of populations (BIOGEAST)” project within the 7th European Framework Programme, Marie Curie Actions, that allowed sampling to take place in Russia.
Availability of data and material
The produced DNA sequences have been deposited at NCBI: From KX870126 to KX870182 and from KX898005 to KX898020.
EMG and TG designed and supervised the overall research with an initial input from MT. DM developed sequence capture and produced the mitogenome data. SG produced the HVR data. TG, EMG, DM, and SG analyzed and interpreted the data. TG, EMG, and EAB wrote the manuscript with inputs from DM and SG. JPB provided samples, discussed the data and corrected the manuscript. MT, RMA, GBa, GBo, JCC, SD, SM, OP, NS, and HPU provided samples and feedbacks on the manuscript. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Not applicable since no data from individuals were used.
Ethics approval and consent to participate
Not applicable since we did not use biological material.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Hewitt GM. Genetic consequences of climatic oscillations in the Quaternary. Phil Trans RL Soc Lond B. 2004;359:183–95.View ArticleGoogle Scholar
- Markova AK, Puzachenko AY, van Kolfschoten T, Kosintsev PA, Kuznetsova TV, Tikhonov AN, Bachura OP, Ponomarev DV, van der Pflicht J, Kuitems M. Changes in the Eurasian distribution of the musk ox (Ovibos moschatus) and the extinct bison (Bison priscus) during the last 50 ka BP. Quatern Int. 2015;378:99–110.View ArticleGoogle Scholar
- Guthrie RD. Frozen fauna of the Mammoth Steppe. The story of Blue Babe. Chicago and London: The University of Chicago Press; 1990.Google Scholar
- Brugal JP. Le Bison (Bovinae, Artiodactyla) du gisement Pléistocène moyen ancien de Durfort (Gard, France). Bull Mus Nat Hist Nat Section C: Sciences de la Terre. 1995;16(2-4):349–81.Google Scholar
- Brugal JP, David F, Enloe J, Jaubert J. Le Bison, gibier et moyen de subsistance des hommes du Paléolithique aux Paléoindiens des grandes plaines. Antibes: APDCA Editions; 1999.Google Scholar
- Flerow CC. Systematics and Evolution. In: Sokolov VE, editor. European Bison –Morphology, Systematics, Evolution, Ecology. [in Russian]. Moscow: Nauka (USSR Ac. of Sc.); 1979. p. 9–127.Google Scholar
- 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–5.PubMedView ArticleGoogle Scholar
- Groves CP. Systematic relationship in the Bovini (Artiodactyla, Bovidae). Z Zool Syst Evolutionforsch. 1981;19(4):264–78.View ArticleGoogle Scholar
- Hassanin A, An J, Ropiquet A, Nguyen TT, Couloux A. Combining multiple autosomal introns for studying shallow phylogeny and taxonomy of Laurasiatherian mammals: Application to the tribe Bovini (Cetartiodactyla, Bovidae). Mol Phylogenet Evol. 2013;66(3):766–75.PubMedView ArticleGoogle Scholar
- Kowalski K. The evolution and fossil remains of the European Bison. Acta Theriol. 1967;12(21):335–8.View ArticleGoogle Scholar
- Krasińska M, Krasiński ZA. European Bison: The Nature Monograph. 2nd ed. Heidelberg: Springer; 2013.View ArticleGoogle Scholar
- Cooper A, Turney C, Hughen KA, Brook BW, McDonald HG, Bradshaw CJ. Paleoecology. Abrupt warming events drove Late Pleistocene Holarctic megafaunal turnover. Science. 2015;349(6248):602–6.PubMedView ArticleGoogle Scholar
- Boeskorov GG, Potapova OR, Mashchenko EN, Protopopov AV, Kuznetsova TV, Agenbroad L, Tikhonov AN. Preliminary analyses of the frozen mummies of mammoth (Mammuthus primigenius), bison (Bison priscus) and horse (Equus sp.) from the Yana-Indigirka Lowland, Yakutia, Russia. Integr Zool. 2014;9(4):471–80.PubMedView ArticleGoogle Scholar
- Kirillova IV, Zanina OG, Chernova OF, Lapteva EG, Trofimova SS, Lebedev VS, Tiunov AV, Soares AER, Shidlovskiy FK, Shapiro B. An ancient bison from the mouth of the Rauchua River (Chukotka, Russia). Quatern Res. 2015;84(2):232–45.View ArticleGoogle Scholar
- DegerbØl M, Iversen J. The Bison in Denmark. A zoological and geological investigation of the finds in Danish Pleistocene deposits. Copenhagen: Danm. Geol. UndersØgelse; 1945.Google Scholar
- Geist V. The relation of social evolution and dispersal in ungulates during the Pleistocene, with emphasis on the old world Deer and the genus Bison. Quat Res. 1971;1(3):283–21.View ArticleGoogle Scholar
- Rautian GS, Kalabushkin BA, Nemtsev AS. A new subspecies of the European Bison, Bison bonasus montanus ssp. nov. (Bovidae, Artiodactyla). Dokl Biol Sci. 2000;375:636–40.PubMedView ArticleGoogle Scholar
- Bocherens H, Hofman-Kaminska E, Drucker DG, Schmolcke U, Kowalczyk R. European bison as a refugee species? Evidence from isotopic data on Early Holocene bison and other large herbivores in northern Europe. PLoS One. 2015;10(2), e0115090.PubMedPubMed CentralView ArticleGoogle Scholar
- Spassov N, Stovtechev T. On the origin of the wisent, Bison bonasus (Linnaeus, 1758): presence of the wisent in the upper Palaeolithic art of Eurasia. Advances in Paleontology “Hen to Panta”, papers in honour of C. Radulescu and P.M. Samson. Bucharest; 2003. pp. 125–30Google Scholar
- Potter BA, Gerlach SC, Gates CC, Boyd DP, Oetelaar GA, Shaw JS. History of Bison in North America. In: Gates CC, Freese CH, Gogan PJP, Kotzman M, editors. American Bison: Status Survey and Conservation Guidelines 2010. Gland: IUCN; 2010. p. 5-12.Google Scholar
- Verkaar EL, Nijman IJ, Beeke M, Hanekamp E, Lenstra JA. Maternal and paternal lineages in cross-breeding bovine species. Has wisent a hybrid origin? Mol Biol Evol. 2004;21(7):1165–70.PubMedView ArticleGoogle Scholar
- Zeyland J, Wolko L, Lipinski D, Wozniak A, Nowak A, Szalata M, Bocianowski J, Slomski R. Tracking of wisent-bison-yak mitochondrial evolution. J Appl Genet. 2012;53(3):317–22.PubMedPubMed CentralView ArticleGoogle Scholar
- Marsolier-Kergoat MC, Palacio P, Berthonaud V, Maksud F, Stafford T, Begouen R, Elalouf JM. Hunting the extinct Steppe Bison (Bison priscus) mitochondrial genome in the Trois-Freres paleolithic painted cave. PLoS One. 2015;10(6), e0128267.PubMedPubMed CentralView ArticleGoogle Scholar
- Lisiecki LE, Raymo ME. A Pliocene-Pleistocene stack of 57 globally distributed benthic δ18O records. Paleoceanography. 2005;20(1):1–17.View ArticleGoogle Scholar
- Andersen KK, Azuma N, Barnola JM, Bigler M, Biscaye P, Caillon N, Chappellaz J, Clausen HB, Dahl-Jensen D, Fischer H, et al. High-resolution record of Northern Hemisphere climate extending into the last interglacial period. Nature. 2004;431(7005):147–51.PubMedView ArticleGoogle Scholar
- Van Andel TH, Tzedakis PC. Palaeolithic landscapes of Europe and environs, 150,000-25,000 years ago: an overview. Quat Sci Rev. 1996;15(5–6):481–500.View ArticleGoogle Scholar
- Guiot J, Pons A, de Beaulieu JL, Reille M. A 140,000 year climatic reconstruction from two European pollen records. Nature. 1989;338:309–13.View ArticleGoogle Scholar
- Palombo MR. Deconstructing mammal dispersals and faunal dynamics in SW Europe during the Quaternary. Quat Sci Rev. 2014;96:50–71.View ArticleGoogle Scholar
- Julien M-F, Bocherens H, Burke A, Drucker DG, Patou-Mathis M, Krotova O, Péan S. Were European steppe bison migratory? 18O, 13C and Sr intra-tooth isotopic variations applied to a palaeoethological reconstruction. Quarternary Int. 2012;271:106–19.View ArticleGoogle Scholar
- Quiles A, Valladas H, Bocherens H, Delque-Kolic E, Kaltnecker E, van der Plicht J, Delannoy JJ, Feruglio V, Fritz C, Monney J, et al. A high-precision chronological model for the decorated Upper Paleolithic cave of Chauvet-Pont d’Arc, Ardeche, France. Proc Natl Acad Sci U S A. 2016;113(17):4670–5.PubMedView ArticlePubMed CentralGoogle Scholar
- Boyd DP, Wilson GA, Gates CC. Taxonomy and nomenclature. In: Gates CC, Freese CH, Gogan PJP, Kotzman M, editors. American Bison: Status Survey and Conservation Guidelines 2010. Gland: IUCN; 2010. p. 13–8.Google Scholar
- Yudin NS, Kulikov IV, Gunbin KV, Aitnazarov RB, Kushnir AV, Sipko TP, Moshkin MP. Detection of mitochondrial DNA from domestic cattle in European bison (Bison bonasus) from the Altai Republic in Russia. Anim Genet. 2012;43(3):362.PubMedView ArticleGoogle Scholar
- Halbert ND, Derr JN. A comprehensive evaluation of cattle introgression into US federal bison herds. J Hered. 2007;98(1):1–12.PubMedView ArticleGoogle Scholar
- Martínez-Navarro B, Rook L, Papini M, Libsekal Y. A new species of bull from the Early Pleistocene paleoanthropological site of Buia (Eritrea): parallelism on the dispersal of the genus Bos and the Acheulian culture. Quatern Int. 2010;212(2):169–75.View ArticleGoogle Scholar
- Gautier M, Moazami-Goudarzi K, Leveziel H, Parinello H, Groh C, Riall S, Kowalczyk R, Flori L. Deciphering the wisent demographic and adaptive histories from individual whole-genome sequences. Mol Biol Evol. 2016;33(11):2801–14.PubMedPubMed CentralView ArticleGoogle Scholar
- Wecek K, Hartmann S, Paijmans JLA, Taron U, Xenikoudakis G, Cahill JA, Heintzman PD, Shapiro B, Baryshnikov G, Bunevich AN, et al. Complex admixture preceded and followed the extinction of wisent in the wild. bioRxiv. 2016. doi: http://dx.doi.org/10.1101/059527.
- Bennett EA, Massilani D, Lizzo G, Daligault J, Geigl EM, Grange T. Library construction for ancient genomics: single strand or double strand? BioTechniques. 2014;56(6):289–90. 292–86, 298, passim.PubMedView ArticleGoogle Scholar
- Cote NM, Daligault J, Pruvost M, Bennett EA, Gorge O, Guimaraes S, Capelli N, Le Bailly M, Geigl EM, Grange T. A new high-throughput approach to genotype ancient human gastrointestinal parasites. PLoS One. 2016;11(1), e0146230.PubMedPubMed CentralView ArticleGoogle Scholar
- Guimaraes S, Pruvost M, Daligault J, Stoetzel E, Côté NM-L, Bennett EA, Nicolas V, Lalis A, Denys C, Geigl E-M, et al. A cost-effective high-throughput metabarcoding approach powerful enough to genotype ~44,000 year-old rodent remains from Northern Africa. Mol Ecol Res. 2016. Ahead of print. doi:10.1111/1755-0998.12565.
- Pruvost M, Grange T, Geigl E-M. Minimizing DNA contamination by using UNG-coupled quantitative real-time PCR on degraded DNA samples: application to ancient DNA studies. BioTechniques. 2005;38(4):569–75.PubMedView ArticleGoogle Scholar
- Champlot S, Berthelot C, Pruvost M, Bennett EA, Grange T, Geigl EM. An efficient multistrategy DNA decontamination procedure of PCR reagents for hypersensitive PCR applications. PLoS ONE. 2010;5(9), e13042.PubMedPubMed CentralView ArticleGoogle Scholar
- Pruvost M, Schwarz R, Bessa Correia V, Champlot S, Braguier S, Morel N, Fernandez-Jalvo Y, Grange T, Geigl E-M. Freshly excavated fossil bones are best for amplification of ancient DNA. Proc Natl Acad Sci U S A. 2007;104(3):739–44.PubMedPubMed CentralView ArticleGoogle Scholar
- Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C, et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28(12):1647–9.PubMedPubMed CentralView ArticleGoogle Scholar
- Gorge O, Bennett EA, Massilani D, Daligault J, Pruvost M, Geigl EM, Grange T. Analysis of ancient DNA in microbial ecology. Methods Mol Biol. 2016;1399:289–315.PubMedView ArticleGoogle Scholar
- Rueckert A, Morgan HW. Removal of contaminating DNA from polymerase chain reaction using ethidium monoazide. J Microbiol Methods. 2007;68(3):596–600.PubMedView ArticleGoogle Scholar
- Lopez-Barragan MJ, Quinones M, Cui K, Lemieux J, Zhao K, Su XZ. Effect of PCR extension temperature on high-throughput sequencing. Mol Biochem Parasitol. 2011;176(1):64–7.PubMedView ArticleGoogle Scholar
- Gnirke A, Melnikov A, Maguire J, Rogov P, LeProust EM, Brockman W, Fennell T, Giannoukos G, Fisher S, Russ C, et al. Solution hybrid selection with ultra-long oligonucleotides for massively parallel targeted sequencing. Nat Biotechnol. 2009;27(2):182–9.PubMedPubMed CentralView ArticleGoogle Scholar
- Renaud G, Stenzel U, Kelso J. leeHom: adaptor trimming and merging for Illumina sequencing reads. Nucleic Acids Res. 2014;42(18):e141.PubMedPubMed CentralView ArticleGoogle Scholar
- Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol. 2010;59(3):307–21.PubMedView ArticleGoogle Scholar
- Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30(9):1312–3.PubMedPubMed CentralView ArticleGoogle Scholar
- Drummond AJ, Rambaut A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007;7:214.PubMedPubMed CentralView ArticleGoogle Scholar
- Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9(8):772.PubMedPubMed CentralView ArticleGoogle Scholar
- Soares P, Ermini L, Thomson N, Mormina M, Rito T, Rohl A, Salas A, Oppenheimer S, Macaulay V, Richards MB. Correcting for purifying selection: an improved human mitochondrial molecular clock. Am J Hum Genet. 2009;84(6):740–59.PubMedPubMed CentralView ArticleGoogle Scholar
- Martínez-Navarro B, Ros-Montoya S, Espigares MP, Palmqvist P. Presence of the Asian origin Bovini, Hemibos sp. aff. Hemibos gracilis and Bison sp., at the early Pleistocene site of Venta Micena (Orce, Spain). Quatern Int. 2011;243(1):54–60.View ArticleGoogle Scholar
- Qui Z, Deng T, Wang B. Early Pleistocene mammalian fauna from Longdan, Dongxian, Gansu. China Paleontologia Sinica. 2004;191(C27):1–198.Google Scholar
- Airvaux J. Le site des Plumettes à Lussac-les-Châteaux (Vienne). In: Préhistoire de Poitou-Charentes. vol. Actes du 111é congrès des Sociétés Savantes. Paris: CTHS; 1987. pp. 193–200.Google Scholar
- Baryshnikov G. Late Pleistocene brown bear (Ursus arctos) from the Caucasus. Russ J Theriol. 2010;9(1):9–17.Google Scholar
- Baryshnikov G, Hoffecker JF, Burgess RL. Palaeontology and Zooarchaeology of Mezmaiskaya Cave (Northwestern Caucasus, Russia). J Arch Sci. 1996;23(3):313–35.View ArticleGoogle Scholar
- Baryshnikov GF. Pleistocene Felidae (Mammalia, Carnivora) from Paleolithic site in Kudaro caves in the Caucasus. Proc Zool Inst RAS. 2011;315(3):197–226.Google Scholar
- Beauval C, Morin E. Les repaires d’hyènes du Lussacois (Lussac-les-Châteaux,Vienne, France). Apport des sites des Plumettes et des Rochers de Villeneuve. In: Préhistoire entre Vienne et Charente Hommes et sociétés du Paléolithique. Buisson-Catil J, Primault J, editors. Chauvigny: Association des Publications Chauvinoises; 2010. pp. 175–90.Google Scholar
- Bitard B, Madelaine S. La grotte de La Berbie (Castels – Dordogne). Spéléo. 1994;AN8:31–40.Google Scholar
- Boeskorov GG, Bakulina NT, Davydov SP, Shchelchkova MV, Solomonov NG. Study of Pollen and Spores from the Stomach of a Fossil Woolly Rhinoceros Found in the Lower Reaches of the Kolyma River. Doklady Biol Sci. 2011;436:23–5.View ArticleGoogle Scholar
- Castel J-C, Boudadi-Maligne M, Ducasse S, Renard C, Chauvière F-X, Kuntz D, Mallye J-B. Animal exploitation strategies in Eastern Aquitaine (France) during the Last Glacial Maximum. In: Foulds FWF, Drinkall HC, Perri AR, Clinnick DTG, Walker JWP, editors. Wild Things: Recent Advances in Palaeolithic and Mesolithic Research. Durham Symposium, March 24–25, 2012. Oxford: Oxbow; 2014. p. 160–74.Google Scholar
- Castel J-C, Coumont M-P, Boudadi-Maligne M, Prucca A. Rôle et origine des grands carnivores dans les accumulations naturelles. Le cas des loups (Canis lupus) de l’Igue du Gral (Sauliac-sur-Célé, Lot, France). Revue de Paléobiologie Genève. 2010;29(2):411–25.Google Scholar
- Castel J-C, Coumont M-P, Brugal J-P, Laroulandie V, Camus H, Chauvière F-X, Cochard D, Guadelli J-L, Kuntz D, Martin H, et al. La fin du Paléolithique supérieur en Quercy: l’apport de l’Igue du Gral (Sauliac-sur-Célé, Lot). In: Jaubert J, Ortega I, Bordes J-G, editors. Les sociétés du Paléolithique dans un Grand Sud-Ouest: nouveaux gisements, nouveaux résultats, nouvelles methods. Journées de la Société Préhistorique Française, Bordeaux, 24–25 Novembre 2006, vol. Mémoire XLVII. Nanterre: Société Préhistorique Française; 2008. p. 335–53.Google Scholar
- Coumont M-P, Brugal J-P, Castel J-C, Costamagno S, Jaubert J, Mourre V. Les avens-pièges à faible indice de fréquentations humaines: caractérisation paléoécologique, taphonomique et anthropologique. In: Modalité d’occupation et exploitation des milieux au Paléolithique dans le Sud-Ouest de la France: l’exemple du Quercy Actes de la session C67, XVème Congrès mondial de l’UISPP, Lisbon, September 2006. Jarry M, Brugal J-P, Ferrier C, editors. vol. n° 4: Paleo, supplément; 2013. pp. 181–96.Google Scholar
- Discamps E. Hommes et hyènes face aux recompositions des communautés d’ongulés (MIS 5-3): Eléments pour un cadre paléoécologique des sociétés du Paléolithique moyen et supérieur ancien d’Europe de l’Ouest. PhD Dissertation. Bordeaux: University of Bordeaux; 2011.Google Scholar
- Gamberi Almendra de Carvalho L, Argant A, Argant J, Barth P, Boudadi-Maligne M, Boulbes N, Brugal J-P, Caramelli D, Condemi S, Crégut E, et al. L’aven de l’Arquet - Barjac (30) : étude d’un aven piège. Ardèche Archéologie. 2011;28:3–10.Google Scholar
- Golovanova LV, Hoffecker JF, Kharitonov VM, Romanova GP. Mezmaiskaya cave: a Neanderthal occupation in the Northern Caucasus. Curr Anthropol. 1999;40(1):77–86.View ArticleGoogle Scholar
- Lioubine VP. L’Acheuléen du Caucase, ERAUL 93. Liège: Université de Liège; 2002.Google Scholar
- Madelaine S. La Berbie (Castels). Bilan scientifique de la région aquitaine 2000. Paris: Ministère de la Culture et de la Communication, Sous-Direction de l’Archéologie; 2001. p. 24.Google Scholar
- Merceron G, Madelaine S. Molar microwear pattern and palaeoecology of ungulates from La Berbie (Dordogne, France): environmental context of the last Neanderthal populations. Boreas. 2006;35:272–8.View ArticleGoogle Scholar
- Pétrequin P. Les sites littoraux néolithiques de Clairvaux-Les-Lacs et de Chalain (Jura) III Chalain station 3, 3200–2900 av. J.C. Paris: Maison des Sciences de l’Homme; 1997.Google Scholar
- Pétrequin P, Arbogast R-M, Bourquin-Mignot C, Lavier C, Viellet A. Demographic growth, environmental changes and technical adaptations: responses of an agricultural community from the 32nd to the 30th centuries BC. World Archaeol. 1998;30(2):181–92.PubMedView ArticleGoogle Scholar
- Putelat O. Les relations homme-animal dans le monde des vivants et des morts, Étude archéozoologique des établissements et des regroupements funéraires ruraux de l’Arc jurassien et de la Plaine d’Alsace, de la fin de l’Antiquité tardive au premier Moyen Âge. PhD Dissertation. Paris: University of Paris 1 Panthéon-Sorbonne; 2015.Google Scholar
- Skinner AR, Blackwell BAB, Martin S, Ortega A, Blickstein JIB, Golovanova LV, Doronichev VB. ESR dating at Mezmaiskaya Cave, Russia. Appl Radiat Isot. 2005;62(2):219–24.PubMedView ArticleGoogle Scholar
- Tournepiche JF. Les grands mammifères pléistocènes de Poitou-Charente. Paléo. 1996;8:109–41.View ArticleGoogle Scholar