Research ArticlePALEOECOLOGY

Climate, CO2, and the history of North American grasses since the Last Glacial Maximum

See allHide authors and affiliations

Science Advances  25 Mar 2016:
Vol. 2, no. 3, e1501346
DOI: 10.1126/sciadv.1501346


The spread of C4 grasses in the late Neogene is one of the most important ecological transitions of the Cenozoic, but the primary driver of this global expansion is widely debated. We use the stable carbon isotopic composition (δ13C) of bison and mammoth tissues as a proxy for the relative abundance of C3 and C4 vegetation in their grazing habitat to determine climatic and atmospheric CO2 controls on C4 grass distributions from the Last Glacial Maximum (LGM) to the present. We predict the spatial variability of grass δ13C in North America using a mean of three different methods of classification and regression tree (CART) machine learning techniques and nine climatic variables. We show that growing season precipitation and temperature are the strongest predictors of all single climate variables. We apply this CART analysis to high-resolution gridded climate data and Coupled Model Intercomparison Project (CMIP5) mean paleoclimate model outputs to produce predictive isotope landscape models (“isoscapes”) for the current, mid-Holocene, and LGM average δ13C of grass-dominated areas across North America. From the LGM to the present, C4 grass abundances substantially increased in the Great Plains despite concurrent increases in atmospheric CO2. These results suggest that changes in growing season precipitation rather than atmospheric CO2 were critically important in the Neogene expansion of C4 grasses.

  • stable isotopes
  • last glacial maximum
  • C4 grasses
  • paleoecology
  • bison
  • mammuthus


C4 photosynthesis currently contributes to ~20 to 25% of global terrestrial productivity (1), but the conditions under which this photosynthetic pathway evolved and spread have remained unclear (24). C4 photosynthesis is an adaptation that increases the efficiency of carbon fixation at low atmospheric CO2 concentrations and under high-light and high-temperature conditions, whereas C3 photosynthesis remains more efficient under conditions of higher atmospheric CO2 and under lower temperatures (2, 5). Given these climatic and atmospheric controls on photosynthesis, changes in atmospheric CO2 and/or temperature were originally suggested as the cause of the increase in global abundance of C4 grasses in the late Neogene (6). However, reconstructions of atmospheric CO2 concentrations (7) as well as paleotemperatures (8) show stable conditions during the Miocene-Pliocene transition, implying a different driver such as hydrologic change (4, 9).

Previous studies have tried to determine the climatic control on the distribution of C3 and C4 vegetation in North America using statistical methods, modern vegetation abundances, and δ13C of soil organic matter (SOM) (1012). However, approaches that focus solely on modern climate-vegetation relationships cannot quantify the effect of changing atmospheric CO2 concentrations through time because of differences in CO2 responses between C3 and C4 vegetation. Additionally, while process-based models have been shown to accurately predict the distribution of biomes from the modern through the Pleistocene (13, 14), many models, including the latest Coupled Model Intercomparison Project (CMIP5) Earth system models, have yet to reach a consensus in predicting the modern and past distribution of C3 and C4 vegetation at a range of scales (1, 15, 16), and both statistical and process-based model results are highly variable. Here, we use the δ13C of modern and fossil bison (Bison bison), and fossil mammoth (Mammuthus columbi and Mammuthus primigenius) hair, tooth enamel, and bone collagen from the Last Glacial Maximum (LGM) to the present as a proxy for the abundance of C3 and C4 vegetation in the animal’s habitat. Differences in photosynthetic isotope fractionation between C3 and C4 plants (17, 18) result in average values of ~−27‰ and −12‰ (19) for C3 and C4 plant leaf isotope composition, respectively. Bison and mammoths are obligate grazers, and bison diets correlate closely with local abundance of C3 and C4 grasses (20). Mammoth diets appear to display a similar correlation with grass abundances (fig. S1). Additionally, low δ13C variability among individuals from the same site suggests that many prehistoric bison in the Great Plains had limited foraging ranges (21, 22). The average δ13C of grasses is represented in their tissues with known discrimination values (23), and thus, we can use the δ13C of the grazer’s diet to estimate the C4 abundance in their habitat to within 10% (20, 24). Grazers integrate the isotopic composition of vegetation over tens of square kilometers, a much larger area than can be incorporated by spot sampling of soils or vegetation abundances, and their diets should average ecosystem heterogeneity. The LGM to the modern is an ideal time period to study the influence of both climate and CO2 on vegetation distributions because of the abundance of modern and fossil data available compared to the Miocene and Pliocene.

We combine these isotope data with high-resolution modern climate data (25), CMIP5 paleoclimate model outputs (26), and atmospheric CO2 concentrations through time to constrain C3/C4 climatic controls and produce a multivariate predictive isotope landscape model (“isoscape”) of the distribution of C3 and C4 grasses across western North America for the modern, mid-Holocene, and LGM. We improve on previous statistical methods based solely on modern data by using both modern and fossil isotope data to allow us to incorporate changes in atmospheric CO2 concentrations into our model. We use three different classification and regression tree (CART) methods of recursive partitioning [random forest (RF) (27), conditional forest (CF) (28), and boosted regression tree (BRT) (29)] to infer the climate and CO2 controls on the distribution of C3 and C4 grasses (see Materials and Methods for further discussion). Decision trees and CART analyses have advantages over conventional, parametric statistical regressions because they do not require normally distributed data, have superior predictive capability, easily incorporate both categorical and continuous data in the same analysis, capture interactions between environmental predictor variables, and can distinguish threshold values and responses of predictor variables (27, 30).


The isotope data set consists of 632 bison and mammoth individuals from 115 locations across North America across a time period that spans the LGM through the modern (Fig. 1). Bison specimens were dated using radiocarbon analysis or association with archeological remains. See individual references (table S1) for dating methods of particular fossil specimens. Accounting for trophic enrichment in 13C in different tissues and changes in the δ13C of atmospheric CO2 (see Materials and Methods), the δ13C values of bison and mammoth dietary vegetation range from −11.3 to −29.8‰, encompassing the entire range of δ13C values from both C3 and C4 photosynthesis. Each sample was assigned the following climate variables for its respective time period: growing season temperature (GST) and growing season precipitation (GSP), defined as the mean temperature and sum of precipitation in months that have mean temperatures exceeding 5°C while also receiving a minimum of 25 mm of precipitation, mean annual temperature (MAT), mean annual precipitation (MAP), and the proportion of precipitation falling in both the winter (WP) and the summer (SP) (10). Each sample was also assigned a specific number of months exceeding a crossover temperature [COT; a physiological modeled temperature, dependent on atmospheric CO2 concentration (2), in which C4 grasses are predicted to outcompete C3 grasses because of higher C4 photosynthetic rates], while receiving a minimum of 25-mm precipitation in that month (table S1) (1, 15). GSP, GST, and COT were calculated from high-resolution interpolated observations for modern specimens and the mean of six CMIP5 climate model outputs for the fossil specimens (see Materials and Methods). Given that our study uses mean conditions during each climatic interval (the mean of ~100 CMIP5 model years as well as mammal δ13C, which incorporate months to years of time), we cannot assess the influence of drought events on the distribution of grasses. Instead, we include aridity [PET (potential evapotranspiration)/MAP] as an indicator of water deficit as a climate variable in the CART analysis. We have not included fire return interval in these analyses because of the lack of high-resolution continental-scale fire data from the Holocene and Pleistocene. The modern vegetation isoscape was produced using 1901–1950 average climate data to simulate natural vegetation abundance under climatic conditions similar to those occurring before the spread of agriculture (that is, before 1800 for North America). The mid-Holocene and LGM isoscapes were produced using CMIP5 mid-Holocene and LGM mean climate model outputs to determine the effect of changing atmospheric CO2 on vegetation distributions.

Fig. 1 Location of bison and mammoth samples used in the CART analyses.

White circles represent modern localities and red diamonds represent fossil localities. The references for these sites are listed in table S1. There are 632 individual specimens in total; 281 of these samples are identified as modern (dated to within the last 300 years), 256 specimens are identified as Holocene, and 95 specimens are dated to the LGM.

Variable importance

The predictive CART regressions explain up to 92.1% of the variability in the δ13C of bison and mammoth tissues, with an average error of ±1.55‰ on the validation (out-of-bag) data set. For the BRT and RF CART methods, the number of months above the COT at which C4 outcompetes C3 vegetation (which varies with atmospheric CO2) was the most important predictor of the distribution of C4 grasses across North America. Because this COT variable includes temperature, precipitation (with a threshold of 25 mm per month for C4 growth), and atmospheric CO2, it is necessary to exclude it from the variable importance calculations to determine the importance of single climate variables. When COT is excluded from the analysis, GSP is the most important variable in the BRT and RF methods, followed by GST and MAT (Fig. 2). In both the BRT and RF methods, atmospheric CO2 ranks fifth in importance. The CF results do not match those of the BRT and RF, which ranks MAT as the most important variable, followed by GST, COT, atmospheric CO2, and GSP. When the COT is removed from the variable importance analysis, CF still ranks MAT first, followed by GST and atmospheric CO2 (Fig. 2). Although each method’s variable importance results do not exactly match, BRT and RF results are very similar and have exactly the same ranking for the top three most important variables. This similarity gives us confidence that these results are the most accurate and reliable.

Fig. 2 Variable importance.

Variable importance in determining bison and mammoth δ13C as inferred by the three different CART methods is shown. Because COT is a compound variable that incorporates temperature, precipitation, and atmospheric CO2, two analyses were executed to assess the overall importance of all variables and single climate variables. We only display the analysis excluding COT to show the importance of single climate variables. WP, proportion of precipitation falling in the winter; SP, proportion of precipitation falling in the summer; MSE, mean standard error.

The BRT and RF methods return the partial dependence of each variable on the overall isotopic variability of bison and mammoth tissues. The partial dependence is calculated as the variance explained by a variable when all others are held constant. The compound COT variable explains ~3.5‰ of the total isotopic variability (~25% of variability in C4 abundance; Fig. 3A). When the COT variable is removed from the analysis, GSP explains ~3.5‰ of the total isotopic variability (Fig. 3B), GST explains ~2.5‰ of the total isotopic variability (Fig. 3C), and atmospheric CO2 explains ~0.7‰ of the total isotopic variability (equivalent to ~5% variability in C4 abundance; Fig. 3D). COT, GSP, and GST are positively correlated with δ13C and C4 abundance (higher C4 abundances with higher temperature and precipitation), whereas atmospheric CO2 is negatively correlated with δ13C (higher abundances of C4 at lower levels of atmospheric CO2). Although we are confident that the model reliably captures the effect of CO2 at preindustrial to modern atmospheric concentrations, the CART regressions may be underestimating the importance at glacial concentrations of CO2 because of the smaller number of LGM specimens in the overall data set (~15%).

Fig. 3 Partial dependence of individual variables.

(A to D) Partial dependence of the COT variable (A), GSP (B), GST (C), and atmospheric CO2 (D) on the δ13C of bison and mammoth diets.

Modern isoscape

Figure 4 presents the mean isoscapes of the three CART prediction methods for the three time periods (modern, mid-Holocene, and LGM). Regression results are similar for each CART method, although the CF method infers different importance rankings. The BRT method predicts higher abundances of C4 grasses than RF and CF in all time periods, but the pattern of grass distributions and general percentages of C4 abundance agree between each method. These results as well as error maps can be found in figs. S2 to S6. Because this predictive map is only applicable to regions with a climate conducive to the growth of grasses, the forested regions are not displayed (see Materials and Methods for the discussion of vegetation masks). For the modern (Fig. 4, A and B), the regions with high δ13C values (that is, dominated by C4 grasses) occur in the southern Great Plains, with a swath of the highest δ13C values extending from northern Mexico into the tallgrass prairie peninsula region of Iowa and Missouri. The modeled north-south transition between C4-dominated and C3-dominated regions occurs across South Dakota (~44°N). West of the Great Plains in the Rocky Mountains, C3 vegetation is typically dominant due to higher elevations and lower GSTs and precipitation. In the southwest United States, the CART regressions predict C3-dominated vegetation at high elevations but 40 to 70% C4 in the Mojave and Sonoran deserts, likely due to the summer precipitation from the North American monsoon. The predicted C4 abundances are in broad agreement with measured C4 relative abundance from sites across North America (circles in Fig. 4B) (10).

Fig. 4 Carbon isoscapes.

Top: Predicted modern average δ13C of vegetation consumed by bison across western North America. The filled area represents the predicted δ13C values based on the mean of the three CART regression methods, and the dots are the differences between the measured and predicted δ13C of diet values at sample locations. (A, C, and E) Modern (A), mid-Holocene (C), and LGM conditions (E), respectively. Middle: Corresponding predicted modern relative abundance of C4 grasses converted from δ13C values based on typical C3 and C4 isotopic composition end members. (B, D, and F) Modern (B), mid-Holocene (D), and LGM conditions (F), respectively. Overlaying the modern C4 abundance prediction map (B) is the measured relative abundance of C4 vegetation for a number of sites across North America from a separate study (10). Gray squares in (E) and (F) represent areas covered by ice during the LGM. Bottom: Changes in the concentration of atmospheric CO2 from the LGM to the present (64). yBP, years before present.


For the mid-Holocene (Fig. 4, C and D), the models predict generally similar patterns in the δ13C and distribution of C4 grasses to the modern isoscape, except for higher abundances of C4 grasses in the northern Great Plains and slightly lower abundances in the central and southern Great Plains in the mid-Holocene compared to the modern. From the mid-Holocene to the modern, the CART models predict moderate decreases in the abundance of C4 grasses in the northern Great Plains as well as slight decreases in northern Mexico, the western half of Texas, and New Mexico (Fig. 5A). The mid-Holocene had a generally stable climate that was slightly warmer than modern (31), likely caused by changes to orbital forcing (32). This warmer climate as well as slightly lower atmospheric CO2 increased the abundances of C4 grasses in North America.

Fig. 5 Change in the fractional abundance of C4 grasses.

(A) From the mid-Holocene to the present. (B) From the LGM to the present. The colored areas represent the intersection of modern and past grass-dominated areas.

Last Glacial Maximum

The CART regressions explain 91% of the variability in the LGM isotope data and predict low δ13C and low abundance of C4 grasses in the central and northern Great Plains, driven mainly by low temperature during the LGM (Fig. 4, E and F). However, the isoscapes predict high abundances of C4 grasses similar to the modern along the Texas coast, southern Arizona, and northwest Mexico. The rest of the continent is predicted to have been dominated by C3 vegetation during the LGM. C4 grasses are dominant during LGM only in regions with an increase in mean summer temperature less than ~6° to 8°C from the LGM to the modern (fig. S7). This result is consistent with physiological models that estimate a decrease in the crossover GST of C3/C4 dominance of ~8°C between 300 and 200 ppm (2, 15), which would predict an increase in C4 grasses during the LGM in areas where the change in temperature was less than ~8°C. Figure 5B shows the modeled fractional change in abundance of C4 vegetation from the LGM to the present. Similar to previous work (24), the models predict little change in C4 abundances along the Gulf Coast and an absolute increase of 30 to 50% from the LGM to the present in C4 abundance in the central Great Plains (Fig. 5B). This increase suggests that temperature and precipitation changes primarily controlled the distribution of C3 and C4 grasses at lower latitudes during the LGM except in the southern Great Plains, where smaller temperature changes allowed for CO2 changes to modulate vegetation composition.


To validate the paleo-isoscapes, we compare the mid-Holocene and LGM predictions to independent paleovegetation proxies from soils. The δ13C of bulk SOM and pedogenic carbonates record the mean δ13C of all vegetation growing at the time of soil formation (17). By contrast, the δ13C of bison tissue only records the mean δ13C of herbaceous vegetation, excluding woody vegetation. Therefore, our isoscape and the soil δ13C proxies are only exactly comparable when vegetation is 100% herbaceous. The soil δ13C should otherwise be more negative than bison δ13C values in areas of mixed woody and herbaceous vegetation. The soil-based proxy reconstructions can only be compared to our isoscape generally and cannot be used as an approximation of CART model accuracy because of these differences in the aspect of the environment that each proxy captures. A number of soil δ13C measurements from the mid-Holocene are in agreement with our predictions (3335). In particular, Nordt et al. (35) find that the δ13CSOM of many mid-Holocene soils in the central and southern Great Plains were more negative than modern δ13CSOM, suggesting lower abundances of C4 grasses. Our isoscape predicts a 10 to 30% increase in C4 abundance in the central and southern Great Plains from the mid-Holocene to the present (Fig. 5A), consistent with the findings from δ13CSOM. For the LGM, there are few records of δ13CSOM as well as δ13C of packrat middens in western North America, but those that do exist are in agreement with our isoscape. In central Kansas, both SOM (36) and our isoscape predict vegetation δ13C of ~−22‰, substantially more negative than ~−15‰ of present-day values. In eastern Colorado (37) and southeastern Arizona (38), SOM and packrat middens [which record the diets of Neotoma sp. (39)] reconstruct similar vegetation δ13C values to our isoscape, values that are slightly more negative than present-day values. In southern Texas, δ13CSOM records are more negative than our isoscape predictions (40). This divergence is likely due to the presence of moderate abundances of oak, pine, and alder in the area during the LGM (41, 42).

There are a number of independent sources of error, and ways to assess these errors, in our CART-based prediction maps for the distribution of C4 grasses in North America. One possible source of error is that these isoscapes combine bison δ13C values from three different tissue types, which each record diet at different times in the bison’s or mammoth’s life. The turnover rates for bone collagen are slow, and so, the carbon isotopic composition of bone collagen reflects the dietary inputs of C3 and C4 vegetation during bone growth. If the bone growth is seasonal, then that collagen may be reflective of seasonal vegetation composition rather than annual average composition that would be recorded in the tooth enamel. Additionally, hair records the previous few weeks of the diet of the animal (43). If the hair grows continuously, the length of time of the diet that is represented by the hair is proportional to the length of the hair sampled. Most bison tooth enamel data for this analysis come from the data set of Hoppe et al. (20), which sampled roughly 1 year of growth from the third molar of the bison, which forms after weaning (44, 45). Table S2 displays the average isotope composition of bison diet derived from collagen, hair, and tooth enamel for eight sites with analyses of multiple types of bison tissue. There is no systematic enrichment or depletion in any particular tissue type, and in all but one site, the range of δ13C values between tissues is less than 1‰, indicating that differing time scales of tissue formation are unlikely to substantially influence the model.

Error is assessed in the fit of the CART-based regression to the data. The RF regression cross-validates the results by randomly selecting a subset of the data and training a tree on that subset. The fit of that tree is then tested on the excluded data. This is repeated for the number of trees that are included in the regression (1000 to 5000, depending on the method). The error associated with the prediction map from the CART regressions is ±1.55‰, and the regression explains 92% of the variability in the isotopic data. Because the modern isoscape is generated using one observationally based high-resolution global climate data set, error is attributed only to variable uncertainty and the regression. For the paleo-isoscapes, there is also variability associated with the CMIP5 paleoclimate model outputs used to create the prediction maps. We used the mean of six paleoclimate models, and to assess how the variation in those models would affect our isoscapes, we created isoscapes using ±1 standard deviation (SD) (σ) around the mean of the six climate models for the LGM and the mid-Holocene. Figures S5 and S6 show the differences in the isoscape predictions when incorporating variability from the climate models. For the LGM, this variation is fairly low, and the models still predict high abundances of C4 grasses in the southern Great Plains region regardless of the variation in climate data. The variability in the mid-Holocene is higher, likely due to the high sensitivity of vegetation to temperature changes in the northern Great Plains. Figure S8 shows the SD for months above COT, GSP, and GST for the six models simulating the LGM. The SD and model variability for these most important climate variables are low for much of the area across North America for which we have predicted δ13C values. This low SD contributes to the low predicted variability in the LGM isoscapes. We also calculated variable importance rankings using ±1σ around the mean of the six climate models for the LGM and the mid-Holocene to assess the impact of model uncertainty on our climate variable importance results (fig. S9). Varying the climate data does not change the importance rankings substantially, suggesting that the result of greater importance of precipitation and temperature over atmospheric CO2 is robust.

This work represents the first prediction for changes to the distribution of C3 and C4 grasses from the LGM to the present based on a combination of isotope data, advanced statistical methods, and paleoclimate model outputs. Our CART analyses and predictive isoscapes show that, where climate changes are small to moderate, low atmospheric CO2 concentrations such as those observed during the LGM will favor increases in the abundance of C4 grasses. This agrees with physiological simulations, which predict the greatest sensitivity to atmospheric CO2 concentrations below ~300 ppm (15). These results suggest that low atmospheric CO2 concentrations only favor C4 grasses at lower latitudes where temperatures remain warm. However, the large increases in C4 grass abundance in the central and northern Great Plains from the LGM to the present as atmospheric CO2 rose from ~190 to ~280 ppm in the Holocene imply that climate, instead of atmospheric CO2, has a much stronger influence on the balance between C3 and C4 vegetation at higher latitudes.

These results have important implications for the Neogene expansion of C4 grasses. Whereas atmospheric CO2 decreased during the time at which the C4 photosynthetic pathway evolved, there is no evidence for a drop in atmospheric CO2 during the global spread of C4 grasses from the late Miocene through the early Pliocene (7, 8). Our CART models suggest that atmospheric CO2 primarily exerts a control on the distribution of C4 grasses when atmospheric CO2 concentrations are below preindustrial values, levels that were not reached until long after the initial spread of C4 grasses in the glacial periods of the Pleistocene (8). Therefore, these results are consistent with the broader conclusions that atmospheric CO2 alone could not have caused the spread of C4 grasses in the late Neogene (9). Global average temperature was also stable through this ecological transition (8), and so it is unlikely that a global warming event drove C4 expansion. Because GSP is the most important single variable in the model, this work instead suggests that hydrologic changes, more specifically increases in warm-season precipitation, may have driven the late Neogene global expansion of C4 grasses, and independent proxy evidence in the Himalayas (4, 9) and in east Africa (46) supports this scenario. Because both atmospheric CO2 and climate continue to change in the future as a result of anthropogenic warming, these results also suggest that temperature and precipitation changes, rather than CO2 fertilization, will drive future changes to the distribution of C3 and C4 grasses.


We compiled an extensive data set of 632 bison and mammoth individuals from 115 locations across North America, including new measurement and previously published data from the literature, to train the CART models in predicting the δ13C of grass vegetation. Mammoths (M. columbi and M. primigenius) were also relatively indiscriminant grazers and are used to reconstruct past vegetation changes in North America (24, 47). We include measurements of mammoth bone collagen and tooth enamel δ13C to increase the representation of fossil specimens in the overall data set. To determine whether mammoth enamel and collagen can be used interchangeably with bison tissue in these analyses, we plotted the mammoth collagen and enamel δ13C against bison enamel and collagen δ13C for animals of the same age from the same locality. The data set contains nine sites from North America with both fossil bison and mammoth of the same age, and we also added two more Russian sites for comparison (48, 49). Bison and mammoth tissue δ13C are extremely well correlated and plot close to the 1:1 line (R2 = 0.99; fig. S1), demonstrating that co-occurring mammoth and bison fossils record similar diets and that mammoth can be used to augment the paleo–data set.

Bison and mammoth samples were assigned Global Positioning System (GPS) coordinates based on locality information. For modern bison specimens, the average climate of the specimen’s location was assigned using the high–spatial resolution (30–arc sec) Global Climate Data ( (25). Because of the changes in atmospheric CO2 and climate in the last century (50), the bison data set was split into specimens that lived before and after 1950. Specimens dated to 1700–1950 CE were assigned climatic variables averaged from interpolated observations between 1901 and 1950. The vast majority of specimens from after 1950 lived between 1990 and 2010, and so for this subset of the data, bison were assigned climatic variables from the averaged 1990–2010 observations. Although we do not have gridded climate observations before 1901, preindustrial specimens are assigned an early 20th century average climate before significant anthropogenic warming began (50).

The isotope data set is composed of samples from hair, horn sheath, tooth enamel, and bone collagen from modern specimens through the LGM. New bison hair, horn sheath, enamel, and bone collagen samples were analyzed for this study; additional enamel, horn sheath, and bone collagen δ13C values and location data were compiled from the literature. Hair from bison skins and horn sheaths with accompanying location data were gathered from museum collections across the United States (table S1). Whenever possible, tail or mane hair was sampled because of its continual growth throughout the life of the bison (43). Bison hair samples were washed in a 2:1 chloroform/methanol solution and chopped into ~2-mm sections, and then ~1 to 2 cm of hair were loaded into tin capsules for isotopic analysis (51). The samples were analyzed on Costech Elemental Analyzer attached to a Finnigan MAT 252 isotope ratio mass spectrometer at the University of Utah. Powdered enamel samples representing 1 year of growth were collected from bison third molars (20, 44). Samples were pretreated according to the methods described by Koch et al. (52). Bone collagen samples represent small chips (<150 mg) of compact bone that were removed with a cutoff wheel attached to a variable-speed Dremel tool set at low speed. Bone chips were demineralized by soaking in 0.5 N HCl at 4°C for 24 to 48 hours, then treated with chloroform/methanol and lyophilized. Collagen and enamel samples were analyzed on a Finnigan DELTAplusXP isotope ratio mass spectrometer equipped with Costech Elemental Analyzer and a GasBench II sample preparation device, which was located at Stanford University. All carbon isotope compositions are reported in units of per mil (‰) compared to the international standard Vienna Pee Dee Belemnite (VPDB). Differing mammal tissues undergo varying isotopic enrichment during formation, with an enrichment of 3‰ for hair (23, 53), 5‰ for bone collagen (54), and 14.6‰ for tooth enamel (55). To interpret δ13C values from different materials, each sample was corrected with the respective isotopic enrichment to the δ13C of dietary vegetation carbon. Grazers from this study span a large time range, from the LGM to 2012 CE, and the δ13C of the atmosphere and, subsequently, average vegetation has shifted by approximately −2‰ as a result of the combustion of fossil fuels over that time period (5658). The δ13C values of mammal tissues were corrected to preindustrial atmospheric δ13C values (−6.3‰) (57) based on collection year or date for fossil specimens using the high-resolution atmospheric δ13C record from the Law Dome (56) to account for changes to the atmosphere.

Modern bison were also assigned an atmospheric CO2 level based on the collection year using the annual National Oceanic and Atmospheric Administration (NOAA) atmospheric CO2 measurements at Mauna Loa (58) ( and the Law Dome glacial record (57). For fossil bison and mammoth, GPS coordinates were used to assign climate variables using CMIP5 paleoclimate data (51). For this study, we defined the LGM as between 28,000 and 16,000 calendar years before present (BP) (59) and the mid-Holocene as 8000 to 2000 calendar years BP (31). Paleo-atmospheric CO2 levels were assigned using the European Project for Ice Coring in Antarctica (EPICA) Dome C composite record (56) ( CMIP5 climate simulations for the LGM, the mid-Holocene, and the past 1000 years were obtained from the Earth System Grid Federation ( The six models that had been run for all three paleoclimate time periods used in this study (CCSM4, GISS-E2-R, IPSL-CM5A-LR, MIROC-ESM, MPI-ESM-P, and MRI-CGCM3) were averaged to create one mean grid for each climatic variable for the LGM, the mid-Holocene, and the past 1000 years. Multiple CMIP5 models were used to incorporate a range of simulations and represent uncertainty in the simulations (16). Holocene (<11,700 years BP) bison specimens not dated to the climatically stable time period between 8000 and 2000 years BP (47) were assigned climatic variables from averaged model outputs for the past 1000 years. Modern aridity index (PET/MAP) values were assigned to modern specimens using the CGIAR-CSI (Consultative Group for International Agricultural Research–Consortium for Spatial Information) aridity index data set (60). Aridity index values for the LGM and the mid-Holocene were calculated using CMIP5 data, following the same methods as those used to calculate the modern aridity index outlined by Trabucco and Zomer (60). Assignment of climatic variables and atmospheric CO2 concentrations by time period is summarized in Table 1.

Table 1 Summary of climate and CO2 data.

Summary of the climate and CO2 data assigned to bison and mammoth specimens determined by specimen age.

View this table:

Statistical analyses

Statistical analyses were performed using R software version 3.0.1. The relative importance of each climatic variable on bison δ13C and resulting isoscapes was produced by recursive partitioning (a machine learning technique) using the three different CART analyses: RF [randomForest package (27)], CF [party package (28)], and BRT [gbm package (29)]. In contrast to other statistical methods, machine learning assumes that the ecological processes driving the data distribution are complex and unknown, and works to learn the response by finding patterns and increasing predictive ability. Recursive partitioning is a method of nonlinear multivariable analysis that is useful for determining predictive models when independent variables interact in nonlinear ways, which is represented visually by the decision tree. Single decision trees typically have low predictive value, but incorporation of multiple trees increases the predictive power of the model. Regression trees are advantageous over linear models because they are nonparametric and do not require previous transformation of the data, can incorporate continuous and categorical data in a single model, do not require elimination of outliers, and incorporate interactions between variables (29). The final isoscape is a mean of the predictions from each of the following three CART methods. We have used the three standard approaches to RF-based predictive models to assess how varying the CART methods influences the predictive results.


The randomForest (RF) method builds a model by bootstrapping the data set, uses a subset of data to train each tree, and validates the prediction using the unused portion of the data (out of bag examples). This process is carried out for a set number of decision trees (in this case, 1000 trees) and averages them to produce a stable and single best regression result. Because a portion of the data is used to validate each of the 1000 trees, each data point is used in both the training and the validation data sets, and a priori splitting of the data or additional cross-validation is not necessary. Variable importance is calculated as the increase in prediction error in the out-of-bag samples when a given variable is randomly permuted.

Conditional forest

The CF tree regression method is a CART machine learning technique similar to RF. However, it has been shown that RF can overpredict the importance of variables that are highly correlated (28). Whereas the RF method calculates only the marginal importance, or importance of a single variable without reference to other variables, CF calculates the conditional importance, taking into consideration importance contingent on other variables. Therefore, the CF method should more accurately determine variable importance ranking when variables are highly correlated (28).

Boosted regression tree

The BRT method combines two statistical algorithms: the regression tree as described above and gradient boosting. Gradient boosting is a technique that produces a number of simple models (here, the decision tree) and builds a more complex model in a stepwise fashion. In the BRT method, the model is optimized by the progressive addition of regression trees that best reduce the error on the prediction. The major difference between BRT and the CF and RF methods is that, instead of averaging each tree to produce a single predictive model, BRT combines the trees, adaptively producing a collection of predictive models. The BRT analysis follows the procedure outlined by Elith et al. (29). The relative importance values are calculated as the number of times the decision trees are split on the basis of a particular variable, weighted by the squared improvement to the predictive model at each split. These values are then normalized so that the sum of all importance values equals 100. This method of calculating variable importance also eliminates issues with covariance that arise in the RF method because importance is calculated in consideration of all other variables.

Variable importance

We have calculated variable importance using all three methods above (Fig. 2). Strobl et al. (28) explain how the RF calculation method of variable importance can artificially attribute higher importance to correlated variables. Table S3 displays the Pearson correlation coefficients for each of the variables used in the CART analysis, showing that some of the variables are correlated. We report each of the three variable importance rankings below. Because the method of calculation is different for each, the exact importance scores are not comparable between methods, and only the rankings should be compared among the three.

Predictive isoscapes

In addition to inferring the importance of the climate and CO2 predictors described above, we assessed the importance of soil order in the distribution of modern C4 grasses. Each modern bison sample was assigned a dominant soil order (61) to assess the influence of soil order on vegetation [using the Natural Resources Conservation Service U.S. General Soil Map (STATSGO2), as well as the Soil Landscapes of Canada database from the Canadian Soil Information Service (]. However, soil order was found to be the least important variable for the modern data set, and because of the lack of continental-scale paleosol maps for the mid-Holocene and the LGM, we ultimately dropped this variable from the final analyses. We also assessed the importance of fire for the modern distribution of C4 grasses by assigning each bison specimen a fire return interval using the high-resolution mean fire return interval data set available at Fire was not among the top predictors for each of the CART methods, and the inclusion of this variable did not change the overall predictive capability or error estimates for the regressions. No high-resolution estimate of fire frequency currently exists for the Holocene and Pleistocene, and the fire frequency model outputs from the CMIP5 Earth system models are highly variable, so we also dropped this variable from the final analysis.

To model how the distribution of C3 and C4 grasses has changed in relation to changes in climate and atmospheric CO2 through the Pleistocene, the entire bison and mammoth data set (that is, fossil plus modern samples) was used to incorporate a wide range of climate and CO2 conditions in the CART analyses. The CART regressions were trained on the entire data set spanning the LGM to the modern with the corresponding climate and CO2 values, and then the resulting regressions were applied to the high-resolution 30–arc sec climate grids for the average climatology from 1901 to 1950 (25) across North America and a CO2 level of 300 ppm using the “raster” packages within R to create the modern (pre-1950) isoscape. The resulting isoscape maps (Fig. 4) were produced with ArcGIS version 10.2. The predictive isoscapes are only applicable to regions dominated by herbaceous cover, and the forested and desert regions of the continent have been removed from the map. Modern regions dominated by herbaceous cover were identified using the ecoregion map of the U.S. Environmental Protection Agency (62). The same method as described above for the modern isoscape was used to produce paleo-isoscapes for the LGM and the mid-Holocene, applying the CART regressions to the mean CMIP5 paleoclimate model outputs for the LGM and mid-Holocene scenarios using CO2 levels of 180 and 270 ppm, respectively. Carbon isotopic compositions were converted to abundances of C4 grasses using a two-component mixing model with a pure C3 end member of −25.3‰ and a pure C4 end member of −11.3‰, representing the average isotopic composition of plants under a modern preindustrial atmosphere of −6.3‰. For these paleo-isoscapes, forested regions were removed by creating a mask using the CMIP5 Earth system model outputs for LGM and mid-Holocene percent grass cover. Grid cells covered by less than 35% grasses were excluded.


Supplementary material for this article is available at

Fig. S1. Carbon isotope composition of co-occurring contemporaneous bison and mammoth tissues.

Fig. S2. Individual model predictions for the modern.

Fig. S3. Individual model predictions for the mid-Holocene.

Fig. S4. Individual model predictions for the LGM.

Fig. S5. Model results incorporating ±1σ in mean CMIP5 paleoclimate model output for the mid-Holocene.

Fig. S6. Model results incorporating ±1σ in mean CMIP5 paleoclimate model output for LGM.

Fig. S7. Average summer temperature difference map from LGM to the modern (1901–1950).

Fig. S8. SD maps for LGM COT, GSP, and GST.

Fig. S9. Variability in importance rankings using ±1σ in mean CMIP5 paleoclimate model data.

Table S1. Bison and mammoth isotope data and climate data used in CART analyses.

Table S2. Comparison of hair, bone, and enamel of bison from the same locality.

Table S3. Pearson correlation coefficients for climate variables used in CART analysis.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank the National Museum of Natural History, the American Museum of Natural History, the Harvard Museum of Comparative Zoology, the Berkeley Museum of Vertebrate Zoology, the Alaska Museum of the North, the Denver Museum of Nature and Science, the Field Museum, the University of Michigan Museum of Zoology, the Sam Noble Oklahoma Museum of Natural History, the Yale Peabody Museum, Kansas State University, and the University of Kansas Museum of Natural History for access to bison samples. We also thank K. Chritz and S. Blumenthal for isotope analyses and discussions, S. Harrison for helpful suggestions on an earlier version of the manuscript, and anonymous reviewers who have greatly improved this manuscript. Funding: Funding for this research was provided by the NSF MacroSystems Biology program (grant 1137336). Author contributions: J.M.C., C.J.S., and T.E.C. designed the study. J.M.C. and K.A.H. performed the sample analysis. J.M.C., T.M.M., C.J.S., and T.E.C. performed the data analysis. J.M.C. wrote the manuscript with contributions from all other authors. Competing interests: The authors declare they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article