Incomplete lineage sorting rather than hybridization explains the inconsistent phylogeny of the wisent

Abstract

The wisent or European bison is the largest European herbivore and is completely cross-fertile with its American relative. However, mtDNA genome of wisent is similar to that of cattle, which suggests that wisent emerged as a hybrid of bison and an extinct cattle-like species. Here, we analyzed nuclear whole-genome sequences of the bovine species, and found only a minor and recent gene flow between wisent and cattle. Furthermore, we identified an appreciable heterogeneity of the nuclear gene tree topologies of the bovine species. The relative frequencies of various topologies, including the mtDNA topology, were consistent with frequencies of incomplete lineage sorting (ILS) as estimated by tree coalescence analysis. This indicates that ILS has occurred and may well account for the anomalous wisent mtDNA phylogeny as the outcome of a rare event. We propose that ILS is a possible explanation of phylogenomic anomalies among closely related species.

Fig. 1

Phylogenetic relationship between wisent and related species. a Genetic distances between taurine cattle and wisent, taurine cattle and bison, and wisent and bison, respectively, for a 500 kb sliding window across the genome. The trace for taurine cattle vs wisent lines largely coincides with the trace for taurine cattle vs bison. b Maximum likelihood trees of 4,278,251 nuclear fourfold degenerate sites and whole MT genomes with divergence time estimates. The scale indicate in green gives the estimates based on the human mtDNA mutation rate

Fig. 2

Phylogenetic relationships and inference of gene flow between six species. a Tree of allele-sharing distances between 36 individuals based on SNP calling. b Log total length of identical-by-descent haplotypes shared by individuals, the dark red corresponding to e = 454 Mbp. c ABBA/BABA statistic D with the indicated populations. All D values shown are significant (Z-score > 3, and p value < 0.00135). Positive values indicate that species X is more related to the first species than to the second species (in the first panel taurine cattle and aurochs, respectively). The respective wisent populations differ in their degree of gene flow with cattle lineages. d Pairwise genome-wide D_xy genetic distances between taurine cattle and zebu cattle, yak, American bison, and other wisent populations. X indicates the respective species on the X-axis. e D_xy distances between taurine cattle and zebu, different wisent populations, or American bison along chromosome 1. The Caucasian wisent population has the lowest D_xy in most regions, but not in the regions from 60 to 80 Mb; American bison is remarkably close to taurine cattle, which is not observed in other chromosomes

Fig. 3

Mitochondrial genome replacement model. a Phylogeny of the mtDNA genomes of 12 species and per branch the number of mutations, estimated using RAxML. b Hypothetical crossing of two species indicated by red and blue, respectively, proposed to have led to the replacement of mtDNA (small circles). c The ghost species introgression model, leading to the transfer of mtDNA from the ghost species to wisent (red line). d Histogram of wisent-specific mutations. The red line above the histogram represents the distributions generated by 100,000 binomial sampling based on prior parameters (Methods) and is not expected to fit the observed distribution if the wisent genome contains sequences originating from the ghost species

Fig. 4

Incomplete lineage sorting model. a Gene trees with high bootstrap support. Blue lines indicate the most frequent nuclear phylogeny (Fig. 1b). b The empirical frequencies of all 105 possible gene tree topologies ordered according to (1) ΔRF (>0 for gene trees more similar to the nuclear genome tree; 0 for trees equally similar to the nuclear and mtDNA tree; and <0 for trees more similar to the mtDNA tree) and (2) frequency. The four classes of tree topologies are indicated by different colors

Fig. 5

Validation of coalescent simulations. a Comparison of simulated frequencies based on coalescent probabilities of gene topologies of different classes with the corresponding observed frequencies. The mtDNA topology is indicated. b Comparison of the distributions of Robinson–Foulds distances between the species tree and the observed gene trees (blue numbers and solid blue points) and between the species tree and the gene trees expected on the basis of coalescent simulations (violin plot). The data of violin plot are from 1000 separated simulations (each simulation contains 15,836 gene trees). c ILS model showing fixation of different MT genomes in the branches leading to wisent and bison, respectively. d Simulated gene trees with the same topology as mtDNA. e Observed gene trees with the same topology as the mtDNA tree. f Observed mtDNA tree. By scaling all trees in (d, e, f) have the same total length