Bison phylogeography constrains dispersal and viability of the Ice Free Corridor in western Canada

To infer the chronology of the corridor linking Beringia and interior North America, we generated radiocarbon dates from 78 North American bison fossils, 49 of which were recovered from the corridor region (Table 1 and Dataset S1). Sites included Clover Bar and Charlie Lake Cave in western Canada (Fig. 1), both of which previously yielded dates in the time frame of interest (11,500–13,500 cal y BP) (10). The Charlie Lake Cave fossils were also associated with archaeological materials (28). We generated mitochondrial haplotypes from 45 of these dated fossils, including from 18 of the 22 from Clover Bar and Charlie Lake Cave that fall within the time frame of interest. We then used these and previously published haplotypes to estimate a mitochondrial genealogy for a total of 192 late Pleistocene, Holocene, and present-day North American bison, including 37 from the corridor region and within the time frame of interest (Fig. 1, Fig. S1, Table 1, and Dataset S1). To facilitate discussion, we divide the mitochondrial genealogy into two major clades, clade 1 and clade 2, which are based on the most deeply diverging lineages within the mitochondrial tree. We also highlight several well-supported subclades within these two major clades.

Our recovered genealogy is similar to previously published mitochondrial genealogies for bison (11, 29), in which the most striking feature of the tree is the clustering of all present-day bison into clade 1a, with a maternal common ancestor that postdates the LGM (Fig. 2A). The increased density of new postglacial bison included in our analysis refines the pattern of extinction of clade 2, within which the latest surviving lineages tend to cluster together both phylogenetically and geographically (Fig. 2A). Clade 2a comprises a late-surviving population of bison that appears to be geographically restricted to southern Yukon and interior Alaska and includes a fossil that dates to as recently as ∼325–490 cal y BP. Clade 2b is geographically isolated to the region of the postglacial corridor and includes a skull recovered from Banff National Park that dates to ∼305–430 cal y BP (Dataset S1). This increased phylogeographic clustering during postglacial times, compared with before the LGM, probably reflects a trend toward physical isolation of bison populations as open habitats were largely replaced by spruce forest and increasing paludification across northwestern Canada (30).

The pattern of bison dispersal into the postglacial corridor provides increased resolution of both the timing of establishment of an ecosystem that can support grazing herbivores and the capacity for the animals to traverse this region. Previously published dates for a horse from Vauxhall, Alberta (10), and a bison from Tsiigehtchic, Yukon (29), indicate a postglacial corridor route had begun to form at both the northern and southern ends by at least 13,500 cal y BP. Bison from clade 1a were present in the southern corridor near present-day Calgary by 13,430–12,875 cal y BP, and bison from clade 2b were in northeastern British Columbia by 13,175–13,045 cal y BP. The postglacial corridor was fully open for dispersals by 13,000 cal y BP, when multiple overlapping radiocarbon dates belonging to bison from all three late-surviving clades (Fig. 2; clades 1a, 2a, and 2b) are found at the Clover Bar site (Fig. 1 and Table 1). The fact that bison from all three late-surviving clades are present during this same interval is crucial; although we cannot determine with confidence the direction of dispersal of bison from clade 2b, our analyses indicate a northern origin for bison from clade 2a (Fig. 2B). Bison from clades 1a and 2b were present in stratigraphic subzones IIb–d at Charlie Lake Cave for ∼1,000 y before the end of the Pleistocene, from ∼12,500 to 11,500 cal y BP (28) (Fig. 1 and Table 1). These corridor bison were part of a diverse megafaunal community that included American lion (Panthera leo atrox), horse, western camel (Camelops hesternus), caribou (Rangifer tarandus), tundra muskoxen (Ovibos moschatus), helmeted musk oxen (Bootherium bombifrons), and mammoth (1, 31), further supporting the notion that the area was productive habitat at the end of the Pleistocene. Based on present understanding of the biogeography of many of these mammal taxa, there is a general hypothesis that, apart from the clade 2a bison, most of the taxa that recolonized the corridor region came from the south, including horses, lions, camels, and muskoxen (10).

Data presented here and elsewhere demonstrate the critical role that continental ice sheets and the postglacial corridor played in biogeographic patterning among species and populations. For many mammalian taxa, paleontological and genetic data reveal distinct separation of populations north and south of the ice. Intriguingly, despite our conclusion that bison disperse into the postglacial corridor from both the north and south, we find only limited evidence of dispersal beyond the region of the corridor. In fact, the only evidence of dispersal completely through the corridor is the occurrence of a bison from clade 1a at Liard River, Northwest Territories, at 12,390–12,035 cal y BP, which corresponds to a northward dispersal (Fig. 2B and Table 1). Similarly to bison, several other taxa with Eurasian origins, such as caribou and American lion, probably traversed the corridor region before the LGM and established populations in the southern interior of North America (10). However, the responses of these species to deglaciation and environmental change at the end of the Pleistocene varied. For example, although bison from the south dispersed as far north as the Liard River, there is no mitochondrial evidence of further northward expansion into Alaska and Yukon during the Holocene. In contrast, mitochondrial data suggest that all present-day wolves (Canis lupus) in the North American subarctic and arctic are descended from a population that was south of the continental ice sheets during the LGM and that dispersed northward during the postglacial period (32).

These results have two key implications for the role of the postglacial corridor as a pathway for biotic exchange between Beringia and interior North America. First, the opening of the postglacial corridor may have favored south to north, rather than north to south dispersal. Detailed biome reconstructions indicate that southern and central portions of the deglaciating corridor in this time range featured potentially more productive grasslands, open spruce woodlands, and boreal parkland, whereas northern portions of the corridor were marked by alpine, herb, and shrub tundras (30). Southward dispersal may also have been limited for biological reasons, for example, if southern bison were better adapted than northern bison to the expanding grasslands within the corridor region (27). Second, the interval of time during which it is feasible to transect the corridor may have been limited. For bison, the barrier to further dispersal may have been the relatively quick replacement of grasslands at the northern end of the corridor by increasingly closed spruce forests, which are difficult for grazing herbivores such as bison to transect (26, 33).

The expansion of bison into the corridor region provides proxy evidence for when this route was viable for human populations and, in doing so, allows further refinement of New World human settlement scenarios. Human genetic and archaeological evidence indicate that eastern Beringia and parts of the Americas well south of the ice sheets were populated by 14,000 cal y BP, suggesting that migration out of Beringia probably began more than 15,000 cal y BP ago (15, 34–36). Our chronology for the opening of the postglacial corridor indicates that a fully habitable corridor connected Beringia and interior North America by ∼13,000 cal y BP. This timing precludes the postglacial corridor as a southward route for initial human dispersal into the Americas, the corollary being that the first indigenous peoples leaving Beringia probably took a coastal route or potentially moved through western North America before glacial coalescence (37, 38).

We find that a bison belonging to the northern clade (2a) reached the Edmonton area by 13,000 cal y BP. It is therefore possible that established northern human populations also reached the central corridor by this time. Evidence from the archaeological record supports this hypothesis. For example, Alaskan archaeological sites including Swan Point, Mead, Broken Mammoth, Tuluaq, and Dry Creek, which were occupied from ∼14,000 to 11,500 cal y BP, feature a variety of projectile technologies, sometimes associated with microblade industries (39). Similar microblade technologies are present at Vermilion Lakes (Banff National Park) and Charlie Lake Cave by ∼11,500 cal y BP (28, 40, 41). In addition, human genetic data from Upward Sun River, Alaska, show founding New World mitochondrial haplotypes B2 and C1b in Alaska at ∼11,500 cal y BP. Small, isolated groups of people may therefore have continued to disperse from Beringia to interior North America well after the corridor region opened (14, 16, 41, 42).

Our bison data also suggest that biotic conditions favored northward rather than southward movements through the corridor, paralleling archaeological data involving fluted point technology (13). Bison clade 1a, which originates south of the ice sheets during the period of coalescence, predominates in our corridor sample, with one instance occurring as far north as the Liard River by ∼12,200 cal y BP. The oldest recognizable Clovis complex sites in North America are estimated to range from 13,000 to 12,600 cal y BP (8, 36), which slightly postdates our chronology for the opening of the southern end of the corridor. The Anzick child burial in Montana that dates to ∼12,600 cal y BP (43), and slightly earlier evidence for human hunting of western camels and horses at Wally’s Beach, 420 km north of the Anzick site (17, 43), document the presence of people at the southern end of the corridor. Fluted point density maps indicate that this technology diminishes in frequency from the south to the north in the corridor region, consistent with the northward spread of this technology (1, 12, 42). By 12,500 cal y BP, fluted points are present at Charlie Lake Cave in British Columbia (along with clades 1a and 2b bison) and in sites in Alaska (13, 28).

Collagen was extracted from bone and tooth samples using ultrafiltration methods outlined in ref. 44 and radiocarbon dated using accelerator mass spectrometry (AMS) at the Keck laboratory–University of California (UC) Irvine (UCIAMS), the Center for AMS (CAMS) at Lawrence Livermore National laboratory, or Beta Analytics (Table 1 and Dataset S1). AMS dates were calibrated using the IntCal13 calibration curve (45) in OxCal v4.2 (https://c14.arch.ox.ac.uk/oxcal/OxCal.html) and are reported with 1 standard deviation (SD). For samples that were redated at UCIAMS (SFU 1848, SFU 1849, SFU 3429, and SFU 15004; Dataset S1), we used only the new dates for the calibration. All dates reported in the text are in calendar years before present, unless otherwise noted.

To facilitate comparison with previously published data from bison, we isolated ∼600 bp of the hypervariable portion of the mitochondrial control region (CR) from 45 Canadian bison bone and tooth samples, dated to the late Pleistocene and Holocene (Table 1 and Dataset S1). We performed DNA extraction, library preparation, and PCR setup in dedicated ancient DNA facilities at the Pennsylvania State University (PSU) and UC Santa Cruz (UCSC) that were physically isolated from modern molecular biology research. Depending on the sample, we extracted DNA from 100 to 250 mg bone powder using one of two methods that are highly optimized for the recovery of ancient DNA molecules (46, 47), performing one negative extraction control for every five to eight processed samples.

To generate CR sequence data, we used a mixture of four approaches: (i) direct Sanger sequencing of PCR products, (ii) Sanger sequencing of cloned PCR products, (iii) Illumina amplicon sequencing of PCR products, and (iv) mitochondrial target enrichment followed by Illumina sequencing. For approaches i–iii, we amplified the target CR fragment either in a single PCR amplification or as a series of overlapping fragments, depending on the preservation of the sample (primer combinations are provided in Table S1) (11). We performed PCR in 25-μL reactions with the following components: 20 μg rabbit serum albumin, 0.25 mM dNTPs, 1× High Fidelity buffer, 1 U Platinum Taq High Fidelity (Life Technologies), 2.4 mM MgSO₄, 0.4 μM of each primer, and 1 μL DNA extract, with the following cycling conditions: 12 min at 94 °C, 30 s at 94 °C, 45 s at variable annealing temperature (Table S1), 45 s at 68 °C, and 1 min at 68 °C, with the middle three steps repeated for 50 cycles. We cleaned PCR products using either Millipore μ96 plates or Sera-Mag SpeedBeads (ThermoScientific) in 18% (wt/vol) PEG-8000, the latter of which followed the bead-based reaction clean-up protocol of ref. 48.

For approaches i–iii above, we then used one or more of these approaches to assess the accuracy of the resulting PCR amplicon sequences: (i) bidirectional direct sequencing on AB3730xl genetic analyzers at the PSU Genomics Core Facility or the UC Berkeley DNA Sequencing Facility, using BigDye v3.1 chemistry (Life Technologies); (ii) cloning PCR products using the TOPO-TA cloning kit (Life Technologies) according to the manufacturer’s instructions, followed by selection and PCR amplification of six to eight colonies following ref. 11; and/or (iii) pooled sequencing of barcoded PCR products using the Illumina MiSeq platform, in which PCR products derived from the same sample were pooled in equimolar ratios and turned into Illumina DNA libraries using ref. 49 with modifications from ref. 50. These indexed libraries were then pooled in equimolar ratios and sequenced on the Illumina MiSeq platform using v2 150-bp paired-end chemistry, following the manufacturer’s instructions.

For the fourth data generation approach described above, we constructed Illumina DNA libraries as above directly from the DNA extract. We then enriched these libraries for the whole bison mitochondrial genome using biotinylated RNA baits (MYbaits v2; MYcroarray), following the manufacturer’s instructions. These enriched libraries were then sequenced on the MiSeq as described above, but using v3 75-bp paired-end chemistry.

For data generation approaches i and ii, we assessed sequence quality manually and called consensus sequences using Lasergene v9 (DNASTAR) or Geneious v6.1.6 (Biomatters). For approaches iii and iv, we binned the short-read data by index, and then removed adapters and merged paired-end reads using SeqPrep (https://github.com/jstjohn/SeqPrep). For approach iii, we removed primer sequences from merged and remaining unmerged reads using in-house scripts and mapped each read to both Bison bison (GenBank: NC_012346) and B. priscus (AY748705) CR reference sequences using BWA v0.6.1 (51), resulting in two datasets per sample. We called consensus sequences using Geneious with the base agreement threshold set to 75% and the minimum coverage set to 50×. For approach iv, we aligned merged reads to the reference B. bison mitochondrial genome sequence using an iterative assembler MIA (52). We called bases that had >3 times coverage and >67% agreement. For analysis, we then trimmed the resulting consensus sequence to the ∼600-bp target. We combined data from all four approaches to create robust consensus sequences for each sample. If we identified conflicts between sequences resulting from the different approaches, we either performed additional analyses or coded conflicting sites as ambiguous. Unsequenced regions are considered missing data.

We aligned the 45 new CR sequences to a dataset of 147 previously published, radiocarbon dated or present-day, North American bison CR sequences (Dataset S1), using Se-Al (Sequence Alignment Editor; v2.0a11). We performed a Bayesian phylogeographic analysis using BEAST v1.8.3 (53). We assumed a generalized time-reversible evolutionary model, with gamma distributed rate variation and a proportion of invariable sites (GTR+G+I), a strict molecular clock with a rate calibrated using the median calendar ages of each radiocarbon dated specimen, and the flexible skygrid model of the coalescent process (54). To infer the timing and directionality of movement between Beringia and interior North America, we adopted the discrete phylogeographic model described in ref. 55, assigning each sample to either the North or South population based on whether the sample originated to the north or south of 60° N, which archaeological data and simulation studies indicate was the final barrier to a corridor (56). To simplify interpretation of the results, the two present-day bison that were sampled at locations just north of this cutoff (at 60° and 61.4°) were classified as southern.

To learn about ancestral movements between north and south, we estimated the posterior distribution of the time of north-to-south and south-to-north movements using Markov Jumps (57). This technique exploits dynamic programming and tricks from numerical analysis to efficiently compute the expected number and timing of specific transitions within a continuous-time Markov chain (CTMC) that conditions only on the directly observed end states at the tips of phylogenetic tree on which the CTMC acts. The resulting posterior distribution of transition times naturally incorporates uncertainty in the tree and estimated CTMC rates. We ran two MCMC chains for 50 million states each, sampling the posterior states of all model parameters and trees every 5,000 states. We discarded the first 10% of samples from each run as burn-in, combined the remainder using LogCombiner v1.8.3, and computed posterior means, posterior SDs, and 95% highest posterior density intervals in Tracer v1.6. We summarized the combined set of posterior trees and identified the maximum clade credibility (MCC) tree using TreeAnnotator v1.8.3, which we visualized using Figtree v1.4.2.

The input BEAST file is available as Dataset S2. Novel control region mitochondrial DNA sequences have been deposited in GenBank, with accession numbers KU705765–KU705809. All fossil specimens used in this study are curated in the repositories listed in Dataset S1.