Deciphering the Wisent Demographic and Adaptive Histories from Individual Whole-Genome Sequences

Abstract

As the largest European herbivore, the wisent (Bison bonasus) is emblematic of the continent wildlife but has unclear origins. Here, we infer its demographic and adaptive histories from two individual whole-genome sequences via a detailed comparative analysis with bovine genomes. We estimate that the wisent and bovine species diverged from 1.7 × 10⁶ to 850,000 years before present (YBP) through a speciation process involving an extended period of limited gene flow. Our data further support the occurrence of more recent secondary contacts, posterior to the Bos taurus and Bos indicus divergence (∼150,000 YBP), between the wisent and (European) taurine cattle lineages. Although the wisent and bovine population sizes experienced a similar sharp decline since the Last Glacial Maximum, we find that the wisent demography remained more fluctuating during the Pleistocene. This is in agreement with a scenario in which wisents responded to successive glaciations by habitat fragmentation rather than southward and eastward migration as for the bovine ancestors. We finally detect 423 genes under positive selection between the wisent and bovine lineages, which shed a new light on the genome response to different living conditions (temperature, available food resource, and pathogen exposure) and on the key gene functions altered by the domestication process.

Keywords: European bison; adaptation; cattle; demography; domestication; evolution; genetics; systems biology.

Fig. 1

Population size histories inferred from the two wisent and the four bovine genomes under the PSMC model. Backward in time (in YBP) estimates of the effective population sizes derived from the psmc analyses of the BBO_3569 and BBO_3574 individual wisent genomes and the AAN_0037, HOL_0101, JER_0009, and SIM_0043 bull bovine genomes and assuming ${\mu }_y=2.10\mathrm{ }\times \mathrm{ }{10}^{-9}$ (Liu et al. 2006) and g = 6 (Keightley and Eyre-Walker 2000; Gautier et al. 2007). At the bottom of the figure, the timing of cattle domestication and the zebu (BIN) and taurine (BTA) divergence (Ho et al. 2008) are indicated by a dark red triangle while the timing interval of the bovine and wisent divergence (see the main text) is indicated by dark red arrows. Similarly, European ice ages (Würm, Riss, Mindel, and Günz glacial stages) according to the Penck and Bruckner work (Penck and Brückner 1901–1909) cited in Elias (2013) are indicated by blue intervals at the top of the figure. Finally, the LGM, approximately 20,000 YBP (Clark et al. 2009), is indicated by a vertical blue dotted line.

Fig. 2

Characterization of the divergence between the BBO and BTA lineages under the isolation (I-) and the isolation-with-migration (IM-) models. (A) Estimation of the time ${\tau }_D$ of divergence between the BBO and BTA lineages under the I model (that assumes a clean split). The two violin plots show the distribution of ${\tau }_D$ estimates obtained for each 10 Mb nonoverlapping segments from the alignments of the BBO_3569 (n = 123 segments) and BBO_3574 (n = 177 segments) whole-genome sequence onto the UMD 3.1 bovine genome assembly. Time scale was translated into YBP assuming ${\mu }_y=2.10\mathrm{ }\times \mathrm{ }{10}^{-9}$ (Liu et al. 2006). (B) Estimation of starting (${\tau }_2$) and ending (${\tau }_1$) time estimates of the divergence between the BBO and BTA lineages under the IM model (that assumes the two ancestral lineages exchange migrants between ${\tau }_2$ and ${\tau }_1$). The violin plots show the distribution of ${\tau }_2$ and ${\tau }_1$ time estimates obtained for each 10 Mb nonoverlapping segments from the alignments of the BBO_3569 (n = 123 segments) and BBO_3574 (n = 177 segments) whole-genome sequence onto the UMD 3.1 bovine genome assembly. Time scale was translated into YBP assuming ${\mu }_y=2.10\times {10}^{-9}$ (Liu et al. 2006). (C) Estimation of the average number of migrant per generation ($2N_{\mathrm{e}}m$) during the divergence of the BBO and BTA lineages under the IM model. The two violin plots show the distribution of the migration rate estimates obtained for each 10 Mb nonoverlapping segments from the alignments of the BBO_3569 (n = 123 segments) and BBO_3574 (n = 177 segments) whole-genome sequence onto the UMD 3.1 bovine genome assembly. (D) Model comparisons between the I and the IM models. The two violin plots show the distribution of the difference between the Akaike Information Criterion for the isolation (AIC_I) and for the IM (AIC_IM) obtained for each 10 Mb nonoverlapping segments from the alignments of the BBO_3569 (n = 123 segments) and BBO_3574 (n = 177 segments) whole-genome sequence onto the UMD 3.1 bovine genome assembly. Because the model with the smallest AIC should be preferred, the distributions provide strong support in favor of the IM-model (AIC_IM is always lower than AIC_I).

Fig. 3

Representation of the genes with a K_a/K_s >1 connected to their key functions. Each global function and its links to the corresponding key genes are differently colored, that is, in purple for nervous system, blue for by-products of domestication, green for functions related to embryonic, organ development and hereditary disorders, pink for immune and inflammatory responses, orange for reproduction, gray for hair development and thermogenesis, yellow for lipid metabolism, and white for gustation. Global functions contain genes belonging to several Ingenuity Pathway Analysis functional categories. Gene symbols are colored in red and color intensity is correlated to the K_a/K_s value.