Rise and Fall of the Beringian Steppe Bison

See allHide authors and affiliations

Science  26 Nov 2004:
Vol. 306, Issue 5701, pp. 1561-1565
DOI: 10.1126/science.1101074


The widespread extinctions of large mammals at the end of the Pleistocene epoch have often been attributed to the depredations of humans; here we present genetic evidence that questions this assumption. We used ancient DNA and Bayesian techniques to reconstruct a detailed genetic history of bison throughout the late Pleistocene and Holocene epochs. Our analyses depict a large diverse population living throughout Beringia until around 37,000 years before the present, when the population's genetic diversity began to decline dramatically. The timing of this decline correlates with environmental changes associated with the onset of the last glacial cycle, whereas archaeological evidence does not support the presence of large populations of humans in Eastern Beringia until more than 15,000 years later.

Climatic and environmental changes during the Pleistocene epoch [from 2 million years ago (Ma) to 10,000 years before the present (ky B.P.)] played an important role in the distribution and diversity of modern plants and animals (1, 2). In Beringia, local climate and geology created an ice-free refugium stretching from eastern Siberia to Canada's Northwest Territories (3). Periodic exposure of the Bering Land Bridge facilitated the exchange of a diverse megafauna (such as bison, mammoth, and musk ox) supported by tundra-steppe grasses and shrubs (3, 4). Humans are believed to have colonized North America via this route, and the first well-accepted evidence of human settlement in Alaska dates to around 12 ky B.P. (5). The latest Pleistocene saw the extinction of most Beringian megafauna including mammoths, short-faced bears, and North American lions. The reasons for these extinctions remain unclear but are attributed most often to human impact (6, 7) and climate change associated with the last glacial cycle (8).

Pleistocene bison fossils are abundant across Beringia and they provide an ideal marker of environmental change. Bison are believed to have first entered eastern Beringia from Asia during the middle Pleistocene [marine oxygen isotope stages (MISs) 8 to 6, circa (ca.) 300 to 130 ky B.P.] and then moved southward into central North America during the MIS 5 interglacial period (130 to 75 ky B.P.), where they were distributed across the continental United States (9). During this time, Beringian and central North American bison populations may have been periodically separated by glacial ice that formed over most of Canada (10, 11). The timing and extent of genetic exchange between these areas remain unclear (2).

The abundance and diversity of bison fossils have prompted considerable paleontological and archaeological research into their use as stratigraphic markers. Extensive morphological diversity, however, has complicated discrimination between even the most accepted forms of fossil bison, and the lack of stratigraphy in Beringian sites has prevented the development of a chronological context. These complications create a complex literature of conflicting hypotheses about bison taxonomy and evolution (9, 12). After a severe population bottleneck, which ocurred only 200 years ago (13), two subspecies survive in North America: Bison bison bison, the plains bison, and B. b. athabascae, the wood bison (9, 13).

To investigate the evolution and demographic history of Pleistocene bison, we collected 442 bison fossils from Alaska, Canada, Siberia, China, and the lower 48 United States (14). We used ancient DNA techniques to sequence a 685–base pair (bp) fragment of the mitochondrial control region (14). Accelerator mass spectrometry radiocarbon dates were obtained for 220 samples, which spanned a period of >60 ky (14).

The association of radiocarbon dates with DNA sequences enables the calibration of evolutionary rates within individual species (15). Bayesian phylogenetic analyses produced an evolutionary rate estimate for the bison mitochondrial control region of 32% per million years (My) [95% highest posterior density (HPD): 23 to 41% per My] (14). This estimate is independent of paleontological calibrations but agrees with fossil-calibrated rates for cattle of 30.1% per My (16) and 38% per My (17). This rate was used to calculate the ages of key nodes in the bison genealogy (14). The most recent common ancestor (MRCA) of all bison included in this analysis lived around 136 ky B.P. (95% HPD: 164 to 111 ky B.P.). In the majority (66%) of estimated trees, Eurasian bison cluster into a single clade, with a MRCA between 141 and 89 ky B.P. Although these two estimates overlap, the age of the MRCA of Eurasian bison was the same as that of the root in 4.8% of 135,000 posterior genealogies (with a Bayes factor of 20.83 that the Eurasian MRCA is not also the MRCA of all clades), suggesting that the Eurasian clade is not the oldest in the tree. This suggests that late Pleistocene bison from the Ural Mountains to northern China are descendants of one or more dispersals from North America. Several North American lineages fall within the Eurasian clade, indicating subsequent asymmetric genetic exchange, predominantly from Asia to North America.

Figure 1A depicts inferred gene flow between bison populations in Beringia and central North America during MIS 3 (∼60 to 25 ky B.P.), which is the interstadial period before the Last Glacial Maximum (LGM, ca. 22 to 18 ky B.P.). Bison were continuously distributed from eastern Beringia southward into central North America during this period, before the formation of the Laurentide (eastern) and Cordilleran (western) ice sheets created a barrier to north-south faunal exchange. Although any coalescence between these ice masses was brief (11), the absence of faunal remains aged 22 to12 ky B.P. (Fig. 1B) (18) indicates that the area was uninhabitable by large mammals during this time. Bison fossils in central North America during the LGM are sparsely distributed across the continent (9). DNA could be retrieved only from two specimens from this period, both from Natural Trap Cave, Wyoming (20,020 ± 150 and 20,380 ± 90 ky B.P.). These specimens are not closely related (14), indicating that populations south of the ice retained some genetic diversity until the LGM.

Fig. 1.

Distribution of bison in Beringia and central North America through time. (A to C) Double-headed arrows show gene flow between regions. Black arrows indicate colonization events. Circles in maps (D) and (E) designate either northern (red) or southern (blue) ancestry and the number of samples from that location.

The ice sheets began to retreat around 14 ky B.P., forming an ice-free corridor (IFC) through which dispersal between Beringia and North America could occur. The first observed bison haplotypes in the IFC are southern in origin (Fig. 1, C and D), with the oldest specimen being in southern Alberta by 11.3 ky B.P., and others near Athabasca, northern Alberta, by 10.4 ky B.P. This finding is consistent with evidence that the first faunal assemblages and archaeological presence in the IFC were southern in origin (1820). The opening of the northern end of the IFC saw a limited southward dispersal of Beringian bison, with a subset of the northern diversity found near the Peace River (northwestern British Columbia) by 11.2 to 10.2 ky B.P. (Fig. 1C) (14). Southern bison are also found in this area around 10.5 ky B.P., making it the only location where post-LGM northern and southern clades occurred at the same time. Subsequent genetic exchange between Beringia and central North America was limited by the rapid establishment of spruce forest across Alberta around 10 ky B.P. (21) and by the widespread development of peatland across western and northwestern Canada (22). North of these ecological barriers, grasslands were reduced by invading trees and shrubs, yet despite the decrease in quality and quantity of habitat (3), bison persisted in eastern Beringia until a few hundred years ago (14, 23).

It has been hypothesized that modern bison descended from Beringian bison that moved south through the IFC after the LGM (9, 19) and have since undergone a decline in diversity due to over-hunting and habitat loss (13). In contrast, our data show that modern bison are descended from populations that were south of the ice before the LGM and that diversity has been restricted to at least 12 ky B.P., around the time of the megafaunal extinctions. All modern bison belong to a clade distinct from Beringian bison. This clade has a MRCA between 22 and 15 ky B.P., which is coincident with the separation of northern and southern populations by the western Canadian ice barrier. This clade diverged from Beringian bison by 83 to 64 ky B.P. and was presumably part of an early dispersal from Beringia, as indicated by the long branch separating it from Beringian bison (14). If other remnants of these early dispersals survived the LGM, they contributed no mitochondrial haplotypes to modern populations.

Coalescent theory is used to evaluate the likelihood of a demographic history, given plausible genealogies (24). Under a coalescent model, the timing of divergence dates provides information about effective population sizes through time. To visualize this for bison, a technique called the skyline plot was used (14, 25). The results showed two distinct demographic trends since the MRCA, suggesting that a simple demographic model, such as constant population size or exponential growth, was insufficient to explain the evolutionary history of Beringian bison. We therefore extended the Bayesian coalescent method (26) to a two-epoch demographic model with exponential population growth at rate rearly, until a transition time, ttrans, after which a new exponential rate, rlate, applies until the present effective population size, N0, is reached (Fig. 2A). In this model, both the early and late epochs can have positive or negative growth rates, with both the rates and the time of transition estimated directly from the data.

Fig. 2.

(A) The two-epoch demographic model with four demographic parameters: N0, rearly, rlate, and ttrans. The effective population size is a compound variable considered linearly proportional to census population size. (B) Log-linear plot describing the results of the full Bayesian analyses. Smoothed curves provide mean and 95% HPD (blue-shaded region) values for effective population size through time. Dashed vertical lines and gray-shaded regions describe mean and 95% HPD for the estimated time of the MRCA (111 to 164 ky B.P.), transition time (32 to 43 ky B.P.), and the earliest unequivocal reported human presence in eastern Beringia (∼12 ky B.P.) (5). The stepped line is the generalized skyline plot derived from the maximum a posteriori tree of the exponential growth analysis. The bar graph shows the number of radiocarbondated samples in bins of 1000 radiocarbon years. No relation is apparent between the absolute number of samples and the estimated effective population size or transition time.

The analysis strongly supported a boom-bust demographic model (Table 1), in which an exponential expansion of the bison population was followed by a rapid decline, with a transition around 37 ky B.P. (Fig. 2B). At the height of the boom, the population size was around 230 times (95% HPD: 71 to 454 times) that of the modern population. When this model is applied to the modern clade alone, a growth period peaks around 1000 years ago (95% HPD: 63 to 2300 yr B.P.) and is followed by a rapid decline (14), which is consistent with historical records of a population bottleneck in the late 1800s (13). These results illustrate the power of this method to recover past demographic signals.

Table 1.

Results of Bayesian analyses assuming constant population size, exponential growth, and a two-epoch model for the full analysis of 191 bison associated with finite radiocarbon dates (14). Model parameters are as defined in (26). The large difference between the mean goodness-of-fit statistics [ln(posterior)] indicates that under either the Akaike information criterion or Bayesian information criterion tests, the two-epoch model is a significantly better fit to the data than the simpler models.

Constant size Exponential growth Two epoch
Lower Mean Upper Lower Mean Upper Lower Mean Upper
Age estimates (yr B.P.)
    Root height 117,000 152,000 189,000 113,000 146,000 181,000 111,000 136,000 164,000
    Modern/southern clade 20,200 28,000 36,600 18,600 26,400 35,000 15,400 23,200 32,200
    Eurasian clade 85,000 116,000 151,000 83,000 112,000 144,000 89,000 114,000 141,000
Model parameters
    Mean ln(posterior) -6530.795 -6517.35 -6394.568
    Mutation rate (substitutions/site/year) 2.79 × 10-7 3.78 × 10-7 4.85 × 10-7 2.30 × 10-7 3.20 × 10-7 4.13 × 10-7 2.30 × 10-7 3.20 × 10-7 4.13 × 10-7
    Kappa 19 27 37 19 27.4 37 19 27 37
    Shape parameter 0.22 0.35 0.49 0.22 0.35 0.49 0.22 0.35 0.5
    Proportion of invariant sites 0.33 0.45 0.56 0.33 0.45 0.56 0.34 0.45 0.56

The effects of population subdivision and patch extinction and recolonization on coalescence patterns have not been fully characterized, yet they can influence demographic estimates such as skyline plots (27). To test for the effect of population subdivision on our models, the two-epoch analysis was repeated first without the Eurasian bison and then without both Eurasian and central North American bison. The results of these analyses were consistent with those for the entire data set (14), suggesting that the assumption of panmixia does not affect the analysis. These results suggest that the major signal for the boom-bust scenario came from the well-represented eastern Beringian population.

The timing of the decline in Beringian bison populations (Fig. 2B) predates the climatic events of the LGM and events at the Pleistocene-Holocene boundary. The bison population was growing rapidly throughout MIS 4 and 3 (∼75 to 25 ky B.P.), approximately doubling every 10,200 (95% HPD: 7500 to 15,500) years. The reversal of this doubling trend at 42 to 32 ky B.P. and the subsequent dramatic decrease in population size are coincident with the warmest part of MIS 3, which is marked by a reduction in steppe-tundra due to treecover reaching its late Pleistocene maximum (28). Modern boreal forests serve as a barrier to bison dispersal because they are difficult to traverse and provide few food sources (3). After the interstadial, cold and arid conditions increasingly dominated, and some component of these ecological changes may have been sufficient to stress bison populations across Beringia. Previous reports of local extinction of brown bears (29) and hemionid horses (8) in Alaska around 32 to 35 ky B.P. support the possibility of a larger scale environmental change affecting populations of large mammals.

These results have considerable implications for understanding the end-Pleistocene mass extinctions, because they offer the first evidence of the initial decline of a population, rather than simply the resulting extinction event. These events predate archaeological evidence of significant human presence in eastern Beringia (5), arguing that environmental changes leading up to the LGM were the major cause of the observed changes in genetic diversity. If other species were similarly affected, differences in how these species responded to environmental stress may help to explain the staggered nature of the megafaunal extinctions (7, 30). However, it is possible that human populations were present in eastern Beringia by 30 ky B.P., with reports of human-modified artifacts as old as 42 to 25 ky B.P. from the Old Crow basin in Canada's Yukon Territory (31). Although the archaeological significance of these specimens is disputed and the number of individuals would be low, the specimens are consistent with the timing of the population crash in bison. This emphasizes that future studies of the end–Pleistocene mass extinctions in North America should include events before the LGM.

Ancient DNA is a powerful tool for studying evolutionary processes such as the response of organisms to environmental change. It should be possible to construct a detailed paleoecological history for late Pleistocene Beringia using similar methods for other taxa. Almost none of the genetic diversity present in Pleistocene bison survived into Holocene populations, erasing signals of the complex population dynamics that took place as recently as 10,000 years ago.

Supporting Online Material

Materials and Methods

SOM Text

Figs. S1 to S5

Tables S1 to S4


References and Notes

View Abstract

Navigate This Article