13

Bison Phylogeography Constrains Dispersal and Viability of the Ice-Free Corridor in Western Canada

PETER D. HEINTZMAN,* DUANE FROESE^†§§§JOHN W. IVES,^‡ANDRÉ E. R. SOARES,* GRANT D. ZAZULA,^§BRANDON LETTS,^||THOMAS D. ANDREWS,^#JONATHAN C. DRIVER,** ELIZABETH HALL,^§P. GREGORY HARE,^††CHRISTOPHER N. JASS,^‡‡GLEN MACKAY,^#JOHN R. SOUTHON,^§§MATHIAS STILLER,^{|| ||}ROBIN WOYWITKA,^†MARC A. SUCHARD,^##***^††† AND BETH SHAPIRO*^{‡‡‡§§§}

__________________

^* Department of Ecology and Evolutionary Biology, University of California, Santa Cruz, CA 95064; ^†Department of Earth and Atmospheric Sciences, University of Alberta, Edmonton, AB T6G 2E3, Canada; ^‡Institute of Prairie Archaeology, Department of Anthropology, University of Alberta, Edmonton, AB T6G 2H4, Canada; ^§Yukon Palaeontology Program, Government of Yukon, Whitehorse, YT Y1A 2C6, Canada; ^||Department of Biology, The Pennsylvania State University, University Park, PA 16802; ^#Prince of Wales Northern Heritage Centre, Government of the Northwest Territories, Yellowknife, NT X1A 2L9, Canada; **Department of Archaeology, Simon Fraser University, Burnaby, BC V5A 1S6, Canada; ^††Yukon Archaeology Program, Government of Yukon, Whitehorse, YT Y1A 2C6, Canada; ^‡‡Royal Alberta Museum, Edmonton, AB T5N 0M6, Canada; ^§§Keck-CCAMS Group, Earth System Science Department, University of California, Irvine, CA 92697; ^{|| ||}Department of Translational Skin Cancer Research, German Consortium for Translational Cancer Research, D-45141 Essen, Germany; ^##Department of Biomathematics, David Geffen School of Medicine at University of California, Los Angeles, CA 90095; ***Department of Human Genetics, David Geffen School of Medicine at University of California, Los Angeles, CA 90095; ^†††Department of Biostatistics, University of California, Los Angeles Fielding School of Public Health, Los Angeles, CA 90095; and ^‡‡‡University of California Santa Cruz Genomics Institute, University of California, Santa Cruz, CA 95064. ^§§§To whom correspondence may be addressed. Email: duane.froese@ualberta.ca or bashapir@ucsc.edu.

A central question in New World biogeography and archaeology has been the role of an “Ice-Free Corridor” along the eastern slopes of the Rocky Mountains in facilitating biotic exchange between Beringia—unglaciated Alaska and Yukon—and southern interior parts of the Americas (Johnston, 1933; Mandryk, 1996; Mandryk et al., 2001; Ives et al., 2013). Of central importance is the potential role for the corridor between the Laurentide and Cordilleran ice sheets in defining the pattern and timing of the initial entry of humans into the New World (Fig. 13.1). Early geological models from western Canada suggested that the ice sheets did not coalesce during the Last Glacial Maximum (LGM) (Stalker, 1977) and implied that the corridor was a viable route for dispersal throughout the late Pleistocene (~125,000–11,500 years ago). This view gained support from numerous, misleadingly “old” radiocarbon ages from bulk deposits in the corridor region (MacDonald et al., 1987). The suggestion of an open and viable corridor throughout the LGM gave rise to its prominence as the most likely pathway taken by the first people to colonize the Americas from Beringia (Fladmark, 1986; Mandryk, 1996; Mandryk et al., 2001). Eventually, this notion of a corridor became tightly linked with the “Clovis First” archaeological model, in which widespread fluted projectile point technology, often found in association with remains

of Pleistocene megafauna, was believed to reflect initial human dispersal into the Americas (Goebel et al., 2008).

More recent geological interpretations indicate that the Laurentide Ice Sheet coalesced with the Cordilleran Ice Sheet along the eastern foothills of the Rocky Mountains by ~21,000 calendar years before present (cal y BP, with present defined as CE 1950) (Jackson et al., 1997). Ice sheet coalescence would have formed an impenetrable barrier to terrestrial dispersals between Beringia to the north and interior North America to the south

that lasted at least until glaciers retreated near the end of the Pleistocene. As the ice sheets receded, a corridor formed that opened a route between Beringia and the continental interior, with its most likely route running through Yukon Territory along the rapidly wasting Cordilleran Ice Sheet or perhaps along the Mackenzie Valley (Ives et al., 2013). A gap in the regional vertebrate radiocarbon date record in Alberta (Burns, 2010) and genetic evidence of a strong barrier to gene flow within large mammal populations to the north and south of the ice sheets during this time interval (Shapiro et al., 2004) both support the geological interpretation of coalescence of the Laurentide and Cordilleran ice sheets.

As geological and paleoenvironmental evidence mounted that the corridor was not available during the LGM, a hypothesized Pacific coastal route emerged as the more likely route to account for the earliest people in southerly regions of the Americas (Mandryk et al., 2001; Amick, 2016). This hypothesis was coupled with the suggestion that, even if an early deglacial corridor existed, the landscape was probably ecologically unproductive and unable to support large mammal or human populations, perhaps for millennia following the retreat of glaciers (Mandryk et al., 2001). The Clovis First model lost favor as sites south of the LGM ice sheets, such as Monte Verde (Chile), the Gault locality (Texas), Manis (Washington), and Paisley (Oregon), provided archaeological evidence that humans were present in the Americas at least a millennium before the appearance of Clovis technology (Goebel et al., 2013). At the same time, genetic analyses suggested that New World founding populations began dispersing from Beringia ~15,000–20,000 cal y BP (Tamm et al., 2007; Achilli et al., 2013; Raghavan et al., 2015).

Despite its importance in North American paleoecology and archaeology, interpretations of the corridor are based on data from a sparsely studied, vast region that was impacted significantly by glacial processes. Questions remain about the precise timing of the postglacial opening of the corridor (Amick, 2016), its biological carrying capacity following deglaciation (Mandryk et al., 2001), and when, if at all, it could have served as a viable route for the movement of people and animals between Beringia and the rest of the Americas. Available archaeological data demonstrate that human populations were present near the southern end of the corridor just before deglaciation (Ives et al., 2013; Waters et al., 2015). Unfortunately, the limited and often perplexing archaeological evidence from the corridor region itself makes it difficult to assess directly when people, technologies, or both traversed the corridor after deglaciation (Ives et al., 2013). More recently, archaeological interest in the corridor has focused on its potential as a route for northward postglacial dispersal, most notably to explain the appearance of fluted projectile point technology in Alaska (Goebel et al., 2013).

BISON PHYLOGEOGRAPHY TO INFER THE CHRONOLOGY OF THE ICE-FREE CORRIDOR

Advances in paleogenetic approaches and DNA sequencing technologies have enabled the recovery of large datasets of populations sampled over time (Shapiro and Hofreiter, 2014). Although DNA preservation and the chronological range of radiocarbon dating have limited ancient DNA studies largely to the last ~50,000 years, this temporal interval spans two events that were undeniably important in shaping the present-day distribution of biodiversity: the transition into and back out of the LGM, which reached its peak ~26,000–19,000 cal y BP (Clark et al., 2009), and the global expansion of human populations (Koch and Barnosky, 2006). Although some genetic signal of both of these events is retained in living populations, paleogenetic data provide temporal snapshots of phylogeographic structure that may not be otherwise detectable (Ramakrishnan and Hadly, 2009). For example, paleogenomic data from ancient Europeans have revealed present-day Sardinians to be more closely related to the ancestral human population in Europe than to any other present-day population, which probably reflects a genetic replacement associated with the expansion of agriculture (Lazaridis et al., 2014; Skoglund et al., 2014). Paleogenetic data have also been used to test hypotheses about the influence of humans on the phylogeography of other species, for example, via domestication (Larson et al., 2007a; Kistler et al., 2014).

Bison (Bison sp.) are an ideal taxon for assessing when the corridor was open and available for biotic dispersals. Although bison taxonomy is complex, with up to several dozen species named based largely on skull and horn morphology (Guthrie, 1990; Wilson, 1996), the first bison to enter North America from Asia during the Pleistocene are generally referred to as the steppe bison, Bison priscus. Bison are one of the most abundant large mammals recovered from faunas within the western interior of North America, although previous studies have provided only limited radiocarbon and mitochondrial data from the corridor region. Bison are one of the first species for which ancient DNA data were used to reconstruct changing patterns of population structure spanning tens of thousands of years at continental scales (Shapiro et al., 2004). In addition, unlike many of their Pleistocene counterparts such as mammoths (Mammuthus), horses (Equus), and lions (Panthera), bison survived the extinction event at the end of the Pleistocene and persist in North America. Bison genomes, however, contain the signature of a genetic bottleneck at the end of the Pleistocene that can be used to test hypotheses about how ice age environmental changes affected their distribution and abundance (Shapiro et al., 2004).

Previous work using radiocarbon and ancient DNA data demonstrated that bison mitochondrial genetic diversity (but not necessarily their morphological diversity) (Wilson, 1996) is phylogeographically par-

titioned across North America and has been since the LGM (Shapiro et al., 2004). Although bison survived the interval of glacial coalescence both north and south of the continental ice sheets, population bottlenecks and barriers to gene flow affected their mitochondrial diversity. By the time the glaciers began to retreat, bison populations that had been isolated to the south of the continental ice sheets were mitochondrially distinct from their contemporary northern counterparts in Beringia. A mitochondrial clade of southern bison, including the two present-day bison subspecies in North America, the plains bison (B. bison bison) and wood bison (B. b. athabascae), shares a common ancestor dating to the period of glacial coalescence ~15,000–22,000 cal y BP (Shapiro et al., 2004). Thus, the identification of bison from this southern clade within and north of the corridor region can be interpreted as reflecting northward dispersal. Likewise, the appearance of bison from a Beringian mitochondrial clade further south in interior North America can be interpreted as southward dispersal.

Here, we report new radiocarbon dates and mitochondrial haplotype data from fossil bison recovered from the corridor region and adjacent areas in Beringia (Fig. 13.1, Table 13.1, and Dataset S1¹) and use these data to assess when the corridor opened and was available as a dispersal route. In addition to resolving the chronology for when the corridor was open, we use a Bayesian phylogeographic approach to infer the timing and directionality of dispersals for ancient bison between the north and south. The specific evolutionary history of bison mitochondria makes this an ideal genetic marker to assess the timing and viability of the corridor and thereby provides insights relevant not only to bison, but also to other species, including humans, that may have used the corridor to disperse between Beringia and interior North America during the late Pleistocene.

RESULTS AND DISCUSSION

A Refined Chronology for the Opening of the Postglacial Corridor

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 13.1 and Dataset S1). Sites included Clover Bar and Charlie Lake Cave in western Canada (Fig. 13.1), both of which previously yielded dates in the time frame of interest (11,500–13,500 cal y BP) (Burns, 2010). The Charlie Lake Cave fossils were also associated with archaeological materi-

__________________

¹ Supporting information for this chapter, which includes Datasets S1 and S2, Fig. S1, and Table S1, is available online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1601077113/-/DCSupplemental.

als (Driver et al., 1996). 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. 13.1, Fig. S1, Table 13.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 (Shapiro et al., 2004; Zazula et al., 2009), 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. 13.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 phylo-genetically and geographically (Fig. 13.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 (Dyke, 2005).

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 (Burns, 2010), and a bison from Tsiigehtchic, Yukon (Zazula et al., 2009), 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. 13.2; clades 1a, 2a, and 2b) are found at the Clover Bar site (Fig. 13.1 and Table 13.1). The fact that bison from all three late-surviving clades are present during this same interval is crucial; although we can-

TABLE 13.1 Bison Specimens Recovered from the Corridor Region During the Time Frame of Corridor Opening (13,500–11,500 cal y BP)

NOTES: Clade numbers correspond to Fig. 13.2A. AB, Alberta; BC, British Columbia; NWT, Northwest Territories; YT, Yukon Territory.

not determine with confidence the direction of dispersal of bison from clade 2b, our analyses indicate a northern origin for bison from clade 2a (Fig. 13.2B). Bison from clades 1a and 2b were present in stratigraphic subzones IIb–d at Charlie Lake Cave for ~1,000 years before the end of the Pleistocene, from ~12,500 to 11,500 cal y BP (Driver et al., 1996) (Fig.

13.1 and Table 13.1). These corridor bison were part of a diverse mega-faunal 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 (Jass et al., 2011; Ives et al., 2013), 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 (Burns, 2010).

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. 13.2B and Table 13.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 (Burns, 2010). 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 (Leonard et al., 2007).

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 (Dyke, 2005). 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 (Guthrie, 1990). 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 (MacDonald and McLeod, 1996; Wilson, 1996).

Consequences of the Postglacial Corridor Chronology for North American Human Prehistory

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 (Gruhn and Bryan, 2011; Reich et al., 2012; Waters and Stafford, 2013; Raghavan et al., 2015). 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 (Madsen, 2004; Erlandson, 2013).

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 (Goebel and Buvit, 2011). Similar microblade technologies are present at Vermilion Lakes (Banff National Park) and Charlie Lake Cave by ~11,500 cal y BP (Fedje et al., 1995; Driver et al., 1996; Tackney et al., 2015. 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 (Tamm et al., 2007; Achilli et al., 2013; Ives, 2015; Tackney et al., 2015).

Our bison data also suggest that biotic conditions favored northward rather than southward movements through the corridor, paralleling archaeological data involving fluted point technology (Goebel et al., 2013). 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 (Goebel et al., 2008; Waters and Stafford, 2013), 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 (Rasmussen et al., 2014), and slightly earlier evidence for human hunting of western camels and horses at Wally’s Beach, 420 km north of the Anzick site (Rasmussen et al., 2014; Waters et al., 2015), 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 (Ives et al., 2013; Ives, 2015; Amick, 2016). 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 (Driver et al., 1996; Goebel et al., 2013).

CONCLUSIONS

Over the last several decades, phylogeographic research has evolved from being largely qualitative in nature to benefitting from a diversity of statistically rigorous approaches to detect evolutionary structure over both geographic space and, with the addition of paleogenetic technologies, over time. Here, we demonstrate the power of this approach to address a longstanding question in New World biogeography, specifically whether and when an ecologically viable corridor linked Beringia and interior North America during the late Pleistocene.

Until now, much of the research concerning the corridor region has focused on whether or not the corridor was the initial pathway for early human inhabitants of the New World. Our data from fossil bison help shift scientific attention to a different first: the LGM coalescence of the ice sheets marked the first time in the Quaternary that the corridor did not exist as a biogeographic pathway. Once the corridor opened, bison populations that had been isolated for millennia, in Beringia and in interior North America, entered the deglaciated region from both its northern and southern ends. Human populations in eastern Beringia and interior North America had similar possibilities to reengage. In the human case, this involved the added cultural complexity that would accompany societal interactions in an era where exotic raw materials and technological ideas began to

circulate widely in the Americas. Increasing genetic and archaeological evidence continue to support the idea that the corridor was a route for dispersal in both directions, with consequences for biological and technological diversity on both sides of the receding ice sheets.

MATERIALS AND METHODS

Radiocarbon Dating

Collagen was extracted from bone and tooth samples using ultra-filtration methods outlined in Beaumont W et al. (2010) 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 13.1 and Dataset S1). AMS dates were calibrated using the IntCa113 calibration curve (Reimer et al., 2013) 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.

DNA Extraction, PCR, Cloning, Library Preparation, Target Enrichment, and Sequencing

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 13.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 (Rohland et al., 2010; Dabney et al., 2013), 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) (Shapiro et al., 2004). 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 cleanup protocol of Rohland and Reich (2012).

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 Shapiro et al. (2004); 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 Meyer and Kircher (2010) with modifications from Heintzman et al. (2015). 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.

Consensus Sequence Construction

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 (Li and Durbin, 2009), 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 (Green RE et al., 2008). 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.

Phylogeographic Analyses

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 (Drummond et al., 2012). 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 (Gill et al., 2013). To infer the timing and directionality of movement between Beringia and interior North America, we adopted the discrete phylogeographic model described in Edwards et al. (2011), 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 (Gowan, 2013). 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 (Minin and Suchard, 2008). 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 Log-Combiner 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.

Data Availability

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.

ACKNOWLEDGMENTS

We thank James Cahill and Andrew Fields for technical assistance, Mike Waters for providing samples, and two anonymous reviewers for comments that improved the manuscript. P.D.H., A.E.R.S., and B.S. were funded by National Science Foundation (NSF)-Applied Research Center Grants 1203990 and 1417036 and Gordon and Betty Moore Foundation Grant GBMF3804. D.F. was funded by the Canada Research Chairs program and a grant from the Natural Science and Engineering Research Council of Canada. J.W.I. was funded by the Landrex Distinguished Professorship. M.A.S. was funded by NSF-Division of Mathematical Sciences Grant 1264153 and National Institutes of Health Grant R01 HG006139.

This page intentionally left blank.