Зоологический журнал, 2019, T. 98, № 10, стр. 1091-1099

MATERIALS AND METHODS

We extracted ancient DNA from 29 bison samples (collection of the “National Alliance of Shidlovskiy “Ice Age”) from five locations (fig. 1 , tabl. 1) following the protocol described by (Dabney et al., 2013), with in-house modifications to reduce contamination (Korlević et al., 2015). We then prepared double-stranded DNA sequencing libraries from these extracts following (Meyer, Kircher, 2010). We cleaned the resulting libraries with SPRI beads in 18% PEG-8000 solution and sequenced each on an Illumina MiSeq machine using the v3 chemistry and a paired-end 2 × 75 cycle strategy. Most libraries were barcoded using only one index, although some (F-1267–F-1297) were dual-indexed to avoid index swaps and incorrect read identification (Costello et al., 2018). Libraries were pooled so as to recover at least 500 000 reads from each library. To avoid contamination, all protocols preceding indexing PCR were executed in the PCR-free sterile ancient DNA laboratory facility at the Paleogenomics lab, University of California Santa Cruz. After shotgun sequencing, we then assessed the quality and authenticity of each DNA extract to learn (1) endogenous content (proportion of DNA assigned to Bison); (2) library complexity (diversity of molecules recovered during library preparation); (3) average DNA fragment length (expected to be less than 100bp due to DNA fragmentation in samples that are older than 5000 years); and (4) deamination profile (cytosine damage reflecting DNA degradation, fig. 2 ).

Using SeqPrep (https://github.com/jstjohn/SeqPrep), we removed adapters and then merged reads that overlapped by at least 15bp and discarded reads shorter than 28bp to allow confident mapping. We then mapped all non-discarded reads to the bison nuclear genome UMD 1.0 (GenBank ID GCF_ 000754665.1) using BWA ‘aln’ with the seeding option disabled and a minimum mapping Phred quality score of Q>20 (Li, Durbin, 2009). We then checked the authenticity of the aDNA reads (exploring cytosine deamination patterns) and the fragment size distribution using MapDamage2 (Jónsson et al., 2013).

After determining whether each extract was sufficiently well preserved to proceed, we enriched each library for bison mitochondrial DNA using an in-solution hybridization capture protocol with bison mitogenome RNA baits (Arbor Biosciences, Ann Arbor, USA, previously MYcroarray), following the MYbaits recommended protocol versions 2 and 3, except that the hybridisation step was extended to 36 hours. The enriched libraries were cleaned and pooled following the strategy described above, and sequenced on the Illumina MiSeq.

To process sequenced bait-captured libraries we followed the same protocol as described above, but using B. bonasus NC_014044 mitogenome as a reference for mapping (Zeyland et al., 2012). Finally, we assembled complete mitogenomes using the Mapping Iterative Assembler MIA (Briggs et al., 2009), as described in (Heintzman et al., 2016).

We aligned the 16 newly assembled Bison mitochondrial genomes to a set of 29 previously published extinct and extant Bison species and other bovids using Muscle (Edgar, 2004). For phylogenetic reconstruction, we split the alignments into four partitions and identified the most appropriate models of nucleotide substitution for each partition using jModeltest2 (Darriba et al., 2012). The following setup was used according to results of the test: HKY + G for protein coding sequences, HKY + I + G for D-loop, GTR+G for rRNAs, and HKY + G for tRNAs.

Finally, we performed phylogenetic inference on the partitioned dataset using RAxML (Stamatakis, 2014) and MrBayes (Ronquist et al., 2012) as described in (Maschenko et al., 2017), but with the nucleotide models listed above.

RESULTS

We extracted aDNA and prepared next-generation sequencing libraries from four new specimens collected at the mouth of the Rauchua River (North-Western Chukotka), eight specimens from the Middle Indigirka River (North Yakutia), and 14 specimens from the Bilibinsky region, North-Western Chukotka: Maly Anuy River; Ostrovnovskaya tundra – to the West of the mouth of the Rauchua River; the East Siberian Sea coast – 200 km to the West of the Rauchua River mouth (fig. 1 , tabl. 1). In addition, we also re-extracted DNA from three new F-3246 sub-samples (the original Rauchua lineage individual). Most of the bison samples listed in table 1 were bones, with the exclusion of two samples of horn sheath that are discussed below. Assembled mitochondrial genomes of well-preserved specimens were deposited in GenBank (accession numbers MN049905-MN049913).

As expected, aDNA preservation dramatically varied between samples (tabl. 1). DNA preservation is often reported as a fraction of sequenced reads that map to a conspecific or closely related genome (endogenous proportion of DNA) in relation to the fraction of all other (exogenous, i.e. contaminant) reads that do not map to the host genome. Sources of exogenous DNA include co-extracted co-preserved DNA from plants, animals, and microbes that invade the sample after death. One of the subsamples of the the original Rauchua River sample, subsample F-3246/12, was the best preserved among the three tested, with 55% of recovered DNA mapping to the Bison nuclear genome. Bones sampled from the Middle Indigirka River were the least well preserved. From this set of samples, F-1291 and F-1294 (fig. 3 , tabl. 1) yielded about 1% endogenous DNA. Interestingly, both of these are horn sheath tissue, rather than bone. Samples from Maly Anuy River, Ostrovnovskaya tundra, and East Siberian Sea coast show very good DNA preservation (tabl. 1) based on endogenous bison DNA.

The Rauchua bison lineage was an unexpected finding. We decided that, as part of our undertaking to identify additional samples belonging to this lineage, we would return to the original sample and replicate these results, using samples that were sent and processed entirely independently of the original sample. We extracted new DNA from this specimen using three different bones (F-3246/9, F-3246/12, F-3246/13): right humerus, right fore autopodia, and left ulna. All three DNA extracts were successfully enriched, producing high coverage mitogenomes (tabl. 1). None of these mitogenomes showed discrepancies with the original published Rauchua haplotype, suggesting that the Rauchua bison lineage as reported is valid and is not an artifact of mitochondrial genome assembly, contamination, or other sources of error.

Despite poor preservation of bison DNA in half of the newly processed bison samples, we were able to enrich most DNA libraries for sufficient quantities of mitochondrial DNA to assemble mitogenomes to >3x coverage for 16 samples, which we then used for phylogenetic analysis. Both Bayesian and Maximum Likelihood approaches produced consistent tree topologies (fig. 3 ). All newly assembled bison mitochondrial lineages fall within the known diversity of steppe bison mitochondrial genomes. None of the new samples clustered with the original Rauchua lineage, including a new sample from the Rauchua River mouth (F-4171), the same locality where the original Rauchua bison was found.

Fig. 1.

The map of Eurasia with sample locations marked with a circle: Middle Indigirka River (North-Eastern Yakutia), Bilibinsky region (North-Western Chukotka), Rauchua River (North-Western Chukotka).

Fig. 2.

Ancient DNA degradation profile of F-1291 from Middle Indigirka (panels a, b) and F-3246/9, the new subsample of the original Rauchua bison (panel c, d). Deamination profile is shown in the panels a and c. Solid black line – frequency of C to T substitutions at 5` ends of DNA molecules. Numbers on the x-axis are relative nucleotide positions on each DNA fragment mapped to the Bison reference nuclear genome. Dotted line – frequency of G to A substitutions. Fragment length distribution is shown on panels b and d.

Fig. 3.

Maximum likelihood phylogenetic tree reconstructed on complete mitochondrial genomes of various bovids, including extinct and extant Bison species. Support for major divergences is marked below or above the branches: bootstrap calculated with RAxML, and posterior probability (in percent, calculated with MrBayes). Both numbers are shown in cases where corresponding node was supported by both phylogeny reconstruction methods. Divergence of Rauchua lineage is shown with a circle. The new specimen from Rauchua locality is marked with a star. Mitochondrial genomes assembled in this study are in bold. All previously published sequences are shown with corresponding GenBank ID numbers. Numbers at the end of sequence names are median IntCal13-calibrated radiocarbon ages (Reimer et al., 2013).

Two samples, F-1297 (from the Middle Indigirka River) and F-4172 (from the Rauchua River mouth), have very little DNA preserved, and resulted in an average mitochondrial coverage of around 1x. While we were unable to include them in our phylogenetic reconstructions, we were nonetheless able to genotype these individuals from the portions of the mitochondrial genome that were recovered. Using the Rauchua mitochondrial haplotype and a panel of extinct B. priscus as a reference, we identified nucleotide substitutions in the mitochondrial genome that were distinct from other bison and unique to the Rauchua lineage. Most of these substitutions are located in the cytochrome oxidase unit 1, the gene that has been used as a barcode for confident identification of various taxa (Luo et al., 2011). We manually inspected those fragments of the poorly preserved F-1297 and F-4172 mitochondrial genomes that were covered by a minimum of three reads (to ensure that any differences were supported by independent reads and not therefore due to DNA damage or sequencing error) looking for informative substitutions. In both samples we identified 35 informative nucleotide sites that were both covered by at least three reads and span nucleotide positions that would have Rauchua-specific substitutions. Neither F-1297 nor F-4172 had the Rauchua allele at any of these positions, and instead both had that of the main lineage of Siberian bison. These results indicate that both F-1297 and F-4172 are part of the main lineage of Siberian B. priscus, and not the Rauchua lineage. We would like to mention, however, that only F-1297 was dual indexed to avoid incorrect read assignment. Base calling of the single-indexed F-4172 library should be interpreted with caution due to a possibility of index hopping (Costello et al., 2018).

DISCUSSION

The results of this study demonstrate that the Rauchua lineage is even more enigmatic than previously thought. Despite extensive sampling in the northern Russian Far East, covering most of the territory where Eurasian steppe bison were present during the Pleistocene–Holocene transition and the Rauchua locality itself, we did not find any new specimens belonging to the Rauchua bison lineage (fig. 3 ). Instead, our results show that the majority of western Beringian Bison priscus cluster together into a clade within the diversity of Holarctic Bison priscus mitochondrial genomes, as described previously (Shapiro et al., 2004; Heintzman et al., 2016; Froese et al., 2017).

We currently do not know why the population to which the Rauchua bison belonged is so rare in the paleontological record. It is notable that the Rauchua bison itself has a relatively recent radiocarbon date of only 9497 ± 92 years BP. The end of the Pleistocene was marked by megafaunal extinctions across the globe, with most mainland extinctions in Siberia occuring by 14-10 thousand years ago (Stuart et al., 2004; Koch, Barnosky, 2006). The Rauchua bison therefore postdates these extinctions, and is one of the youngest bison thus far dated from the Russian Far East (Van Geel et al., 2014; Markova et al., 2015). It is possible that the Rauchua bison lineage persisted in more southerly locations during the Pleistocene, where bones are either not preserved or have yet to be recovered. The Rauchua lineage may have moved northward as the climate improved after the transition into the Holocene, and established populations, albeit briefly, in the region that was previously occupied by the main lineage of Siberian steppe bison.

Bison is among the most well represented taxon in Late Pleistocene fossil deposits and among the taxa best represented by ancient mitochondrial DNA data (Shapiro et al., 2004; Heintzman et al., 2016; Soubrier et al., 2016; Palacio et al., 2017; Froese et al., 2017; Zazula et al., 2017). This research has revealed several surprises, including the presence of other cryptic bison lineages. For example, (Soubrier et al., 2016) reported a genetically distinct but morphologically indistinguishable bison lineage that may have coexisted with Bison priscus in the Southern Ural region during the Late Pleistocene. This cryptic lineage, which has been referred to as Bb1 or Bison schoetensacki, has mitochondrial data that is more similar to European bison, Bison bonasus, than it is to the steppe bison Bison priscus (Soubrier et al., 2016; Palacio et al., 2017). The discovery and genetic characterization of both Clade X and the Rauchua bison lineage highlight the extensive genetic diversity of this geographically widespread taxon, which further paleontological and genetic research will be necessary to explore. In particular, additional bison samples from southern regions of Eurasia are necessary to shed light on evolution of the Bison genus, and on the divergent Rauchua lineage specifically.

Despite an extremely rich paleontological record, hypotheses about Bison evolution and taxonomy have been conflicting for a long time, because of the extensive morphological diversity within the genus (Grange et al., 2018). While mitochondrial DNA data can provide new insights into the evolution of this particular group of taxa, mitochondria represent only a single, maternally inherited locus. The evolutionary history revealed by analysis of mitochondrial DNA tells, therefore, only a part of the evolutionary story of the species in question. To fully resolve the evolutionary history of bison, and to fully understand the relationship between the Rauchua lineage and other Siberian steppe bison, analysis of nuclear DNA will be necessary. Fortunately, some of the samples presented in this study are sufficiently well preserved to attempt this work.

The results of our study confirm the existence of the deeply divergent Rauchua clade and shed new light on the evolutionary history of bison during the Pleistocene to Holocene transition in Siberia. Future analysis of samples processed for this study will reveal when the Rauchua bison lineage diverged from that of the other bison and whether and to what extent the Rauchua bison population admixed with other Late Pleistocene bovids.