Introduction

The earliest members of the genus Bison appeared at the beginning of the Pleistocene in India and China. Bison then spread from Asia to Europe and America. Archaeological evidence shows that two bison species existed in Europe from the Middle Pleistocene on, the steppe bison and the woodland bison. The steppe bison, Bison priscus (Bojanus, 1827), is very common in Pleistocene deposits but became extinct at the end of the last Ice Age, about 10,000 years ago. It was a formidable animal with long horns and robust legs that may have stood more than 2 meters at the withers and reached a total length of more than 2.7 meters. It occupied cool, steppe-like grasslands. Bison priscus had a wide geographic distribution, called the great Pleistocene bison belt, which extended from England to America and from the archipelago of Nova Zembla in the Arctic Ocean in the north of Russia down to Spain, Caucasia and Mexico. It is generally reckoned that it is the steppe bison that was portrayed by paleolithic artists in painted caves like the Altamira cave in northern Spain, the Chauvet, Lascaux and Trois-Frères caves in France. The woodland bison, Bison schoetensacki (Freudenberg, 1910), which has left more scarce fossil remains than the steppe bison, appeared in the early Middle Pleistocene. Its body size and horns were smaller than those of the steppe bison and it was probably a forest-dweller. It became extinct at the end of the Pleistocene at the latest [1].

Only two bison species survive today: the American bison and the European bison. The American bison, in North America, includes two subspecies, the plains bison, Bison bison bison, and the wood bison, Bison bison athabascae. The European bison, Bison bonasus, is found in Europe and the Caucasus, where it has been reintroduced after its extinction as a wild species. The phylogenetic relationships between the extant and the extinct species are a matter of debate. A population of Bison priscus is believed to have first entered eastern Beringia from Asia during the middle Pleistocene, around 300 to 130 ky ago, and then moved southward into central North America approximately 130 to 75 ky ago [2]. The analysis of mitochondrial sequences of bison specimens from Siberia, China and North America (dated between present time and about 70 ky B.P. [before present]) has suggested that late Pleistocene bison from the Ural mountains to northern China are descendants of one or more dispersals from North America [2]. This study has also indicated gene flow between bison populations in Beringia and central North America between 60 and 25 ky B.P., before ice sheet formation during the Last Glacial Maximum (22 to 18 ky B.P.) hampered north-south faunal exchange [2]. The retreat of ice sheets around 14 ky B.P. created an ice-free corridor through which dispersal between Beringia and North America could occur. However, genetic exchanges between bison in Beringia and in central North America were soon after limited due to the establishment of spruce forest in Alberta and of peatland across western and northwestern Canada. As a result, all modern Bison bison belong to a clade distinct from Beringian bison [2]. The European bison Bison bonasus has a complex descent: its nuclear genome is closely related to that of Bison bison, with which it can produce completely fertile hybrid offspring, but its mitochondrial genome is more similar to the genomes of species belonging to the Bos genus than to the genome of Bison bison [3]. These observations could be explained by lineage sorting implying an early divergence of the American and European bison mitochondrial sequences, or by a scenario according to which the nuclear DNA of an ancient Eurasian Bos population would have been changed by the systematic introgression of Bison bulls (either Bison priscus, Bison bison migrating back to Europe, or Bison schoetensacki). This ancient hybridization would have been facilitated by the social structure of herd species, in which most males are excluded from reproduction by the dominating bulls, and would have created a new species with a nuclear genome of bison type, but with a Bos-like mitochondrial DNA from the maternal ancestor. Interestingly, this scenario would explain the sudden paleontological appearance of Bison bonasus at the Preboreal stage (10.3 to 9 ky B.P.) [3].

In an effort to investigate the evolution and demographic history of Pleistocene bison, partial sequences of Bison priscus mitochondrial genome have already been obtained [2,4]. The mitochondrial control region, also called the D-loop region, has been sequenced in ~ 300 Bison priscus fossils collected in Siberia, China and North America [2,4]. However, these partial sequences are a few hundred bp long, representing less than 5% of the expected 17-kb mitochondrial genome. The purpose of the present study was to fill the gap by providing a complete mitochondrial genome sequence for the extinct steppe bison.

This work presents the analysis of several bone samples that were collected from the Trois-Frères cave (Ariège, France). This cave (43° 1′ N, 1° 12′ E) was named after the three sons of Comte Henri Bégouën who discovered its entrance in July 1914. It contains numerous animal representations made famous by the publication of the Abbé Henri Breuil [5], including the drawing and engraving of at least 170 bison. The importance of the bison for the Magdalenian artists who decorated the cave is obvious not only from the fact that its representations outnumber those of any other species, but also from the drawing of a therianthropic figure with a bison upper part. Besides, the cave of the Trois-Frères is part of a single cave-complex with the Tuc d'Audoubert, which contains two modeled clay bison, a unique masterpiece of the Paleolithic period [6]. DNA extracted from one rib bone sample was sequenced using Illumina technology which, in combination with conventional PCR and Sanger sequencing, allowed the assembly of a complete sequence of Bison priscus mitochondrial genome.

Material and Methods

Results and Discussion

The Trois-Frères cave consists of a series of chambers that spread over 800 meters, from the Enlène cave to the Tuc d’Audoubert cave, from which it is currently disconnected (Fig 1A and 1B). In this cave system the Salle du Grand Éboulis (Chamber of the Large Scree) displays unique features. It functioned as a natural trap for Pleistocene animals, and was sealed before the Holocene [21]. The bottom of the scree contains remains of the extinct cave bear (Ursus spelaeus), one of which was dated to 36,600–34,800 calBP [21]. In the other layers remains of the extinct steppe bison predominate [22]. They correspond to 19,400 to 17,800-year-old specimens [21].

We analyzed four bone samples that were collected in the Salle du Grand Éboulis (SGE) for the purpose of the present study (Fig 1C–1E). Three samples stood on the ground surface, in the vicinity of the pathway delineated since the cave discovery (SGE1, SGE3), or at the bottom of the excavation carried out in the 30’s of the previous century (SGE5). The other sample (SGE2, a rib fragment) stood on a rock, three meters away from the pathway used by modern visitors.

We first analyzed the DNA content of the samples by PCR using primer pair 1 (S1 Table). The forward and reverse primers display a perfect match with the Bos taurus and Bos primigenius reference mitochondrial genome sequences and 0 to 1 mismatch with the Bison bison and Bison bonasus reference mitochondrial genomes. The mitochondrial DNA sequence flanked by these primers (63-bp) is identical between the Bos taurus and Bos primigenius reference genomes. Bison bison and Bison bonasus sequences both display 5 differences (including one transversion) with the two Bos genomes, and differ among each other by two transitions. All samples, but none of the negative controls (mock extracts and PCR blanks) yielded amplification products of the predicted size that were characterized through the sequencing of cloned PCR fragments (SGE1: 1 PCR, 13 clones; SGE2 extract 1: 4 PCR, 47 clones; SGE2 extract 2: 2 PCR, 39 clones; SGE3: 1PCR, 12 clones; SGE5: 1PCR, 14 clones). SGE1 and SGE3 DNA extracts yielded sequences that all displayed a perfect match with Bos (Bos taurus, Bos primigenius) reference mitochondrial genomes, and SGE5 returned a mixture of sequences that perfectly matched either Bos or Bison bison genomes (Fig 2). The PCR data may indicate that SGE1 and SGE3 samples correspond to Bos bone fragments, but the observation that SGE5 contains Bos as well as Bison DNA challenges this interpretation and rather suggests that all samples collected from the ground surface had been contaminated by modern cattle DNA brought by modern human intrusions since the cave discovery. At any rate, none of these three samples could be considered as a relevant candidate for a thorough analysis of bison DNA. The SGE2 bone sample turned out to be much more suitable for this purpose. The two SGE2 DNA extracts yielded sequences that better aligned to Bison bison than to Bos taurus and Bos primigenius mitochondrial genomes (Fig 2 and S2 Table). As a whole, the majority (79%) of these sequences displayed a perfect match with the Bison bison reference mitochondrial genome (5 mismatches with Bos genomes), the other ones bearing 1–3 and 6–8 differences with the Bison bison and Bos genomes, respectively.

To further guarantee the relevance of the SGE2 bone sample to study bison DNA, we performed an additional series of PCR experiments using primers designed from Bison (primer pairs 2 and 7–10) or Bos (primer pairs 3–6) reference mitochondrial genomes. As shown in S2 Table, all these experiments yielded DNA sequences that aligned with higher BLAST scores to Bison than to Bos mitochondrial genomes. Moreover, using primers targeting the portion of the Bison priscus mitochondrial genome made available by previous studies [2], namely the D-loop region, we obtained DNA sequences that clearly exhibited higher alignment scores with this species than with any other one. Thus, the bulk of evidence supported the notion that SGE2 corresponds to a Bison bone sample. We therefore decided to further analyze this rib fragment.

As shown in Fig 3, SGE2 was radiocarbon dated to 15,880 ± 70 B.P. (19,390–18,940 calBP). It is therefore representative of the age previously obtained for Bison priscus bones collected in the Salle du Grand Éboulis [21]. In order to obtain genomic data, we produced a library of DNA fragments from the SGE2 sample and performed Illumina sequencing.

A total of 1,033,460,536 Illumina single-pass reads, each at least 20 nucleotides in length, were produced. To confirm the origin of the SGE2 sample and further rule out contamination by modern DNA, the reads were aligned to the reference mitochondrial genomes of Bison bison, Bison bonasus, Bos primigenius and Bos taurus. The alignments were performed under stringent conditions: only the reads with a minimal length of 30 bp matching perfectly the target sequence were retained. Under these conditions, 1,825 reads mapped to the Bison bison mitochondrial genome. By contrast, only 324, 347 and 331 reads mapped to the mitochondrial genomes of Bison bonasus, Bos primigenius and Bos taurus, respectively. Moreover, the large majority of the reads mapping to the Bison bonasus and Bos genomes also mapped to the Bison bison genome: this was the case for 296 reads out of 324 for the Bison bonasus genome, and for 305 reads out of 347 and for 290 reads out of 331, respectively, for the Bos primigenius and Bos taurus genomes. These results strongly suggest that the reads mapping to the Bison bonasus and Bos genomes correspond to ancient bison DNA and map to regions that are identical in the genome of Bison bison and in the genomes of the other species. The alignment analysis thus corroborated the conclusion, previously drawn from the PCR results, that the SGE2 extract contains Bison rather than Bos DNA.

To evaluate the coverage of the nuclear genome by the Illumina reads, we analyzed the 18S ribosomal RNA gene, which displays approximately 200 copies in mammalian genomes [23]. Out of ~ 10⁹ reads, 110 aligned to the 18S rDNA sequence of Bos taurus. These results correspond to an average read depth of 2.2 for 200 copies and accordingly to ~ 0.01 for single copy sequences. Mapping to the complete nuclear genome of Bos taurus yielded 741,848 unique reads with a mapping quality higher than 20. The sum of these reads lengths amounts to 26,117,728 bp, which represents a 0.01-fold coverage of the 2,7 Gb Bos taurus nuclear genome, in agreement with our previous estimate. These analyses confirmed the poor coverage of the nuclear genome, which precluded its further analysis.

To gain some further insight into the content of the Illumina library, one million reads were analyzed by BLAST, which yielded 99,684 significant hits. The large majority of these hits (89%) corresponded to Bacteria while Eukaryota represented only 9% of the total. In the bacterial metagenome, the predominant classes were Actinobacteria (32% of all bacterial hits), Betaproteobacteria (20%), Gammaproteobacteria (18%) and Alphaproteobacteria (16%). Similar distributions of microbial diversity have been described for ancient DNA extracts (e.g. [24]). Among the Eukaryota, the family most represented was the Poaceae (10% of all eukaryotic hits), followed by the Bovidae (6%) and the Noelaerhabdaceae (3%). The sum of the lengths of the 540 Illumina reads corresponding to the Bovidae hits represents a total of 26,856 nucleotides for one million Illumina reads analyzed, predicting 27,754,600 nucleotides for Bovidae sequences in the whole set of reads.

To analyze the Bison mitochondrial genome sequence, we first retrieved from the library all the reads that aligned to the mitochondrial genome of Bison bison (which, as indicated above, turned out to be more relevant for this purpose than Bison bonasus) or to the mitochondrial control region of Bison priscus isolate BS212 (GenBank accession number AY748539.1, [2]). This alignment was performed under less stringent conditions than the previous ones and allowed for mismatches (see details in Materials and Methods). A provisional consensus sequence was derived from the selected reads and was used in turn to retrieve all reads aligning to it. PCR experiments were carried out to analyze regions where less than two concordant reads were available. The final sequence, termed SGE2seq, was established using 81 consensus sequences of PCR fragments to complement the information derived from the Illumina reads. All the Illumina reads were aligned one last time to SGE2seq, which yielded 3,851 unique reads.

SGE2seq consists of a circular genome, for which we obtained a 10.4-fold coverage, taking into account the 3,851 unique Illumina reads and the 81 PCR sequences. Fig 4A shows the length distribution of the 3,851 unique Illumina reads and Fig 4B, the coverage obtained for each position.

Only 1,443 mismatches out of 166,849 aligned bases are observed between the 3,851 unique Illumina reads and the consensus sequence SGE2seq, which corresponds to 99.1% identity. As shown in S1 Fig, the G to A substitution rate increases at 3’ ends, in agreement with what has already been reported for libraries of ancient DNA constructed with the same procedure [25]. An increased substitution rate at DNA ends is considered as a hallmark authenticating sequences generated from ancient DNA fragments which exhibit inflated cytosine deamination rates at 5'-overhangs [26], responsible for G to A transitions at the 3’ end of the opposite, repaired DNA strand.

SGE2seq corresponds to a mitochondrial genome sequence of 16,318 bp. This length is similar to the lengths of Bison bison mitochondrial genomes, which range from 16,318 to 16,323 bp [27]. The genome consists of 13 protein-coding genes, 22 tDNA genes, 2 rDNA genes, and the D-loop region (S3 Table).

Since part of the D-loop region has been sequenced for a number of Bison priscus fossils [2,4], we analyzed the phylogenetic relationships between the orthologous portion of SGE2seq and of several Bison priscus fossils and Bison bison specimens. We also considered the D-loop region of several reference genomes of Bison and Bos species, including the Bos grunniens (yak) reference mitochondrial genome (GenBank accession number NC_006380). The phylogenetic reconstruction, carried out using the Maximum Likelihood method, positions SGE2seq among Bison priscus fossils sequences, outside the strongly supported clade of modern Bison bison sequences (S2 Fig). This result was observed when taking into account three specimens from each of the four Bison groups defined in [2], and similar results were obtained when including all published D-loop sequences of Bison priscus and Bison bison (data not shown). The position of SGE2seq D-loop region in the phylogenetic tree further corroborates the assumption that SGE2seq corresponds to the mitochondrial genome of a Bison priscus specimen.

Several lines of evidence indicate that we obtained a reliable mitochondrial genome sequence for a Pleistocene Bison priscus specimen. First, the DNA material was retrieved from a 19,390 to 18,940-year-old bone sample (Fig 3). Second, considering the mitochondrial genome portion (600 bp of the D-loop) previously characterized for the extinct steppe bison [2,4], phylogenetic analysis demonstrated that the sequence assembled from our sample positions in the cluster delineated by Bison priscus specimens (S2 Fig). Third, we observed the damage-induced misincorporation pattern expected for ancient DNA (S1 Fig). Fourth, although of limited magnitude, the average unique read depth (10.4-fold) is large enough to overcome the error rate (< 1%, as judged by comparing the reads with the consensus sequence), and positions poorly covered were checked by PCR and sequencing of multiple cloned amplicons. Fifth, we obtained a circular genome that contains the expected number of genes (rRNAs, tRNAS, protein-coding genes) for a mammalian mitochondrial genome, with full-length coding sequences (without internal stop codons) for all 13 protein-coding genes. Thus, it is safe to conclude that we provide an authentic and complete Bison priscus mitochondrial genome sequence.

When single-nucleotide polymorphisms (SNPs) and insertion-deletion events are considered, there are 106–117 differences between SGE2seq and the 33 published Bison bison mitochondrial genomes (Bison bison genomes differ among themselves at 24 positions at most). SGE2seq shows 114 differences (listed in S4 Table) compared to the reference Bison bison mitochondrial genome, including two indels: a 1-bp insertion located in the D-loop region and a 2-bp deletion at the end of the tRNA serine gene (this shorter tRNA serine allele is also present in some Bison bison mitochondrial genomes [27]). The other differences comprise 105 transitions and 6 transversions. A similar, strong transitional bias for mitochondrial genomes has also been reported for Bos species [14]. Twenty-five differences (22%) are located in the D-loop region, which therefore exhibits a mutation rate four times higher than the rest of the genome. Finally, 73 substitutions (65%) are located in protein-coding genes, with 16, 3 and 54 substitutions occurring in the first, second and third codon positions, respectively. These substitutions correspond to 11 amino acid differences. Among the 114 differences, 93 correspond to polymorphisms that are specific to Bison priscus and are not found in any published Bison bison genome (S5 Table). These 93 specific polymorphisms are mostly located in protein-coding genes (69%) and correspond to 7 amino acids differences.

A phylogenetic analysis was performed using all the available complete mitochondrial genomes of Bison bison, Bison bonasus and Bos grunniens along with the reference mitochondrial genomes of Bos primigenius, Bos taurus and Bubalus bubalis (swamp buffalo). Phylogenetic trees were constructed from this dataset with the Maximum Likelihood (ML), Minimum Evolution (ME) and Neighbor-Joining (NJ) methods. As shown in Fig 5 for the ML method, SGE2seq is basal to the modern Bison bison clade, and the set including SGE2seq and all Bison bison sequences forms a distinct and well-supported (100% bootstrap support, 1000 replicates) clade. Similar results were obtained with the ME and NJ methods, which conclusively supports the contention that SGE2seq corresponds to a bona fide mitochondrial sequence of Bison priscus.

In summary, we have presented here the analysis of a 19,000-year-old bone specimen that made it possible to obtain the first complete mitochondrial genome for the extinct steppe bison. Phylogenetic reconstruction unequivocally points to the basal position of the Bison priscus mitochondrial genome as compared to the clade delineated by extant Bison bison genomes. The sequence provided here may be a useful tool for the genomic analysis of other ancient Bovinae specimens. Our study also shows that the Trois-Frères cave system, already famous for Magdalenian art, evidence of human habitation and abundant animal remains, can now be considered for paleogenetic and paleogenomic analyses.

Supporting Information

Acknowledgments

We thank A. Martel (CEA Saclay) for bioinformatics support.