BACKGROUND
during the middle and upper pleistocene, the large bovidae in europe and in northern asia chiefly included the aurochs, bos primigenius , and two bison species, the steppe bison and the woodland bison.
the steppe bison, bison priscus , was very common and exhibited a wide geographic distribution stretching from western europe, through central asia and beringia, and into north america.
it was a formidable animal with long horns and robust legs; it stood up to two meters at the withers and reached a total length of almost three meters.
it occupied cool, steppe-like grasslands.
bison priscus became extinct in europe at the end of the last ice age, about  <dig>  years ago  <cit> .
the woodland bison, bison schoetensacki  appeared in the early middle pleistocene.
its size was almost as large as that of bison priscus but its leg bones and metapodials were slenderer than those of bison priscus  <cit> .
moreover the horns of bison schoetensacki were about 30% shorter than those of bison priscus and displayed a slightly different shape: the initial third part is oriented forward, the second third goes upward, and the final part slightly bends backward.
bison schoetensacki remains are often associated with forest biotopes, hence its vernacular name.
the geographic distribution of bison schoetensacki extended from western europe to the south of siberia, but unlike bison priscus, bison schoetensacki was absent from beringia and north america.
the species was first described in germany where it has been recorded in several lower middle pleistocene sites, notably in mosbach, mauer, süßenborn and voigtstedt, but also in roter berg and quedlinburg .
it has also been found in england , in moldova , in russia , and in north sea deposits.
in france, bison schoetensacki is recorded in several lower middle pleistocene sites  and in late middle pleistocene deposits in châtillon-saint-jean  <cit> .
it seems to be still present during the upper pleistocène at chasse-sur-rhône, poleymieux and in the herm and siréjol caves.
however the persistence of bison schoetensacki during the upper pleistocene is under debate, and siréjol bone remains have been alternatively ascribed to bison schoetensacki  <cit>  and to a bison priscus novel subspecies  <cit> .
bison schoetensacki remains from siréjol display a cranial frontal no larger than 214 mm, short horns of the shape described above with maximal vertical and anteroposterior diameter at their base of  <dig> and 65 mm, respectively  <cit> .fig.
<dig> cave sites and samples.
black dots indicate sites where bison schoetensacki remains have been reported in previous studies.
red characters indicate the two cave sites that yielded the samples analyzed in the present study.
the coprolite from the grotte-aux-ours cave .
the bison schoetensacki cannon bone sample from the siréjol cave
two bison species exist today: the american bison and the european bison.
the genetic history of the american bison  has been extensively investigated through the analysis of mitochondrial dna sequences .
the recent reconstruction of the complete mitochondrial genome of a  <dig> -year-old bison priscus specimen confirmed that the closest extant mitochondrial genome sequences for bison priscus are those of bison bison  <cit> .
these data, in combination with nuclear single nucleotide polymorphism  genotyping  <cit> , have established that bison priscus and bison bison are sister groups.
the european bison, bison bonasus, is found in europe and the caucasus, where it has been reintroduced after its extinction as a wild species in the early 20th century.
the bison bonasus genomes reflect a complex descent.
indeed, the bison bonasus nuclear genome is closely related to that of bison bison , in agreement with morphological evidence and the fact that the two bison species can produce completely fertile hybrid offspring.
however, the mitochondrial genomes of bison bonasus specimens are more similar to aurochs and cattle genomes than to bison bison genomes  <cit> .
these observations could be explained either by incomplete lineage sorting of the mitochondrial genome or by a scenario according to which the nuclear dna of an ancient population probably related to the extinct aurochs  would have been changed by the systematic introgression of bison bulls  <cit> .
the phylogenetic relationships between bison bonasus, bison priscus and bison schoetensacki are subject to debate.
some authors suggest that bison schoetensacki could be the ancestor of bison bonasus  <cit> , but others consider that bison bonasus is derived from an unknown form of bison priscus  <cit> .
the lack of genomic data for bison schoetensacki has so far prevented any conclusive view on this point.
as part of an effort to obtain genetic data on animals consumed by the extinct cave hyena , we initiated the analysis of cave hyena coprolites  <cit> .
in the present study, we found a coprolite containing mitochondrial bovinae dna that was different from bos primigenius, bison priscus and bison bonasus mitochondrial sequences.
because the cave site  that yielded this coprolite is located close to the siréjol cave, we also analyzed dna extracted from a siréjol bison.
schoetensacki bone sample.
we report here on the mitochondrial dna sequences  obtained from the coprolite and the bone remain, on their phylogenetic positions, and on a set of snps of the corresponding nuclear genome sequence.
methods
archeological samples
the cave hyena coprolite was retrieved from the grotte-aux-ours cave .
the cave was discovered in  <dig> and owes its name to the abundant evidence for long-term occupation by the extinct cave bear , including cave bear bones, hibernation nests, claw marks on the walls, and footprints on the paleosurface.
occupation by the extinct cave hyena  is documented by the presence of coprolites and a cave hyena skull.
the bison schoetensacki bone fragment originated from the siréjol cave .
this cave, discovered in  <dig>  is located 15 km away from the grotte-aux-ours cave.
scientific excavations carried out from  <dig> to  <dig> in the siréjol cave yielded bison schoetensacki bone remains  <cit>  that are currently stored in the centre de conservation et d’Étude des collections du musée des confluences .
for the purpose of the present study, we analyzed a bison schoetensacki cannon bone fragment .
contamination issues
to guarantee that genuine ancient dna was recovered from the samples, dna extraction was carried out from material derived from the inside of the samples, and the coprolite and bone dna extracts were obtained on different days.
additional precautions to avoid contamination have been described elsewhere  <cit> .
negative controls included mock extracts and pcr blanks , which always failed to yield any amplification product.
dna extraction
the instruments used to cut and pulverize the samples were cleaned in an ultrasonic water-tank, incubated 1 h at 40 °c in  <dig> % sds, rinsed  <dig> times in ultrapure dnase/rnase free distilled water, and dried using absolute ethanol.
the coprolite and the bone sample were cut using a carbon-free circular saw.
for dna extraction we used  <dig>  g of biological material retrieved from the center of the coprolite or the bone medulla using a sterile single-use scalpel.
the dissected material was pulverized using a mikro-dismembrator s  set to a shaking frequency of  <dig>  rpm for 30 s, transferred into 6 ml of dna extraction buffer ,  <dig> % sds, 65 mm dtt,  <dig>  mg/ml proteinase k), and incubated 16 h at 42 °c under constant agitation.
after centrifugation, the supernatant was recovered, extracted once with one volume of phenol, once with a phenol-chloroform  mixture, and once with chloroform.
the aqueous phase was dialyzed and concentrated using a 2 ml, centrifugal ultrafilter with a 30 kda cutoff , the column was washed  <dig> times with distilled water, and the dna extract was recovered in a volume of  <dig> to 120 μl.
it was purified further using a silica-membrane-based method designed for small dna fragments .
the final extract had a volume of 100 μl.
pcr and sanger sequencing dna analysis
pcr reactions were performed in a 50-μl reaction volume containing  <dig> – <dig>  μl of mock or ancient dna extracts,  <dig> nm of sense and antisense primers, 200 μm dntp,  <dig>  mm mgcl <dig>  5 μl of geneamp 10x pcr buffer ii, and  <dig>  u of amplitaq gold dna polymerase .
the primers used for experiments designed to complete the bovine mitochondrial genome sequence assembled from the grotte-aux-ours coprolite are listed in additional file 1: table s <dig>  the primers used to analyze mitochondrial genome fragments of bison schoetensacki bone sample from the siréjol cave are listed in additional file 1: table s <dig>  primer pair  <dig> from additional file 1: table s <dig> was used for initial screening of the samples for bovinae dna content.
pcr procedures consisted of an enzyme activation step , followed by a single round of  <dig> pcr cycles , 20 s; 70 °c, 1 min) performed in a veriti thermal cycler .
the full reaction volume was loaded onto an 8% polyacrylamide gel stained with sybr green i .
only the coprolite and bone extracts yielded amplification products of the predicted size.
pcr amplicons were eluted from the gel and cloned using the topo ta cloning kit .
sequencing analysis of cloned dna fragments was performed on abi  <dig> or 3130xl automatic dna sequencers using m <dig> forward primer and bigdye  <dig>  terminator chemistry .
all mitochondrial genome fragments analyzed by pcr were characterized through 2– <dig> pcr replicates and a consensus sequence for each replicate was derived from several clones .
for the bovine mitochondrial genome sequence reconstructed from the grotte-aux-ours sample, each individual consensus was added to the illumina sequence data .
for the siréjol bone sample, several individual consensuses of the same fragments were combined to deduce the most relevant base of the  <dig> positions of the bison schoetensacki mitochondrial genome sequence investigated by pcr.
generation and sequencing of a library of dna fragments
a library of dna fragments suitable for high-throughput sequencing with the illumina procedure  <cit>  was generated from the coprolite dna extract using the illumina truseq nano dna lt sample kit fc-121- <dig>  we followed the manufacturer’s recommendations , except for modifications that were introduced for the purpose of analyzing ancient, i.e.
highly fragmented dna fragments.
thus, the initial step  of the library construction process was omitted, and reaction products were purified using qiagen columns or polyacrylamide gel electrophoresis instead of magnetic beads.
the process of library construction consisted of  <dig> steps.
first, 6 μl of coprolite dna was 5’end-phosphorylated and blunt-ended using a mixture containing t <dig> polynucleotide kinase, t <dig> and klenow dna polymerases, and the reaction product was purified on a qiagen  <dig> column.
second, a 3’ adenine residue was added to the blunt-ended dna fragments using klenow 3’ to 5’ exo- polymerase, and the enzyme was heat-inactivated at 70 °c.
third, illumina adapters with an overhanging thymine were ligated to the dna fragments; the reaction product was purified using a qiagen  <dig> column and recovered in a volume of 30 μl referred to as the dna library.
fourth, a 5-μl aliquot of the library was pcr-amplified using phusion dna polymerase .
fifth, the full pcr reaction volume was loaded on an 8% polyacrylamide gel stained with sybr green i, and the  <dig> to 220-bp long dna fragments  were cut off the gel and purified.
dna sequencing was performed at genoscope  on the illumina hiseq  <dig> platform  using hiseq v <dig> chemistry with a read length set to  <dig> nucleotides and analysis on the single read mode.
illumina sequencing data analysis
reads were trimmed for adapter sequences, n’s and low quality stretches on the 3’ end, using a software based on the fastx-toolkit package  and designed by genoscope.
after this step, sequences shorter than  <dig> nucleotides were discarded using an in-house python script, yielding a dataset of  <dig> , <dig> dna reads.
initial analysis of the illumina library content
to confirm the identity of the coprolite producer and gain insight into the animal’s diet, the  <dig> , <dig> illumina reads were aligned simultaneously to a set of  <dig> mitochondrial genomes, including the reference mitochondrial genomes for the extinct cave hyena  and for several potential prey species .
the alignment was performed using bwa version  <dig> .
<dig>  <cit> .
only the reads matching perfectly, without indel or mismatch, to a unique genome were taken into account.
assembling simultaneously the bovine and the cave hyena mitochondrial genomes
the sequence assembly was performed in several steps.
first the illumina reads, ranging in size from  <dig> to  <dig> nucleotides, were aligned simultaneously to the reference mitochondrial genomes of the extant european bison  and extinct cave hyena  using bwa version  <dig> .
<dig>  <cit>  with default options except for -o  and -n , which were set to  <dig> and  <dig> , respectively, allowing for a maximum edit distance that increases with read length, from one mismatch for 20– <dig> nucleotide-long reads to four mismatches for 88– <dig> nucleotide-long reads.
in this first round, the number of reads mapping to the bison bonasus genome  was of the same order of magnitude as the number of reads mapping to the crocuta crocuta genome .
of significance,  <dig>  reads mapped to both genomes, which justified our strategy to simultaneously assemble the two mitochondrial genomes so as to eliminate spurious alignments.
a set of  <dig>  unique reads mapping exclusively on the bison bonasus genome  with a mapping quality higher than  <dig> was selected.
the reads provided a median coverage of  <dig> of the bison bonasus reference sequence and left  <dig> positions without coverage.
a provisional consensus sequence termed bb1seq was then derived from the  <dig>  unique reads.
bb1seq differed from the bison bonasus reference genome at  <dig> positions.
in this study a consensus for a given position was based on at least two concordant unique sequences  with a minimal mapping quality of  <dig> for illumina reads.
the most frequent base was taken as the consensus.
the few cases where more than 25% of the bases from the aligned reads differed from the consensus were manually examined: the large majority of the alternative bases could be ascribed to damage-induced g to a substitutions at the 3' end of the illumina reads and the other cases were checked by sequence data obtained from pcr studies.
the cave hyena mitochondrial genome was assembled in parallel, using the same procedure.
a provisional consensus sequence termed cc1seq was derived from  <dig>  unique reads that aligned specifically to the crocuta crocuta reference mitochondrial genome.
the reads gave a median coverage of  <dig> for the crocuta crocuta reference sequence and left  <dig> positions uncovered.
cc1seq and the crocuta crocuta reference genome differed at  <dig> positions.
in a second step, all illumina reads were realigned simultaneously to the provisional consensus sequences bb1seq and cc1seq and new consensus sequences, termed bb2seq for the bovine genome and cc2seq for the cave hyena genome, were derived following the procedure previously described.
for the bovine sequence, pcr experiments were then carried out to fill in the gaps  and derive a robust sequence at positions where only a single read was available .
the  <dig> pcr primer pairs were designed using the bb2seq sequence .
a total of  <dig> consensus pcr sequences were used to complement bb2seq information and to derive a new consensus sequence, termed bb3seq.
all illumina reads were again simultaneously aligned to cc2seq and bb3seq.
after a last iteration, we were able to establish the final, complete consensus sequence, termed gaoseq_bovinae, for the bovine mitochondrial genome, based on  <dig>  unique illumina reads, and a partial sequence for the cave hyena mitochondrial genome, called gaoseq_crocuta, based on  <dig>  unique illumina reads.
annotation of the gaoseq_bovinae mitochondrial genome
the annotation of gaoseq_bovinae was derived from the annotation of the bison bonasus reference mitochondrial genome.
briefly, gaoseq_bovinae and the bison bonasus reference mitochondrial genome were aligned using the lagan software  <cit>  with the annotation data of the bison bonasus genome .
lagan produced a list of the aligned genomic features with their coordinates in the two genomes, which was manually checked.
the occurrence of a start and stop codon at the beginning and the end of each coding sequence, and the absence of internal stop codons were verified using an in-house python script.
phylogenetic analyses
all phylogenetic analyses were performed using mega  <cit> .
the phylogenetic relationships between gaoseq_bovinae and the complete mitochondrial genomes of specimens related to bison bonasus were inferred using the maximum likelihood , the minimum evolution and the neighbor-joining methods.
the ml method was based on the general time reversible model.
a discrete gamma distribution was used to model evolutionary rate differences among sites .
the analysis involved a total of  <dig>  positions.
a similar approach was used to analyze the dna sequences of the bison schoetensacki bone sample.
sixteen consensus sequences, each one derived from several clones of several pcr replicates, were aligned to the set of bison sequences previously used for the analysis of the complete mitochondrial genomes.
the aligned parts were subsequently concatenated and a phylogenetic tree was constructed using the ml method based on the tamura-nei model  <cit>  and a discrete gamma distribution to model evolutionary rate differences among sites .
analysis of the bovine nuclear sequences
in order to compare the bison nuclear sequences of the coprolite to the genomic data  made recently available for ancient specimens related to bison bonasus  <cit> , the library reads were mapped against the taurine cattle reference umd  <dig>   <cit>  using bwa, and the bovine snps were called using samtools and bcftools v <dig>   <cit> .
the calls were filtered using the same criteria as in  <cit>  .
the final number of snps called for the bovine sequences of the coprolite was  <dig> , which is comparable to the number of snps called for the a <dig> sample , but lower than the number of snps called for the other bison samples .
these data were completed by extracting the sequences of  <dig>  snps from the set of  <dig>  unplaced genomic scaffolds corresponding to the nuclear genome of bison bison bison isolate tamuid  <dig> .
multidimensional scaling analysis was performed using the r software v <dig> .
<dig>  <cit>  and in particular the cmdscale() function.
ams 14c dating
radiocarbon dating was attempted on both the hyena coprolite and the bison schoetensacki bone sample.
only the bone sample yielded ams 14c measurements.
no bone, hair, or macrofloral remains were present in the coprolite, and what constituted the matrix dissolved during weak-alkali  extraction.
consequently, no direct 14c measurements were possible on the coprolite.
part of the siréjol bone sample was broken into 3–4 mm fragments, which were washed in deionized  water with brief sonication to remove adhering sediment.
the bone fragments were subsequently decalcified with repeated changes of  <dig>  n hcl for 2 days at 4 °c and the acid-insoluble collagen washed to neutrality with di water.
the collagen was treated with  <dig>  m koh overnight at 4 °c to remove humic acids.
an aliquot of the alkali-extracted fraction was washed to neutrality with di water, acidified with  <dig>  n hcl, freeze dried and used for 14c dating of the koh-extracted, decalcified collagen fraction .
the remainder of the koh-extracted collagen was used for gelatin extraction, which comprised heating the alkali-extracted collagen to 90 °c in  <dig>  n hcl until dissolution occurred, in 45 min.
the gelatin solution was filtered through a  <dig>  μm durapore filter membrane and freeze dried.
approximately 5 mg of collagen and gelatin were weighed into quartz tubes with cuo and ag, the tubes evacuated to <  <dig> mtorr under ln pumping, sealed with a h2/o <dig> torch and the evacuated tubes combusted at 850 °c for 1 h. the co <dig> was isolated cryogenically, converted into graphite using the h2-fe method  <cit>  and used for 14c dating of the gelatin fraction .
both fractions were dated at the direct ams accelerator facility, seattle, washington usa.
the calibrated age of the sample, expressed in calendar years before present , was calculated using oxcal  <dig>  software  <cit>  and the intcal  <dig> calibration curve  <cit> .
RESULTS
analysis of the coprolite collected from the grotte-aux-ours cave
figure 1b-c shows the coprolite collected in the grotte-aux-ours cave and the bison schoetensacki bone fragment from the siréjol cave.
we obtained for the bone fragment radiocarbon dates of  <dig>  ±  <dig> and  <dig>  ±  <dig> bp using the koh-extracted decalcified collagen and the gelatin fraction of the sample, respectively.
we considered the average value  of these two measurements as the most reliable determination of the sample age, which corresponds to  <dig> – <dig>  cal yr bp .
previous analyses performed using conventional 14c dating of bulk bone material from the same cave sector yielded ages of  <dig>   bp and  <dig>   bp  <cit> .
our ams radiocarbon studies agree with this time frame and greatly increase the accuracy of dating analysis of the siréjol material.
dating analysis of the coprolite failed because the sample completely dissolved during the pretreatment for ams measurement.fig.
<dig> uncalibrated  and calibrated  age of the siréjol bison schoetensacki bone sample.
the calbp data correspond to  <dig> %  confidence interval of the sample age using intcal <dig> calibration curve
pcr analysis of the coprolite dna for the mitochondrial cytb fragments demonstrated the presence of cave hyena dna and of a bovinae sequence different from all of those recorded in genbank up to march  <dig>  the 82-bp bovinae cytb fragment displayed  <dig> to  <dig> mismatches with the bison bonasus, bison priscus, and bos primigenius mitochondrial reference sequences.
analysis of the siréjol bone fragment with the same bovinae primers yielded a sequence identical to the one obtained from the coprolite.
these pcr studies, and subsequent quantitative pcr analyses  also indicated far higher amounts of bovine mitochondrial dna in the coprolite than in the bone sample.
we therefore decided to perform high-throughput sequencing of the coprolite dna extract and to undertake a pcr approach for analyzing the bone sample.
a total of  <dig> , <dig> illumina single-pass reads, each at least  <dig> nucleotides in length, were produced from the coprolite dna library.
the reads were aligned simultaneously to a set of  <dig> mitochondrial genomes, including the reference mitochondrial genomes for the extinct cave hyena and several species likely to be part of a carnivorous diet.
only the reads matching perfectly, without indel or mismatch, to a unique genome were considered.
as shown in additional file 1: table s <dig>  this analysis confirmed that the extinct cave hyena was indeed the producer of the coprolite, since the number of reads  mapping specifically to the cave hyena genome was by far the largest.
a large number of reads  mapped to the european bison  genome.
the number of reads mapping specifically to other genomes then plummeted:  <dig> reads for bos primigenius,  <dig> reads for bison priscus,  <dig> reads for the chamois ,  <dig> reads for the musk ox  and less than  <dig> reads for the other genomes tested.
of importance, only one read mapped specifically to the reference mitochondrial genome of homo sapiens, attesting to the absence of contamination of the library by human dna.
these results thus indicated that the coprolite producer was a cave hyena that had ingested a bovine specimen whose closest known relative is bison bonasus.
assembling the bovine mitochondrial genome
initial analyses had revealed that the coprolite had been produced by a cave hyena that preyed on a bovine related to bison bonasus, hereafter referred to as "the gao bovinae".
we set out to reconstitute the mitochondrial genomes of both specimens by aligning the illumina reads to the reference mitochondrial genomes of bison bonasus and crocuta crocuta.
the reconstruction of the two genomes was carried out in parallel so as to eliminate spuriously aligned reads .
we assembled a  <dig> -bp bovine mitochondrial genome, termed gaoseq_bovinae, from  <dig>  unique illumina reads complemented with  <dig> consensus pcr sequences.
we obtained for gaoseq_bovinae a 32-fold median coverage, taking into account both the illumina reads and the consensus pcr sequences.
as shown in additional file 1: figure s <dig>  coverage by illumina reads was quite variable along the genome and strongly correlated with the sequence gc content, as already observed in other cases for illumina sequencing  <cit> .
the length distribution of the  <dig>  unique illumina reads is shown in additional file 1: figure s3a.
the median read length was  <dig> base pairs, as expected for ancient, highly fragmented dna.
we analyzed the differences between the illumina reads and the consensus sequence gaoseq_bovinae.
only  <dig>  mismatches out of  <dig>  aligned bases were observed between the  <dig>  unique illumina reads and gaoseq_bovinae, which corresponds to  <dig> % identity.
analysis of the positions of these differences with respect to the read ends showed that 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  <cit> .
an increased substitution rate at the read ends is considered as a hallmark for authenticating sequences generated from ancient dna fragments.
such fragments exhibit inflated cytosine deamination rates at 5'-overhangs  <cit> , responsible for g to a transitions at the 3’ end of the opposite, repaired dna strand.
gaoseq_bovinae is a complete circular genome that contains the expected number of genes  for a mammalian mitochondrial genome, with full-length coding sequences  for all  <dig> protein-coding genes .
with respect to the reference mitochondrial genome of bison bonasus, gaoseq_bovinae exhibits  <dig> differences, including  <dig> indels and  <dig> substitutions with a large majority  of transitions.
in comparison, the three available mitochondrial genomes of bison bonasus display less than  <dig> differences between each other.
while our manuscript was under review, two articles were published  <cit>  that describe several mitochondrial genomes from late pleistocene/holocene bison specimens related to bison bonasus.
we therefore performed a phylogenetic analysis including all these ancient mitochondrial genomes along with the mitochondrial genomes of bison bonasus specimens.
we also considered bison priscus specimens from these studies and a previous report from our laboratory  <cit> , as well as the reference mitochondrial genomes of bos primigenius, bison bison, bos grunniens  and bubalus bubalis  .
phylogenetic trees were constructed from this dataset with the maximum likelihood , minimum evolution and neighbor-joining methods that all yielded similar topologies.
as shown in fig.
<dig> for the ml method, gaoseq_bovinae positions in a well-supported clade  that includes all the ancient genomes  forming cladex in  <cit>  and clade bb <dig> in  <cit> .
no species from the fossil record has been proposed up to now for this clade which delineates a sister species for bison bonasus ancient , historical, and extant specimens from clade  <dig> .
mitochondrial genomes of clades  <dig> and  <dig> are phylogenetically closer to extant and ancient bos primigenius than to bison priscus and extant american bison.fig.
<dig> maximum likelihood phylogenetic tree of complete mitochondrial genomes of bison and bos species.
the tree with the highest log-likelihood is shown drawn to scale, with branch lengths established from the numbers of substitutions per site.
the percentages of trees in which the associated taxa clustered together are displayed next to the branches .
the genbank accession numbers of the sequences are listed in additional file 1: table s <dig>  specimen name and country of origin are indicated for ancient and historical samples.
the tree was rooted using the sequence of bubalus bubalis.
refer to text for explanation about clade  <dig> and clade 2
analysis of mitochondrial dna sequences derived from a bison schoetensacki bone sample
we explored the phylogenetic relationship between clade  <dig> and bison schoetensacki sequences by analyzing mitochondrial dna fragments of the siréjol bone sample.
for this we designed a set of  <dig> primer pairs allowing pcr amplification of short mitochondrial dna fragments  that were analyzed through the sequencing of cloned amplicons from several pcr replicates.
concatenation of the consensus sequences of the  <dig> pcr fragments yielded information on 609 bp of the mitochondrial genome.
the phylogenetic tree obtained using this concatenated sequence and orthologous fragments of bison mitochondrial genomes  strongly supports the position of the siréjol bison schoetensacki specimen within clade  <dig> and illustrates the similarity between the bison schoetensacki sequence and the sequences of arq <dig> and gaoseq_bovinae, from which it differs at  <dig> and  <dig> positions, respectively.
by contrast, the bison schoetensacki sequence displays  <dig> to  <dig> differences with the corresponding fragments of the mitochondrial genomes belonging to clade  <dig>  including the modern bison bonasus sequences, and  <dig> and  <dig> differences with the reference mitochondrial genomes of bison priscus and bison bison.
these observations support the conclusion that the bovine dna present in the coprolite originates from a bison schoetensacki specimen and strongly suggest that the clade  <dig> itself  corresponds to the bison schoetensacki species.fig.
<dig> phylogenetic position of the bison schoetensacki mitochondrial dna sequence.
maximum likelihood phylogenetic analysis was performed using the concatenated sequence fragments  of the bison schoetensacki siréjol sample and the orthologous sequences of the indicated bison mitochondrial genomes.
the tree with the highest log-likelihood is shown drawn to scale, with branch lengths established from the numbers of substitutions per site.
the percentages of trees in which the associated taxa clustered together are displayed next to the branches .
the tree was rooted using the sequence of bubalus bubalis.
see legend to fig.
<dig> for additional details
analysis of bison schoetensacki nuclear sequences
soubrier et al.
<cit>  analyzed a set of ~  <dig>  genome-wide bovine snps from several specimens of their bison bonasus-related cladex , from two bison bonasus specimens  and from two bison priscus specimens   <cit> .
in order to compare the bison nuclear sequences of the coprolite to these data, the library reads were mapped against the taurine cattle reference genome and the genotypes of the bovine snps were determined.
the genotypes of the snps for the nuclear genome of a bison bison specimen were also extracted from published genomic data .
analysis by multidimensional scaling  shows that the closest specimen to bison shoetensacki is the cladex specimen a <dig>  <cit> .
this result corroborates the conclusions drawn from the mitochondrial sequence data that position the bison sequences of the coprolite in the cladex/bb <dig> sister group to bison bonasus.fig.
<dig> multidimensional scaling analysis of bovine nuclear snps.
we used the snps genotyped herein in the bison schoetensacki sequences of the coprolite, in a modern bison bison specimen and in the reference genome of bos primigenius  as well as the snps genotyped by soubrier et al.
<cit>  in the cladex a <dig> specimen, in the historical bison bonasus a <dig> specimen, in the ancient bison bonasus a <dig> specimen, and in the two bison priscus a <dig> and a <dig> specimens
CONCLUSIONS
the most ancient specimens of bison schoetensacki in the fossil record date back to the beginning of the middle pleistocene , but the continuing presence of this species up to the upper pleistocene has been a matter of controversy.
our studies of a siréjol bone fragment provide strong support to the notion that this species existed for a long period of time, extending to the upper pleistocene.
moreover, comparison of our mitochondrial and nuclear genome data with those of recent studies suggest that  <dig>  to >  <dig> -year-old bison schoetensacki remains are present in a number of eurasian sites.
the dna data show that bison schoetensacki positions in a bison clade referred to as cladex  <cit>  or bb <dig>  <cit>  which should be renamed accordingly.
scenarios of bison evolution considered that bison bonasus evolved either from bison priscus or from bison schoetensacki lineages.
the genetic data clearly support the latter rather than the former hypothesis.
the divergence date between bison bonasus and bison schoetensacki is nevertheless still unclear since it has been estimated around  <dig>  kya and around  <dig>  kya in the two studies with different datasets  <cit> .
opposite conclusions have also been reached regarding the evolution of bison bonasus as due either to incomplete lineage sorting  <cit>  or to hybridization between bison priscus and ancestors of bos primigenius  <cit> .
progress toward a precise assessment of the timing of bovinae evolution and of the respective contribution of incomplete lineage sorting and introgression into the emergence of new species calls for the analysis of specimens much older than those studied up to now.
ancient dna analysis was long considered to be only feasible for upper pleistocene specimens, but studies carried out on middle pleistocene remains retrieved from permafrost  <cit>  or even sites from temperate climate  <cit>  indicate that genome data for more than  <dig> kya-old bovinae specimens are expected in the near future.
additional file
additional file 1: supplementary text.
supplementary methods, supplementary results, and supplementary references.
figure s <dig>  real-time pcr analysis of the grotte-aux-ours coprolite and the siréjol bone fragment for mitochondrial dna.
pcr was carried out using a taqman assay designed using the grotte-aux-ours bovine mitochondrial genome sequence .
figure s <dig>  correlation between the coverage by illumina reads and the gc content of the gaoseq_bovinae sequence.
coverage for each position of the gaoseq_bovinae sequence by the  <dig>  unique illumina reads.
the gc content of gaoseq_bovinae was calculated at each base pair  using a 100 bp sliding window.
the normalized coverage by illumina reads and the normalized gc content are shown on the same graph.
figure s <dig>  characteristics of the illumina reads used to assemble the mitochondrial genome gaoseq_bovinae.
length distribution of the  <dig>  unique illumina reads.
positions of the mismatches between the  <dig>  unique illumina reads and the gaoseq_bovinae sequence.
figure s <dig>  maximum likelihood phylogenetic tree of the mitochondrial genomes of crocuta crocuta.
table s <dig>  pcr primers used to amplify bovine mitochondrial dna fragments from the grotte-aux-ours coprolite.
table s <dig>  pcr primers used to amplify bison schoetensacki mitochondrial dna fragments from the siréjol bone sample.
table s <dig>  the set of mitochondrial genomes against which the illumina reads were aligned to determine the identity of the coprolite producer and of its prey.
table s <dig>  ams radiocarbon and stable isotope data for bison schoetensaki specimen  <dig>  table s <dig>  annotation of the bison schoetensacki mitochondrial genome sequence.
table s <dig>  the set of mitochondrial genomes used for phylogenetic analyses.

 during the middle and upper pleistocene, the large bovidae in europe and in northern asia chiefly included the aurochs, bos primigenius, and two bison species, the steppe bison and the woodland bison. the steppe bison, bison priscus, was very common and exhibited a wide geographic distribution stretching from western europe, through central asia and beringia, and into north america. it was a formidable animal with long horns and robust legs; it stood up to two meters at the withers and reached a total length of almost three meters. it occupied cool, steppe-like grasslands. bison priscus became extinct in europe at the end of the last ice age, about  years ago. the woodland bison, bison schoetensacki appeared in the early middle pleistocene. its size was almost as large as that of bison priscus but its leg bones and metapodials were slenderer than those of bison priscus. moreover the horns of bison schoetensacki were about 30% shorter than those of bison priscus and displayed a slightly different shape: the initial third part is oriented forward, the second third goes upward, and the final part slightly bends backward. bison schoetensacki remains are often associated with forest biotopes, hence its vernacular name. the geographic distribution of bison schoetensacki extended from western europe to the south of siberia, but unlike bison priscus, bison schoetensacki was absent from beringia and north america. the species was first described in germany where it has been recorded in several lower middle pleistocene sites, notably in mosbach, mauer, süßenborn and voigtstedt, but also in roter berg and quedlinburg. it has also been found in england, in moldova, in russia, and in north sea deposits. in france, bison schoetensacki is recorded in several lower middle pleistocene sites and in late middle pleistocene deposits in châtillon-saint-jean. it seems to be still present during the upper pleistocène at chasse-sur-rhône, poleymieux and in the herm and siréjol caves. however the persistence of bison schoetensacki during the upper pleistocene is under debate, and siréjol bone remains have been alternatively ascribed to bison schoetensacki and to a bison priscus novel subspecies. bison schoetensacki remains from siréjol display a cranial frontal no larger than 214 mm, short horns of the shape described above with maximal vertical and anteroposterior diameter at their base of and 65 mm, respectively.fig. cave sites and samples. black dots indicate sites where bison schoetensacki remains have been reported in previous studies. red characters indicate the two cave sites that yielded the samples analyzed in the present study. the coprolite from the grotte-aux-ours cave. the bison schoetensacki cannon bone sample from the siréjol cave two bison species exist today: the american bison and the european bison. the genetic history of the american bison has been extensively investigated through the analysis of mitochondrial dna sequences. the recent reconstruction of the complete mitochondrial genome of a -year-old bison priscus specimen confirmed that the closest extant mitochondrial genome sequences for bison priscus are those of bison bison. these data, in combination with nuclear single nucleotide polymorphism genotyping, have established that bison priscus and bison bison are sister groups. the european bison, bison bonasus, is found in europe and the caucasus, where it has been reintroduced after its extinction as a wild species in the early 20th century. the bison bonasus genomes reflect a complex descent. indeed, the bison bonasus nuclear genome is closely related to that of bison bison, in agreement with morphological evidence and the fact that the two bison species can produce completely fertile hybrid offspring. however, the mitochondrial genomes of bison bonasus specimens are more similar to aurochs and cattle genomes than to bison bison genomes. these observations could be explained either by incomplete lineage sorting of the mitochondrial genome or by a scenario according to which the nuclear dna of an ancient population probably related to the extinct aurochs would have been changed by the systematic introgression of bison bulls. the phylogenetic relationships between bison bonasus, bison priscus and bison schoetensacki are subject to debate. some authors suggest that bison schoetensacki could be the ancestor of bison bonasus, but others consider that bison bonasus is derived from an unknown form of bison priscus. the lack of genomic data for bison schoetensacki has so far prevented any conclusive view on this point. as part of an effort to obtain genetic data on animals consumed by the extinct cave hyena, we initiated the analysis of cave hyena coprolites. in the present study, we found a coprolite containing mitochondrial bovinae dna that was different from bos primigenius, bison priscus and bison bonasus mitochondrial sequences. because the cave site that yielded this coprolite is located close to the siréjol cave, we also analyzed dna extracted from a siréjol bison. schoetensacki bone sample. we report here on the mitochondrial dna sequences obtained from the coprolite and the bone remain, on their phylogenetic positions, and on a set of snps of the corresponding nuclear genome sequence. methods archeological samples the cave hyena coprolite was retrieved from the grotte-aux-ours cave. the cave was discovered in and owes its name to the abundant evidence for long-term occupation by the extinct cave bear, including cave bear bones, hibernation nests, claw marks on the walls, and footprints on the paleosurface. occupation by the extinct cave hyena is documented by the presence of coprolites and a cave hyena skull. the bison schoetensacki bone fragment originated from the siréjol cave. this cave, discovered in is located 15 km away from the grotte-aux-ours cave. scientific excavations carried out from to in the siréjol cave yielded bison schoetensacki bone remains that are currently stored in the centre de conservation et d’Étude des collections du musée des confluences. for the purpose of the present study, we analyzed a bison schoetensacki cannon bone fragment. contamination issues to guarantee that genuine ancient dna was recovered from the samples, dna extraction was carried out from material derived from the inside of the samples, and the coprolite and bone dna extracts were obtained on different days. additional precautions to avoid contamination have been described elsewhere. negative controls included mock extracts and pcr blanks, which always failed to yield any amplification product. dna extraction the instruments used to cut and pulverize the samples were cleaned in an ultrasonic water-tank, incubated 1 h at 40 °c in % sds, rinsed times in ultrapure dnase/rnase free distilled water, and dried using absolute ethanol. the coprolite and the bone sample were cut using a carbon-free circular saw. for dna extraction we used  g of biological material retrieved from the center of the coprolite or the bone medulla using a sterile single-use scalpel. the dissected material was pulverized using a mikro-dismembrator s set to a shaking frequency of  rpm for 30 s, transferred into 6 ml of dna extraction buffer, % sds, 65 mm dtt,  mg/ml proteinase k), and incubated 16 h at 42 °c under constant agitation. after centrifugation, the supernatant was recovered, extracted once with one volume of phenol, once with a phenol-chloroform mixture, and once with chloroform. the aqueous phase was dialyzed and concentrated using a 2 ml, centrifugal ultrafilter with a 30 kda cutoff, the column was washed times with distilled water, and the dna extract was recovered in a volume of to 120 μl. it was purified further using a silica-membrane-based method designed for small dna fragments. the final extract had a volume of 100 μl. pcr and sanger sequencing dna analysis pcr reactions were performed in a 50-μl reaction volume containing –  μl of mock or ancient dna extracts, nm of sense and antisense primers, 200 μm dntp,  mm mgcl 5 μl of geneamp 10x pcr buffer ii, and u of amplitaq gold dna polymerase. the primers used for experiments designed to complete the bovine mitochondrial genome sequence assembled from the grotte-aux-ours coprolite are listed in additional file 1: table s the primers used to analyze mitochondrial genome fragments of bison schoetensacki bone sample from the siréjol cave are listed in additional file 1: table s primer pair from additional file 1: table s was used for initial screening of the samples for bovinae dna content. pcr procedures consisted of an enzyme activation step, followed by a single round of pcr cycles, 20 s; 70 °c, 1 min) performed in a veriti thermal cycler. the full reaction volume was loaded onto an 8% polyacrylamide gel stained with sybr green i. only the coprolite and bone extracts yielded amplification products of the predicted size. pcr amplicons were eluted from the gel and cloned using the topo ta cloning kit. sequencing analysis of cloned dna fragments was performed on abi or 3130xl automatic dna sequencers using m forward primer and bigdye terminator chemistry. all mitochondrial genome fragments analyzed by pcr were characterized through 2– pcr replicates and a consensus sequence for each replicate was derived from several clones. for the bovine mitochondrial genome sequence reconstructed from the grotte-aux-ours sample, each individual consensus was added to the illumina sequence data. for the siréjol bone sample, several individual consensuses of the same fragments were combined to deduce the most relevant base of the positions of the bison schoetensacki mitochondrial genome sequence investigated by pcr. generation and sequencing of a library of dna fragments a library of dna fragments suitable for high-throughput sequencing with the illumina procedure was generated from the coprolite dna extract using the illumina truseq nano dna lt sample kit fc-121- we followed the manufacturer’s recommendations, except for modifications that were introduced for the purpose of analyzing ancient, i.e. highly fragmented dna fragments. thus, the initial step of the library construction process was omitted, and reaction products were purified using qiagen columns or polyacrylamide gel electrophoresis instead of magnetic beads. the process of library construction consisted of steps. first, 6 μl of coprolite dna was 5’end-phosphorylated and blunt-ended using a mixture containing t polynucleotide kinase, t and klenow dna polymerases, and the reaction product was purified on a qiagen column. second, a 3’ adenine residue was added to the blunt-ended dna fragments using klenow 3’ to 5’ exo- polymerase, and the enzyme was heat-inactivated at 70 °c. third, illumina adapters with an overhanging thymine were ligated to the dna fragments; the reaction product was purified using a qiagen column and recovered in a volume of 30 μl referred to as the dna library. fourth, a 5-μl aliquot of the library was pcr-amplified using phusion dna polymerase. fifth, the full pcr reaction volume was loaded on an 8% polyacrylamide gel stained with sybr green i, and the to 220-bp long dna fragments were cut off the gel and purified. dna sequencing was performed at genoscope on the illumina hiseq platform using hiseq v chemistry with a read length set to nucleotides and analysis on the single read mode. illumina sequencing data analysis reads were trimmed for adapter sequences, n’s and low quality stretches on the 3’ end, using a software based on the fastx-toolkit package and designed by genoscope. after this step, sequences shorter than nucleotides were discarded using an in-house python script, yielding a dataset of, dna reads. initial analysis of the illumina library content to confirm the identity of the coprolite producer and gain insight into the animal’s diet, the, illumina reads were aligned simultaneously to a set of mitochondrial genomes, including the reference mitochondrial genomes for the extinct cave hyena and for several potential prey species. the alignment was performed using bwa version.. only the reads matching perfectly, without indel or mismatch, to a unique genome were taken into account. assembling simultaneously the bovine and the cave hyena mitochondrial genomes the sequence assembly was performed in several steps. first the illumina reads, ranging in size from to nucleotides, were aligned simultaneously to the reference mitochondrial genomes of the extant european bison and extinct cave hyena using bwa version. with default options except for -o and -n, which were set to and, respectively, allowing for a maximum edit distance that increases with read length, from one mismatch for 20– nucleotide-long reads to four mismatches for 88– nucleotide-long reads. in this first round, the number of reads mapping to the bison bonasus genome was of the same order of magnitude as the number of reads mapping to the crocuta crocuta genome. of significance, reads mapped to both genomes, which justified our strategy to simultaneously assemble the two mitochondrial genomes so as to eliminate spurious alignments. a set of unique reads mapping exclusively on the bison bonasus genome with a mapping quality higher than was selected. the reads provided a median coverage of of the bison bonasus reference sequence and left positions without coverage. a provisional consensus sequence termed bb1seq was then derived from the unique reads. bb1seq differed from the bison bonasus reference genome at positions. in this study a consensus for a given position was based on at least two concordant unique sequences with a minimal mapping quality of for illumina reads. the most frequent base was taken as the consensus. the few cases where more than 25% of the bases from the aligned reads differed from the consensus were manually examined: the large majority of the alternative bases could be ascribed to damage-induced g to a substitutions at the 3' end of the illumina reads and the other cases were checked by sequence data obtained from pcr studies. the cave hyena mitochondrial genome was assembled in parallel, using the same procedure. a provisional consensus sequence termed cc1seq was derived from unique reads that aligned specifically to the crocuta crocuta reference mitochondrial genome. the reads gave a median coverage of for the crocuta crocuta reference sequence and left positions uncovered. cc1seq and the crocuta crocuta reference genome differed at positions. in a second step, all illumina reads were realigned simultaneously to the provisional consensus sequences bb1seq and cc1seq and new consensus sequences, termed bb2seq for the bovine genome and cc2seq for the cave hyena genome, were derived following the procedure previously described. for the bovine sequence, pcr experiments were then carried out to fill in the gaps and derive a robust sequence at positions where only a single read was available. the pcr primer pairs were designed using the bb2seq sequence. a total of consensus pcr sequences were used to complement bb2seq information and to derive a new consensus sequence, termed bb3seq. all illumina reads were again simultaneously aligned to cc2seq and bb3seq. after a last iteration, we were able to establish the final, complete consensus sequence, termed gaoseq_bovinae, for the bovine mitochondrial genome, based on unique illumina reads, and a partial sequence for the cave hyena mitochondrial genome, called gaoseq_crocuta, based on unique illumina reads. annotation of the gaoseq_bovinae mitochondrial genome the annotation of gaoseq_bovinae was derived from the annotation of the bison bonasus reference mitochondrial genome. briefly, gaoseq_bovinae and the bison bonasus reference mitochondrial genome were aligned using the lagan software with the annotation data of the bison bonasus genome. lagan produced a list of the aligned genomic features with their coordinates in the two genomes, which was manually checked. the occurrence of a start and stop codon at the beginning and the end of each coding sequence, and the absence of internal stop codons were verified using an in-house python script. phylogenetic analyses all phylogenetic analyses were performed using mega. the phylogenetic relationships between gaoseq_bovinae and the complete mitochondrial genomes of specimens related to bison bonasus were inferred using the maximum likelihood, the minimum evolution and the neighbor-joining methods. the ml method was based on the general time reversible model. a discrete gamma distribution was used to model evolutionary rate differences among sites. the analysis involved a total of positions. a similar approach was used to analyze the dna sequences of the bison schoetensacki bone sample. sixteen consensus sequences, each one derived from several clones of several pcr replicates, were aligned to the set of bison sequences previously used for the analysis of the complete mitochondrial genomes. the aligned parts were subsequently concatenated and a phylogenetic tree was constructed using the ml method based on the tamura-nei model and a discrete gamma distribution to model evolutionary rate differences among sites. analysis of the bovine nuclear sequences in order to compare the bison nuclear sequences of the coprolite to the genomic data made recently available for ancient specimens related to bison bonasus, the library reads were mapped against the taurine cattle reference umd using bwa, and the bovine snps were called using samtools and bcftools v. the calls were filtered using the same criteria as in. the final number of snps called for the bovine sequences of the coprolite was, which is comparable to the number of snps called for the a sample, but lower than the number of snps called for the other bison samples. these data were completed by extracting the sequences of snps from the set of unplaced genomic scaffolds corresponding to the nuclear genome of bison bison bison isolate tamuid. multidimensional scaling analysis was performed using the r software v. and in particular the cmdscale() function. ams 14c dating radiocarbon dating was attempted on both the hyena coprolite and the bison schoetensacki bone sample. only the bone sample yielded ams 14c measurements. no bone, hair, or macrofloral remains were present in the coprolite, and what constituted the matrix dissolved during weak-alkali extraction. consequently, no direct 14c measurements were possible on the coprolite. part of the siréjol bone sample was broken into 3–4 mm fragments, which were washed in deionized water with brief sonication to remove adhering sediment. the bone fragments were subsequently decalcified with repeated changes of  n hcl for 2 days at 4 °c and the acid-insoluble collagen washed to neutrality with di water. the collagen was treated with  m koh overnight at 4 °c to remove humic acids. an aliquot of the alkali-extracted fraction was washed to neutrality with di water, acidified with  n hcl, freeze dried and used for 14c dating of the koh-extracted, decalcified collagen fraction. the remainder of the koh-extracted collagen was used for gelatin extraction, which comprised heating the alkali-extracted collagen to 90 °c in  n hcl until dissolution occurred, in 45 min. the gelatin solution was filtered through a  μm durapore filter membrane and freeze dried. approximately 5 mg of collagen and gelatin were weighed into quartz tubes with cuo and ag, the tubes evacuated to <  mtorr under ln pumping, sealed with a h2/o torch and the evacuated tubes combusted at 850 °c for 1 h. the co was isolated cryogenically, converted into graphite using the h2-fe method and used for 14c dating of the gelatin fraction. both fractions were dated at the direct ams accelerator facility, seattle, washington usa. the calibrated age of the sample, expressed in calendar years before present, was calculated using oxcal software and the intcal calibration curve.  analysis of the coprolite collected from the grotte-aux-ours cave figure 1b-c shows the coprolite collected in the grotte-aux-ours cave and the bison schoetensacki bone fragment from the siréjol cave. we obtained for the bone fragment radiocarbon dates of  ±  and  ±  bp using the koh-extracted decalcified collagen and the gelatin fraction of the sample, respectively. we considered the average value of these two measurements as the most reliable determination of the sample age, which corresponds to –  cal yr bp. previous analyses performed using conventional 14c dating of bulk bone material from the same cave sector yielded ages of bp and bp. our ams radiocarbon studies agree with this time frame and greatly increase the accuracy of dating analysis of the siréjol material. dating analysis of the coprolite failed because the sample completely dissolved during the pretreatment for ams measurement.fig. uncalibrated and calibrated age of the siréjol bison schoetensacki bone sample. the calbp data correspond to % confidence interval of the sample age using intcal calibration curve pcr analysis of the coprolite dna for the mitochondrial cytb fragments demonstrated the presence of cave hyena dna and of a bovinae sequence different from all of those recorded in genbank up to march the 82-bp bovinae cytb fragment displayed to mismatches with the bison bonasus, bison priscus, and bos primigenius mitochondrial reference sequences. analysis of the siréjol bone fragment with the same bovinae primers yielded a sequence identical to the one obtained from the coprolite. these pcr studies, and subsequent quantitative pcr analyses also indicated far higher amounts of bovine mitochondrial dna in the coprolite than in the bone sample. we therefore decided to perform high-throughput sequencing of the coprolite dna extract and to undertake a pcr approach for analyzing the bone sample. a total of, illumina single-pass reads, each at least nucleotides in length, were produced from the coprolite dna library. the reads were aligned simultaneously to a set of mitochondrial genomes, including the reference mitochondrial genomes for the extinct cave hyena and several species likely to be part of a carnivorous diet. only the reads matching perfectly, without indel or mismatch, to a unique genome were considered. as shown in additional file 1: table s this analysis confirmed that the extinct cave hyena was indeed the producer of the coprolite, since the number of reads mapping specifically to the cave hyena genome was by far the largest. a large number of reads mapped to the european bison genome. the number of reads mapping specifically to other genomes then plummeted: reads for bos primigenius, reads for bison priscus, reads for the chamois, reads for the musk ox and less than reads for the other genomes tested. of importance, only one read mapped specifically to the reference mitochondrial genome of homo sapiens, attesting to the absence of contamination of the library by human dna. these results thus indicated that the coprolite producer was a cave hyena that had ingested a bovine specimen whose closest known relative is bison bonasus. assembling the bovine mitochondrial genome initial analyses had revealed that the coprolite had been produced by a cave hyena that preyed on a bovine related to bison bonasus, hereafter referred to as "the gao bovinae". we set out to reconstitute the mitochondrial genomes of both specimens by aligning the illumina reads to the reference mitochondrial genomes of bison bonasus and crocuta crocuta. the reconstruction of the two genomes was carried out in parallel so as to eliminate spuriously aligned reads. we assembled a -bp bovine mitochondrial genome, termed gaoseq_bovinae, from unique illumina reads complemented with consensus pcr sequences. we obtained for gaoseq_bovinae a 32-fold median coverage, taking into account both the illumina reads and the consensus pcr sequences. as shown in additional file 1: figure s coverage by illumina reads was quite variable along the genome and strongly correlated with the sequence gc content, as already observed in other cases for illumina sequencing. the length distribution of the unique illumina reads is shown in additional file 1: figure s3a. the median read length was base pairs, as expected for ancient, highly fragmented dna. we analyzed the differences between the illumina reads and the consensus sequence gaoseq_bovinae. only mismatches out of aligned bases were observed between the unique illumina reads and gaoseq_bovinae, which corresponds to % identity. analysis of the positions of these differences with respect to the read ends showed that 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. an increased substitution rate at the read ends is considered as a hallmark for authenticating sequences generated from ancient dna fragments. such fragments exhibit inflated cytosine deamination rates at 5'-overhangs, responsible for g to a transitions at the 3’ end of the opposite, repaired dna strand. gaoseq_bovinae is a complete circular genome that contains the expected number of genes for a mammalian mitochondrial genome, with full-length coding sequences for all protein-coding genes. with respect to the reference mitochondrial genome of bison bonasus, gaoseq_bovinae exhibits differences, including indels and substitutions with a large majority of transitions. in comparison, the three available mitochondrial genomes of bison bonasus display less than differences between each other. while our manuscript was under review, two articles were published that describe several mitochondrial genomes from late pleistocene/holocene bison specimens related to bison bonasus. we therefore performed a phylogenetic analysis including all these ancient mitochondrial genomes along with the mitochondrial genomes of bison bonasus specimens. we also considered bison priscus specimens from these studies and a previous report from our laboratory, as well as the reference mitochondrial genomes of bos primigenius, bison bison, bos grunniens and bubalus bubalis. phylogenetic trees were constructed from this dataset with the maximum likelihood, minimum evolution and neighbor-joining methods that all yielded similar topologies. as shown in fig.  for the ml method, gaoseq_bovinae positions in a well-supported clade that includes all the ancient genomes forming cladex in and clade bb in. no species from the fossil record has been proposed up to now for this clade which delineates a sister species for bison bonasus ancient, historical, and extant specimens from clade. mitochondrial genomes of clades and are phylogenetically closer to extant and ancient bos primigenius than to bison priscus and extant american bison.fig. maximum likelihood phylogenetic tree of complete mitochondrial genomes of bison and bos species. the tree with the highest log-likelihood is shown drawn to scale, with branch lengths established from the numbers of substitutions per site. the percentages of trees in which the associated taxa clustered together are displayed next to the branches. the genbank accession numbers of the sequences are listed in additional file 1: table s specimen name and country of origin are indicated for ancient and historical samples. the tree was rooted using the sequence of bubalus bubalis. refer to text for explanation about clade and clade 2 analysis of mitochondrial dna sequences derived from a bison schoetensacki bone sample we explored the phylogenetic relationship between clade and bison schoetensacki sequences by analyzing mitochondrial dna fragments of the siréjol bone sample. for this we designed a set of primer pairs allowing pcr amplification of short mitochondrial dna fragments that were analyzed through the sequencing of cloned amplicons from several pcr replicates. concatenation of the consensus sequences of the pcr fragments yielded information on 609 bp of the mitochondrial genome. the phylogenetic tree obtained using this concatenated sequence and orthologous fragments of bison mitochondrial genomes strongly supports the position of the siréjol bison schoetensacki specimen within clade and illustrates the similarity between the bison schoetensacki sequence and the sequences of arq and gaoseq_bovinae, from which it differs at and positions, respectively. by contrast, the bison schoetensacki sequence displays to differences with the corresponding fragments of the mitochondrial genomes belonging to clade including the modern bison bonasus sequences, and and differences with the reference mitochondrial genomes of bison priscus and bison bison. these observations support the conclusion that the bovine dna present in the coprolite originates from a bison schoetensacki specimen and strongly suggest that the clade itself corresponds to the bison schoetensacki species.fig. phylogenetic position of the bison schoetensacki mitochondrial dna sequence. maximum likelihood phylogenetic analysis was performed using the concatenated sequence fragments of the bison schoetensacki siréjol sample and the orthologous sequences of the indicated bison mitochondrial genomes. the tree with the highest log-likelihood is shown drawn to scale, with branch lengths established from the numbers of substitutions per site. the percentages of trees in which the associated taxa clustered together are displayed next to the branches. the tree was rooted using the sequence of bubalus bubalis. see legend to fig.  for additional details analysis of bison schoetensacki nuclear sequences soubrier et al. analyzed a set of ~  genome-wide bovine snps from several specimens of their bison bonasus-related cladex, from two bison bonasus specimens and from two bison priscus specimens. in order to compare the bison nuclear sequences of the coprolite to these data, the library reads were mapped against the taurine cattle reference genome and the genotypes of the bovine snps were determined. the genotypes of the snps for the nuclear genome of a bison bison specimen were also extracted from published genomic data. analysis by multidimensional scaling shows that the closest specimen to bison shoetensacki is the cladex specimen a. this result corroborates the conclusions drawn from the mitochondrial sequence data that position the bison sequences of the coprolite in the cladex/bb sister group to bison bonasus.fig. multidimensional scaling analysis of bovine nuclear snps. we used the snps genotyped herein in the bison schoetensacki sequences of the coprolite, in a modern bison bison specimen and in the reference genome of bos primigenius as well as the snps genotyped by soubrier et al. in the cladex a specimen, in the historical bison bonasus a specimen, in the ancient bison bonasus a specimen, and in the two bison priscus a and a specimens  the most ancient specimens of bison schoetensacki in the fossil record date back to the beginning of the middle pleistocene, but the continuing presence of this species up to the upper pleistocene has been a matter of controversy. our studies of a siréjol bone fragment provide strong support to the notion that this species existed for a long period of time, extending to the upper pleistocene. moreover, comparison of our mitochondrial and nuclear genome data with those of recent studies suggest that to >  -year-old bison schoetensacki remains are present in a number of eurasian sites. the dna data show that bison schoetensacki positions in a bison clade referred to as cladex or bb which should be renamed accordingly. scenarios of bison evolution considered that bison bonasus evolved either from bison priscus or from bison schoetensacki lineages. the genetic data clearly support the latter rather than the former hypothesis. the divergence date between bison bonasus and bison schoetensacki is nevertheless still unclear since it has been estimated around kya and around kya in the two studies with different datasets. opposite conclusions have also been reached regarding the evolution of bison bonasus as due either to incomplete lineage sorting or to hybridization between bison priscus and ancestors of bos primigenius. progress toward a precise assessment of the timing of bovinae evolution and of the respective contribution of incomplete lineage sorting and introgression into the emergence of new species calls for the analysis of specimens much older than those studied up to now. ancient dna analysis was long considered to be only feasible for upper pleistocene specimens, but studies carried out on middle pleistocene remains retrieved from permafrost or even sites from temperate climate indicate that genome data for more than kya-old bovinae specimens are expected in the near future. additional file additional file 1: supplementary text. supplementary methods, supplementary results, and supplementary references. figure s real-time pcr analysis of the grotte-aux-ours coprolite and the siréjol bone fragment for mitochondrial dna. pcr was carried out using a taqman assay designed using the grotte-aux-ours bovine mitochondrial genome sequence. figure s correlation between the coverage by illumina reads and the gc content of the gaoseq_bovinae sequence. coverage for each position of the gaoseq_bovinae sequence by the unique illumina reads. the gc content of gaoseq_bovinae was calculated at each base pair using a 100 bp sliding window. the normalized coverage by illumina reads and the normalized gc content are shown on the same graph. figure s characteristics of the illumina reads used to assemble the mitochondrial genome gaoseq_bovinae. length distribution of the unique illumina reads. positions of the mismatches between the unique illumina reads and the gaoseq_bovinae sequence. figure s maximum likelihood phylogenetic tree of the mitochondrial genomes of crocuta crocuta. table s pcr primers used to amplify bovine mitochondrial dna fragments from the grotte-aux-ours coprolite. table s pcr primers used to amplify bison schoetensacki mitochondrial dna fragments from the siréjol bone sample. table s the set of mitochondrial genomes against which the illumina reads were aligned to determine the identity of the coprolite producer and of its prey. table s ams radiocarbon and stable isotope data for bison schoetensaki specimen table s annotation of the bison schoetensacki mitochondrial genome sequence. table s the set of mitochondrial genomes used for phylogenetic analyses. 
in this paper, we analyzed a -year-old bison schoetensaki bone sample from the siréjol cave and a cave hyena coprolite found in a nearby cave and containing large amounts of bovinae dna. despite extensive genetic studies carried out on the steppe and european bison, no remains from the fossil record morphologically identified as bison schoetensacki has been analyzed up to now. this clade only includes ancient bison specimens without taxonomic identification and displays a sister relationship with the extant european bison. bison schoetensacki remains are documented in eurasian middle pleistocene sites, but their presence in upper pleistocene sites has been questioned. this work provides genetic evidence supporting the continuing presence of bison schoetensacki up to the upper pleistocene. bison schoetensacki turns out to be a sister species of bison bonasus, excluding the steppe bison bison priscus as a direct ancestor of the european bison. the european bison, now found in europe and the caucasus, has been proposed to originate either from the extinct steppe/extant american bison lineage or from the extinct bison schoetensacki lineage. the genetic proximity of bison schoetensacki with specimens from this clade is corroborated by the analysis of nuclear dna single nucleotide polymorphisms.
 the european bison, now found in europe and the caucasus, has been proposed to originate either from the extinct steppe/extant american bison lineage or from the extinct bison schoetensacki lineage. bison schoetensacki remains are documented in eurasian middle pleistocene sites, but their presence in upper pleistocene sites has been questioned. despite extensive genetic studies carried out on the steppe and european bison, no remains from the fossil record morphologically identified as bison schoetensacki has been analyzed up to now.  in this paper, we analyzed a -year-old bison schoetensaki bone sample from the siréjol cave and a cave hyena coprolite found in a nearby cave and containing large amounts of bovinae dna. we show that the bovinae mitochondrial dna sequences from both samples, including a complete mitochondrial genome sequence, belong to a clade recently reported in the literature. this clade only includes ancient bison specimens without taxonomic identification and displays a sister relationship with the extant european bison. the genetic proximity of bison schoetensacki with specimens from this clade is corroborated by the analysis of nuclear dna single nucleotide polymorphisms.  this work provides genetic evidence supporting the continuing presence of bison schoetensacki up to the upper pleistocene. bison schoetensacki turns out to be a sister species of bison bonasus, excluding the steppe bison bison priscus as a direct ancestor of the european bison. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords ancient dnabovinaecave hyenacoprolitecrocuta crocutamitochondrial genomepaleogenomicscea saclayloc1043elaloelalouf jean-marc owls roost enterprise lpissue-copyright-statement© the author 2017 