Early cave art and ancient DNA record the origin of European bison

Abstract

The two living species of bison (European and American) are among the few terrestrial megafauna to have survived the late Pleistocene extinctions. Despite the extensive bovid fossil record in Eurasia, the evolutionary history of the European bison (or wisent, Bison bonasus) before the Holocene (<11.7 thousand years ago (kya)) remains a mystery. We use complete ancient mitochondrial genomes and genome-wide nuclear DNA surveys to reveal that the wisent is the product of hybridization between the extinct steppe bison (Bison priscus) and ancestors of modern cattle (aurochs, Bos primigenius) before 120 kya, and contains up to 10% aurochs genomic ancestry. Although undetected within the fossil record, ancestors of the wisent have alternated ecological dominance with steppe bison in association with major environmental shifts since at least 55 kya. Early cave artists recorded distinct morphological forms consistent with these replacement events, around the Last Glacial Maximum (LGM, ∼21-18 kya).

Figure 1. Cave painting example of steppe bison-like and wisent-like morphs.

(a) Reproduction from Lascaux cave (France), from the Solutrean or early Magdalenian period (∼20,000 kya—picture adapted from ref. 53). (b) Reproduction from the Pergouset cave (France), from the Magdalenian period (<17,000 kya—picture adapted from ref. 54).

Figure 2. Identification of CladeX.

(a) Phylogenetic tree inferred from bovine mitochondrial control region sequences, showing the new clade of bison individuals. The positions of the newly sequenced individuals are marked in red for CladeX. (b) Bovine phylogeny estimated from whole-mitochondrial genome sequences, showing strong support for the grouping of wisent and CladeX with cattle (cow) and zebu. For both trees (a,b) numbers above branches represent the posterior probabilities from Bayesian inference, numbers below branches represent approximate likelihood ratio test support values from maximum-likelihood analysis and scale bars represent nucleotide substitutions per site from the Bayesian analysis. (c) Maximum-clade-credibility tree of CladeX and wisent estimated using Bayesian analysis and calibrated with radiocarbon dates associated with the sequenced bones. Dates of samples older than 50 kyr were estimated in the phylogenetic reconstruction. (d) Map showing all sampling locations, using the same colour code (red for CladeX, orange for wisent and blue for steppe bison).

Figure 3. Genome-wide data comparison of bison.

(a) Maximum-likelihood phylogeny of modern and ancient bison from ∼10,000 genome-wide nuclear sites, showing the close relationship between CladeX and steppe bison. However, a bifurcating phylogeny is not capable of displaying the complex relationships between these taxa (see Supplementary Fig. 8). Numbers above branches represent bootstrap values. (b) D-statistics from the same ∼10,000 nuclear sites, using sheep as outgroup. For three bison populations, assuming two bifurcations and no hybridizations, three possible phylogenetic topologies can be evaluated using D-statistics, with the value closest to 0, indicating which topology is the most parsimonious. The topology being tested is shown on the vertical axis. Error bars are three s.e.'s (from block jackknife) either side of the data point. Data points that are significantly different from zero are shown in grey. The data point representing the topology in a, among a set of three possible topologies, is shown with a black outline. (c) Principal Component Analysis of ∼10,000 genome-wide nuclear sites (ancient wisent not included due to the sensitivity of PCA to missing data, see Supplementary Fig. 10). (d) Proportion of steppe bison and aurochs ancestry in both wisent and CladeX lineages, calculated with f₄ ratios.

Figure 4. Temporal and geographical distribution of bison in Europe.

Individual calibrated AMS dates from the present study and published data are plotted on top of the NGRIP δ¹⁸O record. Age ranges for infinite AMS dates are from molecular clock estimates (Fig. 2c). Greenland interstadials (GIs) are numbered in black and marine isotope stages (MIS) in grey.