Complex Admixture Preceded and Followed the Extinction of Wisent in the Wild Open Access

Introduction

The last known wild wisent, or European bison (Bison bonasus), was shot and killed in 1927, marking the extinction of this species in the wild (Pucek 1991). As a result of an intensive captive breeding program and a series of re-establishments, today the species again occupies part of its former range in Central and Eastern Europe. The total population of free-ranging wisent now stands at over 5,000 individuals (Raczyński 2014), and the International Union for Conservation of Nature has downgraded the conservation status of wisent to threatened (Olech 2008).

In historical times, wisent ranged extensively across semi-open habitats and broadleaved, mixed and coniferous forests in Western Europe, from what is today France in the west to the Volga River and the Caucasus in the east, with the northernmost range limits around 60° north (fig. 1; Kuemmerle et al. 2011; Kerley et al. 2012; Bocherens et al. 2015). However, ongoing habitat fragmentation and overhunting eradicated most populations. By the end of the 19th century, there were only two populations of wisent left in the wild that were assigned to separate subspecies: in Białowieża Forest (Lowland wisent, B. b. bonasus) and in the western Caucasus Mountains (Caucasian wisent, B. b. caucasicus). Finally, even these populations collapsed; the last wild Lowland wisent was shot in Poland in 1919 followed by the last Caucasian animal in 1927 (Pucek 1991).

Fig. 1

Map of Western Europe showing the putative historical range of Lowland wisent (shaded green) and Caucasian wisent (shaded grey) based on bone remains and written records (according to Benecke 2005; Kuemmerle et al. 2011; Tokarska et al. 2011; Bocherens et al. 2015) and sample locations. Black circles indicate contemporary free-ranging modern L line herds and white circles indicate modern LC line herds. Purple and peach circles denote the locations of investigated modern L (MdL1, MdL2) and LC (MdLC) line wisent, respectively, orange squares show the location of the Holocene Lowland wisent (Bb1–3) and blue and yellow triangles indicate historical founding wisent from the Pszczyna population (PLANTA and PLATEN) and the extinct Caucasian wisent (Cc1–3), respectively.

Open in new tabDownload slide

In 1924, the captive population consisted of only 54 individuals (29 males and 25 females). However, the actual founding captive population of wisent was considerably smaller, and is thought to comprise of just 12 individuals (Slatis 1960). All but one of these 12 founders were Lowland wisent, almost half of which came from a population established in 1865 in Pless (now Pszczyna, Poland). The remaining founder was a Caucasian wisent bull named M100 KAUKASUS that represented the last surviving pure Caucasian wisent in captivity. The modern herds that are derived from this founding population are managed as two separate genetic lines. The Lowland line (L) derives from seven Lowland founders (4 males and 3 females), and is thus considered to represent a pure Lowland wisent lineage. The Lowland-Caucasian (LC) line originates from all 12 founders (5 males, 7 females), which included the last remaining Caucasian wisent bull (Slatis 1960). Descendants of the LC line thus represent a mixture of Lowland and Caucasian wisent ancestry.

Although the wisent restitution undoubtedly represents a tremendous conservation success, several factors may limit the long-term viability of the species, many of which are applicable to ex-situ conservation strategies in general. A factor that has received particular attention is that of reduced genetic variability, which may be correlated with a lowered resistance to disease and parasites in wisent (Karbowiak et al. 2014a,b; Majewska et al. 2014; Panasiewicz et al. 2015), and also seriously impacts conservation programs for other threatened species (Altizer et al. 2007). Although genetic variability among living wisent herds has been widely investigated using a variety of genetic markers (Gralak et al. 2004; Luenser et al. 2005; Radwan et al. 2007; Wójcik et al. 2009; Babik et al. 2012; Tokarska, Kawałko, et al. 2009; Tokarska, Marshall, et al. 2009; Tokarska et al. 2015), studies at the level of the complete genome have been far less frequent. A recent study (Gautier et al. 2016) presented moderate-coverage (∼10×) genomic data from two modern L line wisent. They found levels of heterozygosity in these individuals to be similar to that found in modern domestic cow breeds. However, levels of genetic variability occurring prior to extinction of wisent in the wild, and any potential loss of genetic variability following the establishment of the captive breeding population, remains unquantified.

A second threat for wisent is potential hybridization with domestic cattle (Bos taurus). Although F1 hybrid bulls are sterile, hybrid female offspring are not (Basrur 1968), and would therefore have the potential to reintegrate back into the wild population. The bison lineage (including wisent, American bison and the extinct steppe bison) is thought to have diverged from the Bos cattle lineage around 0.85–1.7 YBP, and likely included an extended period of geneflow (Gautier et al. 2016). Analysis of genomic data from modern wisent has suggested, in addition, more recent admixture between wisent and the Bos cattle lineage (Gautier et al. 2016). However, the timing of admixture, in particular whether it preceded or followed the establishment of the captive breeding program, remains uncertain. Such information is critical to assess the magnitude of the threat that domestic cattle admixture represents to the long term viability and integrity of living wisent populations.

Here we present low-coverage whole genome sequencing data from modern and historical wisent, including representatives from both modern genetic lines, from the original founding population, and from the extinct Caucasian wisent subspecies. Using these data, we are able to track changes in genetic variability through the wisent extinction and their subsequent restitution. Furthermore, by inclusion of historical individuals, we not only detect admixture between wisent subspecies and between wisent and the Bos cattle lineage, but also place these events along an absolute timescale of pre- or post-extinction age. Our results demonstrate the huge potential of genomic approaches, in particular applied to historical samples, for studying evolutionary histories and also for the conservation management of endangered species.

Results

Sequencing of Wisent Genomes

We conducted shotgun sequencing of wisent genomes using Illumina technology and mapped the resulting sequence reads to the reference genome assembly of the Asian water buffalo, Bubalus bubalis (GenBank accession no. GCA_000471725.1), which represents an outgroup to the wisent/Bos cattle lineage (supplementary fig. S1, supplementary material online; Hassanin et al. 2012, 2013). Although the reference genome assembly of Bos taurus would provide a more complete assembly and evolutionarily less-diverged reference for mapping, we did not use this reference at it may have led to biased estimates of admixture between wisent and the Bos cattle lineage, in particular at lower levels of genomic coverage. In total, we sequenced seven individuals that we divide into three categories. Detailed sample information, including provenance, is provided in table 1 and sample localities are shown in fig. 1.

Modern wisent—three modern individuals representing both genetic lines. These comprise two individuals from the L line (MdL1, mean read depth 1.59×; and MdL2, mean read depth 1.63×; from the Polish and Belarusian parts of the Białowieża Forest, respectively) and one individual from the LC line (MdLC, from Dydiowa in the Bieszczady Mountains, mean read depth 1.49×).

Founding wisent—two individuals assignable to the Lowland wisent subspecies B. b. bonasus from the initial breeding population originating from Pszczyna, both of which contributed to the establishment of both the L and the LC genetic lines: foundress F42 PLANTA (1904–1931, mean read depth 1.82×) and her male offspring, M158 PLATEN (1926–1933, mean read depth 1.36×), who was fathered by another founder M45 PLEBEJER (1917–1937).

Caucasian wisent—two individuals from the early 1900's representing the now extinct Caucasian wisent subspecies B. b. caucasicus (Cc1, from Kuban Oblast, mean read depth 1.17×; and Cc2, from North Ossetia-Alania, mean read depth 0.92×).

For each individual, we collapsed mapped reads into a single pseudo-haploid genome sequence by randomly selecting a single high quality nucleotide from the read stack at each position of the reference genome, following the procedure described by Cahill et al. (2013, 2015). This procedure disregards heterozygous positions, where only one allele will be sampled, but should not introduce any biases in allele sampling. Ancient DNA fragments frequently contain miscoding lesions resulting from postmortem DNA degradation, the most common of which involves the deamination of cytosine to uracil, which causes C to T substitutions in the resulting data (Dabney, Meyer, et al. 2013). This pattern is present at varying levels in sequence data from our historical samples (supplementary fig. S2–S5, Supplementary Material online), and so we restricted all subsequent analyses to transversion sites only to avoid any confounding effects of DNA damage.

Genomic Divergence

We investigated patterns of nuclear genomic divergence among wisent by conducting pairwise comparisons of the number of transversion differences occurring along a sliding window of 1 Mb, producing a distribution of genomic divergence for each wisent pair. The resulting probability densities showed that nuclear genomic divergence is broadly similar among all modern and founding wisent (fig. 2A). The two founding individuals, PLANTA and PLATEN, are somewhat less diverged from one another than either is from all modern wisent (fig. 2A), reflecting their mother-son relationship. Slightly increased divergence is observed between the modern LC line individual (MdLC) and all modern and founding wisent (fig. 2A), which may reflect the increased component of Caucasian wisent ancestry in this individual resulting from the captive breeding program.

Fig. 2

Pairwise genomic divergence among wisent. A, B, C, show probability densities for pairwise transversion divergence (x axes) along a sliding window of 1 Mb. Individual plots show all pairwise comparisons among modern and founding individuals (A); comparisons of modern and founding individuals and Caucasian wisent Cc1 (B); and comparisons of Caucasian wisent Cc2 and all other individuals (C). Specific comparisons discussed in the text are identified by colors, according to the key at the top right of each plot. Schematic neighbor-joining phylogeny based on whole genome distances (D).

Open in new tabDownload slide

Genomic divergence between Caucasian and both modern and founding wisent greatly exceeds that occurring between the latter two groups (fig. 2A–C). Substantial divergence is also found between the two Caucasian wisent individuals. One of these Caucasian wisent (Cc1) was found to be less diverged from modern and founding wisent than other Caucasian wisent individual (Cc2), suggesting the presence of not only substantial genetic diversity but also substantial population structure in the extinct Caucasian wisent subspecies.

We also investigated mitochondrial genome variability among all individuals subjected to nuclear genome sequencing, in addition to seven other modern wisent (supplementary table S1, Supplementary Material online). Sequence analysis revealed that all investigated modern wisent, both founding wisent, and a single historical Caucasian wisent (Cc1), share a single haplotype. The haplotype occurring in the second historical Caucasian wisent (Cc2) differed from this widely shared haplotype by a single transition site. These results hint at a major loss of mitochondrial haplotype diversity prior to the extinction of wisent in the wild. This inference is supported by additional haplotypes that we recovered from three ancient late Holocene Lowland wisent from Austria (the range of calibrated age is from ca. 1.3 kyr to 1.9 kyr) and one ancient middle Holocene Caucasian wisent from Armenia (ca. 5.7 kyr). All these ancient haplotypes are substantially divergent from those found in modern and historical wisent, suggesting a substantial loss of haplotype diversity, potentially within the last ∼1,300 years.

Neighbor-joining phylogenetic analysis of total nuclear genomic divergence supports paraphyly of Caucasian wisent (fig. 2D), as does analysis of mitochondrial haplotypes (supplementary fig. S1, Supplementary Material online). We further investigated the population history of wisent by dividing aligned nuclear genome sequences into non-overlapping 1 Mb blocks and subjecting each block to maximum-likelihood phylogenetic analysis. For this analysis, we included both Caucasian wisent and both modern L line wisent as representatives of the Lowland wisent subspecies, with water buffalo as outgroup. Founding wisent and the modern LC line wisent were not included to avoid any confounding effects of direct ancestor-descendent relationships and documented Caucasian wisent introgression (Slatis 1960) on phylogenetic interpretation, respectively. We found that 57% of the investigated genomic blocks support reciprocal monophyly of Caucasian and Lowland wisent (fig. 3). We therefore conclude that this most likely represents the true population history. All alternative topologies occur, individually, at a much lower frequency. Nevertheless, almost half of the genome sequence alignment of these individuals contradicts the true population history.

Fig. 3

Population history of Lowland and Caucasian wisent, estimated using two representatives of each. The pie chart shows the percentage of 1 Mb genomic blocks supporting each alternatively rooted tree topology. This indicates fraction of genome blocks returning both wisent subspecies as reciprocally monophyletic. Dark and light blue colors show the next most frequently encountered topologies in which Caucasian wisent are paraphyletic (dark blue: Cc2 most divergent, and light blue: Cc1 most divergent, the first of which is compatible with estimates pairwise genomic divergence; see fig 2D).

Open in new tabDownload slide

In order to interpret wisent genomic divergence in the context of total species genetic diversity, we obtained data from the NCBI Short Read Archive for seven domestic cattle and seven yak (Bos grunniens; supplementary table S4, Supplementary Material online) and subjected them to the same analysis pipeline. Genomic divergence among modern wisent was found to be similar to that found among domestic cattle, and exceeded that found among yak (fig. 4). We conducted equivalent comparisons for modern wisent and these other bovid species but with the inclusion of transition as well as transversion sites. Interestingly, the distribution of genomic divergence for pairs of wisent was bimodal in all cases (fig. 4), but the relative levels of genomic divergence between species were similar to that measured using only transversion sites.

Fig. 4

Comparison of pairwise genomic divergence within three bovid species: wisent (Bison bonasus), domestic cattle (Bos taurus), and yak (Bos grunniens). Probability densities were calculated along a sliding window of 1 Mb from transversions only (A), and from transitions and transversions (B). For wisent, only modern samples included.

Open in new tabDownload slide

Wisent Geneflow and Admixture

We investigated patterns of admixture among wisent using the D statistic test for admixture (Green et al. 2010). This test identifies any imbalance in the number of derived alleles that either of two closely related individuals share with a candidate introgressor. A significant excess of derived alleles shared between one individual and the introgressor provides evidence of admixture between them (Durand et al. 2011). For all D statistic tests, we used water buffalo (Bubalus bubalis) as outgroup for allele polarization.

We first investigated patterns of derived allele sharing among modern wisent, and found no statistically significant signal of admixture between the modern LC line individual and either modern L line individual (supplementary table S5, Supplementary Material online). Between modern and founding wisent, we found that modern L line wisent share a significantly greater number of derived alleles with founding wisent than the modern LC line individual does (fig. 5A), indicating a greater contribution of the two founding wisent investigated here to the L line, relative to the LC line, which is consistent with pedigrees (Slatis 1960).

Fig. 5

Results of D statistical analysis. Red and grey points show significant and nonsignificant D values (x axis), respectively, and show: the genetic contribution of the founders (F) to the modern individuals (A); Caucasian wisent (Cc) admixture with modern L and LC herds (MdL and MdLC) and one founder, PLANTA, relative to the least Caucasian admixed wisent, founder PLATEN (B); Caucasian wisent admixture among modern wisent (C); apparent cattle (DC) admixture with all investigated wisent (W) relative to aurochs (Aur) (D); variance in cattle/aurochs admixture among wisent (W) compared to PLATEN (E). Detailed D statistic results are provided in supplementary table S5, Supplementary Material online.

Open in new tabDownload slide

We then investigated admixture involving Caucasian wisent. We found a significant excess of derived allele sharing between one founding wisent, PLANTA, and one Caucasian wisent, Cc2, relative to the other founding wisent, her son PLATEN (fig. 5B). This indicates that a proportion of the genome of PLANTA can be attributed to admixture with Caucasian wisent. Furthermore, we can deduce that PLEBEJER, the father of PLATEN, must have possessed a lower level of Caucasian wisent admixture than PLANTA, and that PLATEN himself was likely admixed to some degree through inheritance from PLANTA. The detection of admixture involving one Caucasian wisent (Cc2) but not the other (Cc1) further supports the existence of genetic structure in Caucasian wisent inferred from estimates of genomic divergence (fig. 2B and C).

Next, we investigated evidence of Caucasian wisent admixture among modern wisent. Consistent with expectations, we found that the modern LC line individual (MdLC) shares an excess of derived alleles with one of the Caucasian wisent (Cc1) relative to modern L line individuals (fig. 5C). We did not, however, detect such an excess between the modern LC line individual and the second Caucasian wisent (Cc2), relative to the modern L line individuals. We can therefore infer that the last surviving Caucasian wisent, KAUKASUS, whose living descendants comprise the modern LC line, was more closely related to Caucasian wisent individual Cc1 that to individual Cc2.

We further investigated Caucasian ancestry in the genome of the modern LC line individual by performing phylogenetic analysis of non-overlapping 1 Mb genomic blocks. The fragmented water buffalo reference genome assembly precludes examination of the variance in ancestry along entire chromosomes. To achieve this, we instead mapped reads to the reference genome of the domesticated zebu cattle, Bos indicus, which was itself generated by mapping short read data to the chromosome-level assembly of B. taurus (Canavez et al. 2012). Any potential biases introduced by using a cattle reference are mitigated by selecting this cattle breed originating from a domestication center on the Indian subcontinent (Loftus et al. 1994), which is geographically distant from the historical distribution of wisent. Analysis of aligned 1 Mb genomic blocks involved the modern LC line wisent (MdLC), the founding wisent (PLATEN) that was found to be least admixed with Caucasian wisent, Caucasian wisent (Cc1) and domestic cattle, with water buffalo as outgroup. Of the investigated genomic blocks, we find that 22% return the modern LC line and Caucasian wisent as monophyletic (fig. 6), and may therefore represent introgressed segments of Caucasian wisent ancestry in this modern LC line individual. Around 8% of these blocks are likely to result from incomplete lineage sorting, based on the frequency of occurrence of the opposing topology (fig. 6), producing an overall estimate of 14% of the genome of the modern LC individual that results from Caucasian wisent admixture, most likely inherited from the bull KAUKASUS. In addition to providing an estimate of admixture proportions, our method is also able to accurately map admixed segments of the genome (fig. 6A). Many of these segments span multiple megabase blocks. For example, a contiguous 22 Mb admixed block is observed on chromosome 4, which may span as much as 33 Mb under the assumption that intervening blocks with missing data are linked to adjacent ones. Relative admixture proportions also vary among chromosomes in this individual. For example, chromosome 27 almost entirely lacks Caucasian wisent ancestry whereas around 50% of chromosomes 4 and 15 are likely derived from admixture.

Fig. 6

Genomic admixture map (A) of Caucasian wisent ancestry in the modern LC line individual (MdLC). Colored blocks indicate 1 Mb genomic blocks returning alternative tree topologies, blue blocks are compatible with the species tree; yellow blocks return the monophly of the modern LC line and Caucasian wisent, and likely result from admixture and to a lesser extent incomplete lineage sorting; red blocks return the monophyly of PLATEN and Caucasian wisent and likely result from incomplete lineage sorting. “X” shows blocks with missing data. The pie chart (B) shows the percentage of 1 Mb genomic supporting each tree topology identified by colors, according to the key presented above.

Open in new tabDownload slide

Finally, we investigated evidence of Caucasian wisent ancestry in the modern L line. We found that both modern L line wisent share a significant excess of derived alleles with Caucasian wisent (fig. 5B), relative to founding wisent. Thus, modern L line individuals appear more admixed with Caucasian wisent than the two founding wisent investigated here. This admixture signal could result from either variable admixture proportions among founding individuals, or recent geneflow between the L and LC lines, although D statistic comparison of modern individuals failed to detect the latter (see above). We further investigated these alternative hypotheses by comparing the sizes of putatively admixed genomic blocks. Recent geneflow results in large contiguous genomic blocks derived from admixture in the genomes of the recipient population, which are broken up over time as a result of recombination. Pedigree information provides an approximate date for geneflow from Caucasian wisent into the modern LC line around 90 years ago (15–22 generations), the result of which are many intact multi-megabase genomic blocks derived from Caucasian wisent in the modern LC line individual (fig. 6). We compared the sizes of these blocks with those of putative Caucasian wisent ancestry in a modern L line individual, and found the abundance of large blocks to be considerably lower in the latter (fig. 7). This rejects recent admixture and instead supports variable admixture proportions among the founding herd in explaining the observed signal of Caucasian wisent admixture in this modern L line individual.

Fig. 7

Variation in the sizes of genomic blocks in modern L (blue) and LC (red) likely resulting from Caucasian wisent admixture. Plots show cumulative probability densities calculated at a scale of 1 Mb. Genomic blocks in the LC line wisent (red) result from admixture occurring around 90 years ago; the lower abundance of larger admixed blocks in the modern L line wisent support that this admixture event preceded the former. The plots have been truncated to aid visualization, and single blocks of 18 and 22 Mb in the LC line individual are not shown. The largest block size detected in the modern L line individual was 8 Mb.

Open in new tabDownload slide

Admixture with the Bos Cattle Lineage

We investigated potential admixture between wisent and the Bos cattle lineage using pseudo-haploid sequences generated from short read data of two domestic Holstein cows and an ancient aurochs (Bos primigenius; Park et al. 2015), the extinct species from which cattle were domesticated and that lived sympatrically with wisent up until its extinction around 400 years ago (van Vuure 2005). First, we looked for significant differences in derived allele sharing between cattle and wisent, relative to the aurochs. We found that all investigated wisent share a significant excess of derived alleles with cattle relative to aurochs (fig. 5D). This suggests either admixture between wisent and domestic cattle, or alternatively, admixture with aurochs, if aurochs populations were highly structured and the admixing individuals were from a population more closely related to domestic cattle than the British aurochs used in this analysis.

We also compared derived allele sharing with cattle among the individual wisent. We found that all modern wisent investigated here share a significant excess of derived alleles with cattle relative to any founding or Caucasian wisent (fig. 5E). Variable admixture was also observed among founding wisent. Specifically, PLANTA shares more derived alleles with domestic cattle than either Caucasian wisent or PLATEN do. Since D statistic is a relative test, and PLATEN is the offspring of admixed founder PLANTA, it is reasonable to infer that PLATEN is also admixed to some extent, however.

Finally, using the $f^$ statistic (Durand et al. 2011), we estimated a fraction of 2.4–3.2% of the genomes of the modern wisent that could be attributed to admixture with the Bos cattle lineage, above that occurring in the founding wisent (supplementary table S6, Supplementary Material online).

This increased admixture signal observed in all modern wisent relative to founding wisent could result from either variable admixture proportions among the founding herd from which all modern wisent are derived, or alternatively, from recent admixture with modern cattle. However, the fact that we do not find evidence of any complete genomic 1 Mb blocks resulting from cattle admixture in the modern LC line individual (fig. 6A) argues strongly against recent cattle admixture, and instead supports variable admixture among the founding herd in explaining the excess of derived alleles shared among domestic cattle and modern wisent.

Discussion

Retracing complex population histories can be challenging. In particular, admixture involving populations or species that are now extinct may be impossible based solely on data from living individuals (Hofreiter et al. 2015). Through the use of low-coverage genomic data from modern and historical wisent, including from the now extinct Caucasian wisent subspecies, we have revealed the complexity of wisent evolution. This complex history involved not only admixture resulting directly from the captive breeding program, but also older processes occurring prior to their extinction in the wild, which included admixture with the Bos cattle lineage (fig. 8).

Fig. 8

Schematic diagram showing inferred admixture among wisent, and among wisent and the cattle/aurochs lineage. Arrows indicate the direction of the geneflow. Black lines indicate admixture between wisent and cattle/aurochs, yellow lines/arrow—between Caucasian and founding or modern wisent respectively, and the blue arrow—from founders to modern wisent.

Open in new tabDownload slide

Wisent Evolution and Admixture

The accepted view of wisent evolution is of two distinct subspecies, Lowland wisent and Caucasian wisent, that both underwent dramatic population declines, with the last few surviving individuals serving as founders of the modern L (Lowland only) and LC (Lowland and Caucasian) lines (Pucek et al. 2004). Our results show that this model is an oversimplification. We find evidence of at least two highly differentiated populations within Caucasian wisent, with one of these showing greater pairwise similarities with Lowland wisent than with the second Caucasian population, at the level of the complete genome (fig. 2). However, analysis of aligned nuclear genomic blocks from four individuals returns Caucasian and Lowland wisent as reciprocally monophyletic across slightly more than half of the genomic alignment (fig. 3), providing support that this topology reflects the true history of population divergence. Thus, among any two Caucasian wisent and any two modern (or founding Lowland) wisent, we may expect that any single locus has only around 50% probability of reflecting the true history of population divergence. Moreover, increased sampling of individuals is likely to further reduce this proportion, potentially to such an extent that, at a given sampling level, the true evolutionary history cannot be untangled from the effects of random drift and more recent admixture. This result reinforces the notion that phylogeny-based interpretation may be inappropriate at the level of the complete genome, and that alternative, more flexible models will be required to keep pace with our ability to generate such data (Hofreiter et al. 2015).

A further implied assumption of the traditional view of wisent evolution is that, with the exception of the Caucasian bull KAUKASUS, all founding wisent represented “pure” Lowland wisent (Pucek et al. 2004; Tokarska et al. 2015). On this basis, the modern L line that is derived only from the latter can also be considered as pure Lowland wisent, referable to the subspecies B. b. bonasus. Our results also demonstrate that the notion of wisent subspecies purity is flawed in the sense that founding Lowland individuals were in fact admixed with Caucasian wisent to varying degrees. We demonstrate this, both directly for the founder PLANTA in comparison to another founder, her son PLATEN, and also indirectly, by the elevated signal of Caucasian wisent admixture in modern L line wisent, relative to these founders, most likely a result of inheritance from other founding individuals not included in this analysis (fig. 7). The notion of subspecies purity therefore disregards the fact that admixture between Caucasian and Lowland wisent almost certainly occurred prior to the extinction of wisent in the wild, and such admixture could therefore be regarded as part of the normal population processes and dynamics of this species.

The notion of subspecies purity has driven efforts to ensure that free-living L and LC herds do not come into contact (Pucek et al. 2004), and also motivated genetic investigations of living populations that may have been recipients of geneflow from the opposing genetic line (Tokarska et al. 2015). The latter study investigated the modern L line population that currently inhabits the eastern, Belarusian part of the Białowieża Forest. Some individuals were found to possess a microsatellite allele that was common among Caucasian wisent but absent in all studied Lowland wisent, which these authors interpreted as evidence of recent admixture with the modern LC line. Although the individual from this population (modern L line, MdL2) that we sequenced does not possess this putative Caucasian wisent allele, we nevertheless detected evidence of Caucasian wisent admixture above that occurring in founding wisent for this individual (fig. 5). However, the small size of admixed blocks, in addition to non-significant D comparisons of modern lines, supports variable Caucasian wisent admixture among founding wisent in explaining this result, which may also account for the occurrence of putative Caucasian wisent alleles in other individuals from this population. Future studies of such individuals using the methodology applied here would provide a robust test of these alternative hypotheses.

Wisent Conservation and De-Extinction

The issue of low genetic variability among living wisent is considered as cause for concern, and has been inferred by several population-level studies using various molecular and biochemical markers, such as blood-group systems and blood serum proteins (Sipko et al. 1995; Gębczyński and Tomaszewska-Guszkiewicz 1987; Hartl and Pucek 1994; Sipko et al. 1997), mtDNA (this study; Tiedemann et al. 1998; Burzyńska et al. 1999; Wójcik et al. 2009; Hassanin et al. 2012), nuclear gene sequences (Burzyńska and Topczewski 1995; Radwan et al. 2007; Babik et al. 2012; Hassanin et al. 2013), microsatellites (Gralak et al. 2004; Luenser et al. 2005; Tokarska, Kawałko, et al. 2009; Tokarska, Marshall, et al. 2009) and SNPs (Tokarska, Marshall, et al. 2009, Tokarska et al. 2015). In apparent contrast to these results, whole-genome heterozygosity in modern wisent has been shown to be similar to that found within domestic cow breeds (Gautier et al. 2016). Our results further support this finding using measures of pairwise genomic divergences among modern and founding wisent, which are approximately equal to or greater than that found between pairs of cattle or yak, respectively. However, the bimodal distributions observed for pairwise comparisons based on transitions and transversions for modern wisent suggests that, among any pair of modern wisent, large chromosomal blocks will show high genetic similarity, which is indicative of inbreeding, and has been shown previously using SNP array data (Pertoldi et al. 2009).

Overall, the wisent captive breeding program, based on a founding herd of just 12 individuals, appears to have succeeded in retaining reasonable levels of genetic variability in modern populations. Nevertheless, the extinction of wisent in the wild was clearly accompanied by a major loss of total genetic diversity. Divergent mitochondrial haplotypes detected in ancient wisent appear to have been entirely lost from modern populations, which our results suggest may possess only a single haplotype. Moreover, substantial pairwise nuclear genomic divergences detected between modern and historical Caucasian wisent indicate a huge loss of diversity following the extinction of the latter population.

A fraction of the genepool of extinct Caucasian wisent survives in the genomes of modern individuals. Our results provide not only a direct measure of this admixture in a modern LC line individual, but also allow us to map with relative accuracy chromosomal segments that are inherited from Caucasian wisent. Although sometimes controversial, the concept of de-extinction has generated considerable interest (Sherkow and Greely 2013), and attempts are currently underway to generate animals that, at least superficially, resemble the quagga (Equus quagga; Heywood 2013), an extinct subspecies of plains zebra, and also the aurochs (van Vuure 2005), by careful selective breeding of their living relatives. In both of these cases, selective breeding and the ultimate success of the project are based solely on morphological criteria. Our study demonstrates that, at least in principle, by generating chromosomal admixture maps for multiple living representatives of the LC line, it would be possible to selectively breed an animal that is, at the genomic level, highly similar to a Caucasian wisent.

Admixture with the Bos Cattle Lineage

Hybridization of wild species with their domesticated close relatives is a subject of considerable discussion and concern for conservation management (Ellstrand et al. 2010). Previous studies have found evidence of admixture with the Bos cattle lineage in modern L line wisent (Gautier et al. 2016). We extend this result by demonstrating that such admixture is widespread across wisent, including both modern genetic lines, representatives of the original founding herd, and also extinct Caucasian wisent. This admixture post-dates the common ancestor of English aurochs and taurine cattle, and involved a representative of the latter lineage. However, the precise identity of the introgressor—aurochs or domestic cattle—is less certain, given the lack of knowledge of population structure in aurochs. Testing these two alternatives would require data from additional aurochs populations from within the core distribution of wisent, and would be a valuable direction for future research.

The timing of admixture also has implications for conservation management. Specifically, the removal of individuals resulting from very recent hybridization may be deemed appropriate (Halbert and Derr 2007). The small size of cattle admixed blocks in modern wisent (at least undetectable at a 1 Mb scale) clearly rejects very recent cattle admixture for the individuals investigated here. Instead, admixture must have occurred prior to the establishment of the captive breeding program, and the admixture signal detected in modern wisent results from inheritance from the founders that were admixed with cattle to varying degrees. Thus, based on the current evidence, cattle introgression appears of low concern for wisent conservation for the following reasons: (1) admixture does not appear to have occurred since the establishment of the captive breeding program, although screening of additional individuals may be desirable to further support this generalization; (2) the number of intervening generations separating living wisent from the F1 hybrids is likely sufficient that all living wisent are admixed to some extent (Chang 1999); and (3) our results may in fact reflect admixture with aurochs, rather than domestic cattle, although this hypothesis requires further investigation.

An Exemplar for the Study of Admixture

The ability to detect admixture is of key importance for both evolutionary and applied conservation studies. However, interpretation of a significant signal of admixture, in terms of both evolutionary inference and the formulation of management strategies, may require information its timing. Using seven low-coverage wisent genomes from both modern and historical wisent we have revealed multiple instances of admixture, but moreover, because the approximate age of introgression of Caucasian wisent into the modern LC line is known, through comparisons of the sizes of likely admixed genomic segments we have inferred the relative ages of other admixture events. This unique historical information, coupled with the ability to recover genomic data from historical samples, establish wisent as an exemplary taxon for the study of admixture in wild populations. As new analytical methods for studying admixture are developed, wisent can serve as a valuable empirical test of both their performance and utility.

Materials and Methods

Complete details of all samples and specimens used in this study are shown in table 1.

Laboratory Methods, Modern Samples

DNA was extracted from tissue samples of three modern wisent using either a DNeasy Blood and Tissue Kit (Qiagen) according to the manufacturer’s instructions (sample MdL1) or by phenol/chloroform extraction (Sambrook and Russell 2001). We mechanically sheared the DNA of the modern samples using a Covaris S220 sonicator to an average fragment length of 500 bp and prepared indexed Illumina libraries from 500 ng of each modern DNA extract using a published double-stranded protocol (Meyer and Kircher 2010) with modifications (Fortes and Paijmans 2015). Library molecules from 450 bp to 1000 bp were then selected using a Pippin Prep Instrument (Sage Science).

Laboratory Methods, Historical Samples

DNA extraction from four museum specimens as well as sequencing library preparation steps preceding amplification were performed in a dedicated ancient DNA laboratory (Evolutionary Adaptive Genomics Group, Potsdam University, Germany). DNA extracts were prepared from horn and bone powder obtained by grinding in a mixer mill (MM 400, RETSCH). DNA extraction followed the protocol of Dabney, Knapp, et al. (2013), except for horn samples where we used a different digestion buffer containing 10mM Tris buffer (pH 8.0), 10 mM NaCl, 5 mM CaCl₂, 2.5 mM EDTA (pH 8.0), 2% SDS (Shapiro and Hofreiter 2012). The museum samples were already fragmented due to degradation, so were not sonicated. We used 25 µl of each DNA extract to construct single-stranded indexed Illumina libraries according to the protocol of Gansauge and Meyer (2013).

Sequencing

Final library concentrations and the distribution of insert sizes were determined using a 2200 TapeStation (Agilent Technologies) and Qubit HS-assay (Thermo Fisher Scientific), respectively. Each library was then sequenced using an Illumina NextSeq 500 instrument. For modern libraries we used a High Output Kit (75 bp paired-end sequencing), for libraries obtained from historical horn samples we used High Output Kits (75 bp single-end and 150 bp paired-end) and each library built from historical bone samples was sequenced separately with High Output Kits (75 bp single-end and paired-end). Full details of sequencing results are provided in supplementary tables S2 and S3, Supplementary Material online. Whole genome shotgun sequencing data produced for this study are available in the NCBI Short Read Archive as SAMN05950802–SAMN05950808.

Data Processing, Mapping and Pseudo-Haploidization

For paired-end data, we trimmed adapter sequences and merged overlapping read pairs using SeqPrep (https://github.com/jstjohn/SeqPrep), requiring a minimum read length of 30 bp (-L 30), minimum overlap of 15 bp (-o 15), and a minimum merge quality of 13 (-q 13). Adapters occurring at the 3′ ends of single-end reads were trimmed using cutadapt (Martin 2011), also requiring a minimum length of reads of 30 (-m 30). We then mapped the resulting data to the zebu (Bos indicus; GenBank accession no. GCA_000247795.2) and water buffalo (Bubalus bubalis; GenBank accession no. GCA_000471725.1) nuclear genomes and wisent mitochondrial genome (GenBank accession no. KY055664) using BWA aln version 0.7.8 (Li and Durbin 2009) with default 0.04 mismatch value. We removed duplicate reads likely resulting from PCR amplification using samtools rmdup (Li et al. 2009). Detailed descriptions of the mapping results are provided in supplementary tables S2 and S3, Supplementary Material online. We then generated pseudo-haploid sequences as described by Cahill et al. (2015) and used these for further analysis.

Pairwise Genomic Divergence

Pairwise genomic divergence was calculated by dividing genomic alignments into non-overlapping 1 Mb blocks and calculating the proportion of transversions, or transitions plus transversions (comparisons of modern individuals only), for each pair of individuals, accounting for the presence of missing data. Blocks with > 75% missing data were disregarded. Probability densities were generated by kernel density estimation in R (R Core Team 2014) using default parameters. Full details of comparative data generated for domestic cattle and yak (data from the NCBI Short Read Archive) are provided in supplementary table S4, Supplementary Material online.

Mitochondrial Genome Analysis

In addition to the mitochondrial genomes generated from the three modern and four historical specimens, we obtained mitogenomes from seven other modern wisent with Sanger technology and from four ancient Holocene wisent individuals using hybridisation capture (supplementary table S1, Supplementary Material online). These ancient samples were radiocarbon dated at either the Oxford University Radiocarbon Accelerator Unit (Oxford, UK) using ultrafiltered collagen and accelerator mass spectrometry (Ramsey, Gigham, et al. 2004; Ramsey, Higham, et al. 2004) or the VERA-Laboratorium Institut für Isotopenforschung und Kernphysik (Vienna, Austria). We calibrated radiocarbon dates using the IntCal13 calibration curve (Reimer et al. 2013) in OxCal v4.2 (https://c14.arch.ox.ac.uk/oxcal/OxCal.html).

DNA was extracted from the Holocene B. b. bonasus samples (Bb1-Bb3) at the Henry Wellcome Ancient Biomolecules Centre (Oxford University, UK), following Shapiro et al. (2004). We extracted the Holocene B. b. caucasicus sample (Cc3) in the specialist Paleogenomics facility at UC Santa Cruz, following Rohland et al. (2010). DNA library construction, mitochondrial target enrichment, sequencing, and sequence data processing protocols for the four Holocene samples followed approach four in Heintzman et al. (2016), except that the whole mitochondrial genome consensus sequence was retained. The mean read depth of these Holocene consensus sequences ranged from 14.1 to 165.7×. The consensus sequences for the four Holocene and one historic Caucasian wisent have been deposited in GenBank with accession numbers KX553930–KX553934. The mitogenomic sequence from the remaining modern and historic samples has been submitted to GenBank with accession number KY055664 (supplementary table S1, Supplementary Material online).

We assessed phylogenetic relationships among wisent mitochondrial haplotypes, as well as their placement within the wider bovine (tribe Bovini) tree. Wisent sequences were aligned with those of 12 other bovin taxa, including the extinct steppe bison (Bison priscus) and aurochs (Bos primigenius) (supplementary table S1, Supplementary Material online). We excluded two previously published wisent mitochondrial genomes from the phylogenetic analysis (GenBank accessions: HQ223450, HM045017/NC_014044), as these sequences were considered problematic. Specifically, HQ223450 has multiple insertions totaling 9 bp in the ND4 coding region, and HM045017/NC_014044 has multiple indels and point mutations concentrated in the large rRNA and ND3 coding regions. Sequence alignment, partitioning, model testing, and phylogenetic and associated statistical support methods followed the ordinal-level analyses of Heintzman et al. (2015), except that we used the B. bison reference mitochondrial genome (NC_012346) for partitioning. We selected the following models of molecular evolution for the six partitions: GTR + I+G (CP1, 3803 bp; rRNAs, 2541 bp), GTR + G (CP3, 3803 bp), HKY + I+G (CP2, 3803 bp; tRNAs, 1526 bp; control region, 927 bp). We used the saola (Pseudoryx nghetinhensis; supplementary table S1, Supplementary Material online) as outgroup in both the maximum likelihood and Bayesian analyses, following Bibi (2013).

D Statistic Tests

The D statistic involves four genomes: a genome from each of two sister populations (P1 and P2), a genome from a third population as the potential source of introgression (P3), and an outgroup genome (O) to identify the ancestral state (identified as the A allele). We identified variable positions at which P3 possessed the derived allele (B) and presence of the derived allele is variable among P1 and P2, leading to two possible patterns: either ABBA or BABA. Under the scenario of incomplete lineage sorting without geneflow these patterns should occur with equal frequency and the expected D value will be zero. An excess of ABBA or BABA patterns is interpreted as evidence of admixture. However, it might also arise from nonrandom mating in the ancestral population due to population structure (Eriksson and Manica 2012). To determine the ancestral state we used the water buffalo genome. In all tests involving data mapped to the zebu genome, we took into consideration the autosomes only. We performed a total of 105 comparisons considering all possible combinations of wisent, all wisent with either domestic cattle or aurochs as candidate admixer, and domestic cattle and aurochs with all wisent as candidate admixer. These results are reported in supplementary table S4, Supplementary Material online. $f^$ test (Green et al. 2010; Durand et al. 2011) was used to estimate the proportion of the genome derived from admixture. This test requires two individuals of the candidate introgressor species that are not themselves admixed. For our datasets this was possible only for admixture involving the cattle/aurochs lineage. For both D statistic and $f^$ test, significance was assessed using a weighted block jackknife using 1 Mb blocks (Green et al. 2010; Durand et al. 2011). The weighted block jackknife tests if admixture signals are uniform across the whole genome and therefore reflect the same population history. By removing one at a time blocks of adjacent sites (larger than the extent of linkage disequilibrium) and computing the variance of the D statistic or $f^$ values over the entire genome M times leaving each block of the genome in turn, and then multiplying by M and taking the square root we generated the standard error. The number of standard errors by which D or $f^$ differs from zero is the Z score. The results with Z scores greater than 3 in absolute value were qualified as statistically significant (Green et al. 2010).

Nuclear Genome Phylogenetic Tests

The aligned pseudohaploid sequences, generated by mapping reads to the Bos indicus reference, were divided into non-overlapping blocks of 1 Mb. If each of the five taxa contained no more than 50% gaps within a window, the sequence data were recoded into binary characters to only score transversions (Rs: 0, Ys: 1), otherwise the window was recorded as having insufficient data. A Maximum Likelihood phylogeny under the BINGAMMA model and with the water buffalo as outgroup was then computed for each alignment with sufficient data using RaxML (Stamatakis 2014). The topology of each phylogeny was evaluated using a custom Perl script that made use of the ETE3 software (Huerta-Cepas et al. 2016). The lengths of admixed genomic regions were estimated by counting the number of consecutive 1 Mb blocks returning the respective tree topology. Due to the presence of blocks with insufficient data, these measurements are likely to be underestimates. Evaluation of the lengths of genomic regions was conducted using the empirical cumulative distribution function in R, with default parameters.

Supplementary Material

Supplementary data are available at Molecular Biology and Evolution online.

Author Contributions

K.W., J.M.S., A.B., and M.H. conceived the study; A.B., K.W., J.L.A.P., and P.D.H. designed laboratory experiments; K.W., U.T., G.X., P.D.H., and B.S. performed lab work; K.W., S.H., A.B., J.A.C., J.L.A.P., and P.D.H. coordinated data analysis; K.W., S.H., A.B., and P.D.H. performed data analysis; K.W., A.B., and M.H. coordinated writing of the manuscript; K.W., J.M.S., and M.H. obtained funding; G.B., A.N.B., J.J.C., R.D., N.M., H.O., M.T., S.T.T., J.M.W., and W.Ż. provided samples. All authors read, gave comments and helped revise the final version of the manuscript.

Acknowledgments

We would like to thank Andrzej Imiołczyk and his colleagues from the Specimen Preparation Laboratory of the Upper Silesian Museum in Bytom. We thank Andrew Fields and Joshua Kapp for technical assistance. This work was supported by National Science Center (NCN) grant 2011/03/N/NZ8/00086 to KW and ERC consolidator grant GeneFlow 310763 to MH. The authors declare that they have no competing financial interests.

References