Climatic Controls on C4 Grassland Distributions During the Neogene: A Model-Data Comparison

Fox, David L.; Pau, Stephanie; Taylor, Lyla; Strömberg, Caroline A. E.; Osborne, Colin P.; Bradshaw, Catherine; Conn, Stephen; Beerling, David J.; Still, Christopher J.

doi:10.3389/fevo.2018.00147

ORIGINAL RESEARCH article

Front. Ecol. Evol., 24 September 2018

Sec. Biogeography and Macroecology

Volume 6 - 2018 | https://doi.org/10.3389/fevo.2018.00147

This article is part of the Research Topic Revisiting The Biome Concept With A Functional Lens View all 10 articles

Climatic Controls on C₄ Grassland Distributions During the Neogene: A Model-Data Comparison

$\r\nDavid L. Fox*$ David L. Fox¹^*

Stephanie Pau²

Lyla Taylor³

Caroline A. E. Strömberg⁴

Colin P. Osborne³

Catherine Bradshaw⁵^†

Stephen Conn⁵^†

David J. Beerling³

Christopher J. Still⁶

¹Department of Earth Sciences, University of Minnesota, Minneapolis, MN, United States
²Department of Geography, Florida State University, Tallahassee, FL, United States
³Department of Animal and Plant Sciences, University of Sheffield, Sheffield, United Kingdom
⁴Department of Biology, University of Washington, Seattle, WA, United States
⁵Bristol Research Initiative for the Dynamic Global Environment, School of Geographical Sciences, University of Bristol, Bristol, United Kingdom
⁶Department of Forest Ecosystems and Society, Oregon State University, Corvallis, OR, United States

Grasslands dominated by taxa using the C₄ photosynthetic pathway first developed on several continents during the Neogene and Quaternary, long after C₄ photosynthesis first evolved among grasses. The histories of these ecosystems are relatively well-documented in the geological record from stable carbon isotope measurements (of fossil vertebrate herbivores and paleosols) and the plant microfossil record (pollen and/or phytolith assemblages). The distinct biogeography and ecophysiology of modern C₃ and C₄ grasses have led to hypotheses explaining the origins of C₄ grasslands in terms of long-term changes in the Earth system, such as increased aridity and decreasing atmospheric pCO₂. However, quantitative proxies for key abiotic drivers of these hypotheses (e.g., temperature, precipitation, pCO₂) are still in development, not yet widely applied at the continental or global scale or throughout the late Cenozoic, and/or remain contentious. Testing these hypotheses globally therefore remains difficult. To understand better the potential links between changes in the Earth system and the origin of C₄ grasslands, we undertook a global-scale comparison between observational records of C₄ plant abundances in Miocene and Pliocene localities compiled from the literature and three increasingly complex models of C₄ physiology, dominance, and abundance. The literature compilation comprises >2,600 δ¹³C-values each of fossil terrestrial vertebrates and of paleosol carbonates, which we interpret as primarily proxies for the abundance of C₄ grasses, based on the modern contribution of C₄ grasses to terrestrial net primary productivity. We forced the vegetation models with simulated monthly climates from the HadCM3 family of coupled ocean-atmosphere general circulation models (OAGCMs) over a range of pCO₂-values for each epoch to model C₄ dominance or abundance in grid cells as: (1) months per year exceeding the temperature at which net carbon assimilation is greater for C₄ than C₃ photosynthesis (crossover temperature model); (2) the number of months per year exceeding the crossover temperature and having sufficient precipitation for growth (≥25 mm/month; Collatz model); and (3) the Sheffield Dynamic Global Vegetation Model (SDGVM), which models multiple plant functional types (PFTs) (C₃ and C₄ grasses, evergreen, and deciduous trees). Model-data agreement is generally weak, although statistically significant for many comparisons, suggesting that regional to local ecological interactions, continent-specific plant evolutionary histories, and/or regional to local climatic conditions not represented in global scale OAGCMs may have been equally strong or stronger in driving the evolution of C₄ grasslands as global changes in the Earth system such as decreases in atmospheric pCO₂ and late Cenozoic global cooling and/or aridification.

Introduction

Plants using the C₄ photosynthetic pathway comprise only a small fraction of the species diversity of the modern global flora, and about half of these extant C₄ species are grasses. Despite the limited number of C₄ species, tropical, and sub-tropical ecosystems dominated by C₄ grasses account for almost one-quarter of modern terrestrial net primary productivity (Still et al., 2003; Sage, 2004; Sage et al., 2011). Molecular phylogenies indicate that the C₄ syndrome evolved multiple times independently in grasses (Kellogg, 1999; GPWGII—Grass Phylogeny Working Group II, 2012), and the calibration of these phylogenies to time using the fossil record suggests that the earliest C₄ grasses had evolved by the Oligocene epoch (33–25 Ma), with successive independent origins of C₄ anatomy and biochemistry continuing throughout the Neogene Period (25–2.6 Ma; Christin et al., 2008, 2014; Spriggs et al., 2014).

The increasingly comprehensive species sampling of grass phylogenies has dramatically improved our ability to generate testable hypotheses about the timing of trait evolution and the diversification of C₄ lineages, but cannot reveal the temporal pattern of increasing biomass and the ecological importance of C₄ grasses in local and regional ecosystems. Many grasslands are dominated locally by only a few species in terms of frequency, abundance, and biomass (Smith and Knapp, 2003; Edwards et al., 2010; Griffith et al., 2015), such that species diversification over time does not necessarily predict the ecological success of C₄ grasses generally. However, the geological record provides evidence about both when grasslands emerged and when they became dominated by C₄ species.

The plant fossil record, particularly phytoliths (microscopic intra- and extracellular silica bodies diagnostic to specific clades of grasses), indicates that grass-dominated habitats first appeared on some continents by the late Oligocene to early Miocene epochs. Complementary evidence based on carbon isotope records from fossil teeth of mammalian herbivores and calcareous paleosols (fossilized soils that preserve soil-derived carbonate) shows that C₄-dominated grasslands first appeared on various continents during the late Miocene to early Pleistocene (Jacobs et al., 1999; Tipple and Pagani, 2007; Edwards et al., 2010; Strömberg, 2011; Fox et al., 2012a). However, various terrestrial isotopic records (e.g., Kleinert and Strecker, 2001), as well as marine records of terrigenous materials (e.g., Feakins et al., 2013; Hoetzel et al., 2013), indicate that the temporal pattern of the emergence of C₄-dominated grasslands was not a simple, monotonic increase in C₄ biomass in some regions. Instead, some records (e.g., the terrigenous signal in marine sediments of northeast African vegetation during the Neogene; Feakins et al., 2013) indicate fluctuations through time in the C₃:C₄ ratio at regional scales prior to the emergence of modern, C₄-dominated grasslands. Thus, the ecological dominance of C₄ grasses in many grassland biomes today post-dates the origins of most C₄ grass clades by many millions of years and reflects a complex history of ecosystem assembly. Despite the attention of a diverse community of systematists, ecologists, paleontologists, and geochemists over many years, the specific combinations of environmental and climatic changes responsible for the ecological success of C₄ grasses over the Neogene and Quaternary both globally and in specific regions have remained elusive and debated (see discussions in, e.g., Edwards et al., 2010; Scheiter et al., 2012; Bond, 2015).

C₄ photosynthesis describes a suite of anatomical and biochemical characteristics that together act as a CO₂ pump to reduce photorespiration, thereby increasing the quantum yield, a measure of how efficiently plants fix CO₂ relative to the amount of absorbed photosynthetically active radiation (APAR) received during light-limited photosynthesis (Ehleringer and Björkman, 1977; Collatz et al., 1998). Photorespiration rates are higher at warmer temperatures for a given atmospheric pCO₂, and at lower pCO₂ for a given temperature (Ehleringer et al., 1997; Sage, 2004). The dominance of modern C₄ grasses in warm to hot and mesic to arid regions (e.g., Teeri and Stowe, 1976; Hattersley, 1983; von Fischer et al., 2008) has led to a long-held assumption that C₄ photosynthesis evolved as an adaptation to warm and dry climates. However, ancestral state reconstructions of the climatic origins of C₄ lineages have called this assumption into question, by showing that most C₄ grass lineages derive from C₃ ancestors that already inhabited warm climates (Edwards and Still, 2008; Edwards and Smith, 2010). From these origins in warm climates, C₄ photosynthesis enabled diversification into cooler and drier climates (Osborne and Freckleton, 2008; Watcharamongkol et al., 2018).

Some paleoclimate records associated with estimates of C₄ biomass, such as compound-specific isotope analyses of terrestrial leaf waxes in Pliocene marine sediments off of West Africa (Kuechler et al., 2018), are consistent with a relationship between aridity and C₄ biomass during a time with atmospheric pCO₂ ≥300 ppmV (i.e., above pre-industrial values; Foster et al., 2017). However, similar records from Pleistocene marine sediments off of North Africa (e.g., Kuechler et al., 2013) indicate an inverse relationship between regional aridity in northern Africa and C₄ abundance against a background of low, but fluctuating glacial-interglacial pCO₂ (see also Urban et al., 2015). These contrasting records suggest that the role of aridity in controlling the abundance of C₄ plants in regional ecosystems may be strongly modulated by factors other than CO₂. These observations are broadly consistent with our understanding of C₄ grass biogeography today. For example, independently evolved C₄ grass lineages show contrasting relationships with rainfall (Visser et al., 2012, 2014). In addition, very arid areas are rarely dominated by C₄ species, whereas particular C₄ grass clades may dominate humid regions if the dry season is sufficient to allow frequent fires (Bond, 2015).

Decreasing atmospheric pCO₂ during the Cenozoic was initially proposed as an explanation for the appearance of C₄-dominated grasslands beginning in the late Miocene (e.g., Cerling et al., 1997; Ehleringer et al., 1997). Subsequent work challenged this assertion, demonstrating that atmospheric pCO₂ first decreased to levels that could favor C₄ over C₃ photosynthesis given prevailing climatic conditions during the Oligocene (Pagani et al., 1999; Beerling and Royer, 2011; Foster et al., 2017), not the late Miocene (although note that some CO₂-proxy records indicates a rise and subsequent fall in pCO₂ during the middle-late Miocene; e.g., Kürschner et al., 2008; Bolton and Stoll, 2013). The broad coincidence in the timing of declining atmospheric pCO₂ and the proposed, earliest evolution of C₄ photosynthesis in grasses has therefore pointed to a role for low pCO₂ in the evolution of C₄ photosynthesis (e.g., Christin et al., 2008; Spriggs et al., 2014). In contrast, the asynchronous appearance of C₄-dominated grasslands in different continents from the late Miocene to the early Pleistocene (Fox and Koch, 2003; Edwards et al., 2010) suggests that additional environmental factors, in combinations specific to each region, were of equal or larger importance for increasing the dominance of C₄ grasses in many ecosystems. These factors could include both biotic and abiotic influences such as disturbance by herbivores, changes in fire frequency and intensity, differential patterns of recovery by C₃ trees and grasses and C₄ grasses linked to traits other than photosynthetic pathway (i.e., evolutionary history of different lineages), soil type and characteristics, and region-specific patterns of climate change during the Neogene and Quaternary (e.g., Osborne, 2008; Edwards et al., 2010; Bond, 2015; Griffith et al., 2015; Charles-Dominique et al., 2016). The importance of declining CO₂ as a single factor in promoting C₄ dominance has also been called into question by the results of a 20 year free-air CO₂ enrichment experiment that showed an eventual enhancement in C₄ biomass relative to C₃ biomass growing under experimentally increased CO₂ (Reich et al., 2018). Water cycle feedbacks induced by enhanced CO₂ have also been shown experimentally to cause C₄ grasses to outperform C₃ grasses, with greater water-use efficiency by C₄ grasses driven by interacting effects of higher CO₂ and warming on both soil moisture and plant photosynthetic responses as the proposed mechanism (Morgan et al., 2011).

Fire in particular is recognized as a critical factor in the abundance of C₄ grasses in many modern ecosystems (Bond, 2015), and fire has been proposed as playing an important role in the emergence of C₄-dominated grasslands during the late Neogene (Keeley and Rundel, 2005; Scheiter et al., 2012). The importance of fire for C₄-dominated grasslands has support from micro-charcoal counts in both Pliocene (Hoetzel et al., 2013) and Holocene (Dupont and Schefuss, 2018) marine sediments from cores off of West Africa, but the terrestrial record of charcoal is still too incomplete (e.g., Scott, 2000) to allow for direct comparisons of charcoal occurrences with either paleobotanical or vertebrate and paleosol carbon isotope proxies documenting the emergence of C₄-dominated ecosystems.

A wide range of continental paleoclimate proxies derived from paleontological and geological materials can provide critical information on the abiotic conditions associated with the increase of C₄ biomass in local and regional ecosystems. Some of these, such as the oxygen isotope composition of fossil teeth of mammalian herbivores (e.g., Passey et al., 2002) or paleosol carbonates (e.g., Fox et al., 2012b) are controlled by both temperature and precipitation, so cannot generally be deconstructed to provide unambiguous quantitative estimates of paleoclimatic conditions. Other proxies, such as the relative tooth crown heights of mammalian herbivores (Eronen et al., 2010) or weathering indices based on elemental ratios of paleosols (Sheldon and Tabor, 2009; Stinchcomb et al., 2016), can provide quantitative estimates of mean annual temperature and precipitation. However, such proxies, as well as newer methods such as carbonate clumped isotope paleothermometry (Ghosh et al., 2006; Eiler, 2011), either rely on assumptions about faunal trait evolution that have been called into question (e.g., MacFadden et al., 1999; Strömberg, 2006), or have not yet been applied widely enough to allow for compilations at the continental or global scale with sufficient coverage to allow for meaningful analyses of empirical paleoclimate data in relation to records of the Neogene-to-Quaternary rise in global C₄ biomass. In contrast, global scale, ocean-atmosphere general circulation models (OAGCMs) coupled with vegetation modeling provide a means to assess the record of C₄ biomass on all continents in relation to possible past climatic conditions that ultimately can be tested with paleoclimate proxy data (e.g., see Lunt et al., 2007 for a comparison of late Oligocene and pre-industrial boundary conditions).

In this paper, we undertake a global-scale comparison between a literature compilation of δ¹³C values of mineralized tissues of fossil herbivorous vertebrates (mammal teeth, ratite eggshell) and of paleosol carbonates from Miocene and Pliocene localities as proxies for the abundance of C₄ grasses, and three increasingly complex models of C₄ dominance and abundance driven by output from a global OAGCM. Specifically, we use climate output from the HadCM3LBL-M2.1 and HadCM3BL-M1 OAGCMs (Valdes et al., 2017) over a range of pCO₂ values for the Miocene (Bradshaw et al., 2012) and Pliocene (Bragg et al., 2012; Conn, 2012) to model C₄ dominance or abundance in grid cells as: (1) months per year exceeding the temperature at which net assimilation is greater for C₄ than C₃ photosynthesis (crossover temperature model; Ehleringer et al., 1997); (2) the number of months per year exceeding the crossover temperature and having sufficient precipitation for grass growth (≥25 mm/month in months that meets the crossover temperature criterion; Collatz model; Collatz et al., 1998), and (3) predictions of various measures of C₄ biomass from the Sheffield Dynamic Global Vegetation Model (SDGVM; Woodward and Lomas, 2004), output from which includes a range of vegetation parameters of multiple plant functional types (PFTs), including C₃ and C₄ grasses. To account for the uncertainty in estimates and temporal variation in CO₂ levels suggested by proxy data for this time interval (see Beerling and Royer, 2011), we used a range of pCO₂ values in each analysis. Statistical comparisons of the isotopic databases with the paleoclimate and vegetation model outputs allow us to assess data-model agreement. We propose that close similarity between proxy data and a particular model will provide support for the explanatory factors inherent in the model playing an important role in C₄ dominance globally (or regionally). A mismatch, on the other hand, will suggest that the factors accounted for in the model are not sufficient for explaining the patterns seen in the data, although alternative reasons are possible. Specifically, we discuss the limitations and opportunities with current data and model resolution for inferring the role of abiotic factors in the evolution of the modern C₄ grasslands over the Neogene and Quaternary.

Carbon Isotope Proxies for C₄ Biomass

The basis for using the stable carbon isotope composition (δ¹³C) of carbonate-bearing minerals of vertebrate herbivores and authigenic soil carbonates as proxies for C₄ biomass is the difference in isotope discrimination by C₃ and C₄ plants during fixation of CO₂. Due to physiological differences, C₃ plants discriminate more strongly against ¹³CO₂ than do C₄ plants (O'Leary, 1981), and, consequently, tissues of C₃ plants have much lower (more negative) δ¹³C values than do C₄ plants. Here, we treat the δ¹³C values of modern plants reported by Passey et al. (2002) as end-members, whereby modern C₃ plants have mean δ¹³C (±1 s.d.) of −27.4 ± 1.6% and modern C₄ monocots have mean δ¹³C of −12.7 ± 1.1%. Passey et al. (2002) assumed an average δ¹³C value for atmospheric CO₂ for these plant data of −8.0%, which allowed estimation of apparent fractionation factors during fixation of CO₂ by C₃ and C₄ photosynthesis of −19.6 and −4.7%, respectively.

Authigenic carbonate precipitates in soils in isotopic equilibrium with the CO₂ that is dissolved in soil water (Cerling, 1984). For soils with moderate to high respiration rates and at modern atmospheric pCO₂, all CO₂ below about 30 cm in the profile is derived from the standing plant biomass or respiration of plant-derived components in the soil, and has the carbon isotopic signature of the standing biomass (Cerling et al., 1991). Relative to atmospheric CO₂, pedogenic carbonate is enriched in ¹³C by ca. 14–17% due to faster diffusion of ¹²CO₂ from the soil to the atmosphere and the temperature-dependent fractionation of carbonate precipitation over the range of typical soil temperatures (Cerling et al., 1991). Paleosol carbonates form on timescales of 10²-10⁵ years, depending on sedimentation rates and landscape position (Gile et al., 1966; Birkeland, 1999); however, they sample the landscape on very localized spatial scales assuming minimal horizontal advection of CO₂ in the subsurface (i.e., by soil water flow). Importantly, vegetation, soil type and properties, and presence or preservation of carbonate can all vary on small spatial scales with landscape position (Zamanian et al., 2016). Thus, the δ¹³C of paleosol carbonates preserves a time-integrated signature of plant biomass that can be used to estimate the C₃:C₄ ratio of the local ecosystem on relatively long timescales but small spatial scales.

The δ¹³C values of vertebrate herbivore tissues record the δ¹³C of diet with tissue-specific fractionation factors or apparent enrichment factors (in the case of tissues such as enamel for which the biosynthetic reactions are not reversible, DeNiro and Epstein, 1978; Cerling and Harris, 1999; Passey et al., 2005). For structurally bound carbonate in the hydroxyapatite of tooth enamel of large-bodied mammalian herbivores, the apparent enrichment factor relative to diet is +14.1 ± 0.5% (Cerling and Harris, 1999). The isotopic composition of mammalian tooth enamel is generally considered to be resistant to alteration during fossilization (Wang and Cerling, 1994; Koch et al., 1997; Zazzo et al., 2004), so that fossil teeth of herbivorous mammals faithfully record the dietary proportions of C₃ and C₄ plants during the interval of tooth formation (Kohn and Cerling, 2002). Additionally, Griffith et al. (2017) showed that the δ¹³C-values of grazer tissue (bison, mammoth) faithfully records the C₃:C₄ ratio of grasslands across larger spatial scales. The carbonate of ratite eggshell also records the δ¹³C of the diet of the hen while the egg is mineralizing (Schirnding et al., 1982), and eggshell δ¹³C values also appear to be resistant to alteration during fossilization (Stern et al., 1994; Miller et al., 2005). Fossil deposits are subject to time averaging, such that data for fossils from individual localities can represent a wide range of ages (Behrensmeyer et al., 2000). However, individual teeth and eggshells only represent a few months to years of mineralization; thus, the temporal resolution of teeth is higher than that for paleosols. On the other hand, the animals did not necessarily live where the fossils accumulated (Behrensmeyer et al., 2000), and mammalian home range generally scales positively with body size (Harestad and Bunnel, 1979). For example, a large-bodied mammal like a moose has a home range of ca. 1,500 ha, but a smaller bodied species such as a peccary has a home range of only ca. 135 ha (data from Harestad and Bunnel, 1979). Teeth of large mammalian herbivores, which have been most commonly used for paleovegetation reconstruction, therefore integrate the landscape C₃:C₄ ratio over a much larger area (up to 100s of km²) over a shorter timescale than do paleosol carbonates.

Clades other than Poaceae (grasses) also include lineages that have independently evolved C₄ photosynthesis and may be ecologically dominant, notably Cyperaceae (sedges) (Sage et al., 2011), and some marine and lacustrine isotopic records of C₄ biomass likely reflect local abundance of C₄ plants other than grasses (e.g., Schefuß et al., 2011; Ivory and Russell, 2016). However, the herbaceous component of the spatially extensive modern tropical and temperate grasslands and savannas are overwhelmingly dominated by C₄ grasses (Still et al., 2003; Lehmann et al., 2011). Relatively few studies have compared carbon isotopic data and plant microfossil assemblages from the same terrestrial sites, but those that have (e.g., McInerney et al., 2011; Strömberg and McInerney, 2011; Cotton et al., 2012; Chen et al., 2015; Smiley et al., 2016) suggest that C₄ isotopic signatures in terrestrial sediments are dominated by grasses and not by other clades, such as sedges, that also evolved C₄ photosynthesis. Based on this reasoning, we assume that the terrestrial isotopic records analyzed herein primarily reflect the abundance of C₄ grasses in Neogene ecosystems.

Carbon Isotope Data Compilations

The dataset of paleosol carbonate δ¹³C-values was compiled from publications and datasets in DLF's personal collection and from extensive (but potentially not exhaustive) bibliographic searches using Google Scholar and the GeoRef bibliographic database maintained by the American Geosciences Institute accessed via the University of Minnesota Libraries during 2011–2013. Search terms included combinations of paleosol, pedogenic, carbon*, isotop*, soil, Neogene, Quaternary, Miocene, Pliocene, and Pleistocene, where * is a wildcard operator appropriate for each database. To be included, published datasets had to include specific geographic information for measured sections or individual samples and radiometric ages and/or precise chronostratigraphic assignments in the Miocene (25.0–4.9 Ma) and/or Pliocene (4.9–2.6 Ma) epochs. Most papers included tables of individual δ¹³C values, but for papers in which data were only presented in figures, Data Thief (https://datathief.org/) was used to capture digitally individual data points from figures. Tests of Data Thief on figures with accompanying data tables verified the accuracy of data capture. The paleosol compilation includes data from 44 publications published from 1986 to 2012 (Data Sheet 1), including 1,280 values from 189 Miocene sections and 1,379 values from 101 Pliocene sections in Africa, Eurasia (including southeast Europe), North America, and South America (Figures 1A, 2).

FIGURE 1

Figure 1. Locations of sample localities for Miocene and Pliocene paleosol carbonates (A) and biominerals (B).

FIGURE 2

Figure 2. Summary of the compilation of paleosol carbonate δ¹³C values. Number of localities (A) and δ¹³C values (B) in each 0.5 Myr bin. All δ¹³C values by mid-point of 0.5 Myr age bin (C) and by latitude of sample site (D). Mean and standard deviation of estimated percent C₄ biomass for each locality by midpoint of 0.5 Myr age bin (E) and by latitude of site (F).

The dataset of fossil mammal tooth and ratite eggshell δ¹³C values was initially compiled by Ben Passey (now at University of Michigan, Ann Arbor, MI, U.S.A.), and was completed with the addition of data from publications and datasets in DLF's personal collection and from extensive bibliographic searches using Google Scholar and GeoRef. Search terms included combinations of mammal*, fossil, tooth, enamel, apatite, bioapatite, carbon*, isotop*, Neogene, Quaternary, Miocene, Pliocene, and Pleistocene. Inclusion of papers and data acquisition followed the criteria for paleosol carbonates. The biomineral compilation includes data from 57 papers published from 1994 to 2013 (Data Sheet 2), including 1,853 values from 173 Miocene faunal sites and 792 values from 52 Pliocene faunal sites in Africa, Eurasia, North America, and South America (Figures 1B, 3). Paleosol carbonate and biomineral data were also compiled for Pleistocene localities (not exhaustively), and these are summarized in Figures 2, 3 but not analyzed here.

FIGURE 3

Figure 3. Summary of the compilation of biomineral δ¹³C values. Number of localities (A) and δ¹³C values (B) in each 0.5 Myr bin. All δ¹³C values by mid-point of 0.5 Myr age bin (C) and by latitude of sample site (D). Mean and standard deviation of estimated percent C₄ biomass for each locality by midpoint of 0.5 Myr age bin (E) and by latitude of site (F).

To compare published δ¹³C values to the vegetation model results, we first had to correct the δ¹³C values for secular change in the δ¹³C of atmospheric CO₂ over the Miocene and Pliocene (Passey et al., 2009; Tipple et al., 2010), which causes the end-member C₃- and C₄-values to vary through time. To do this, we binned the data in both compilations into 0.5 Myr intervals (approximating the Miocene-Pliocene boundary as 5.0 Ma and the end of the Pliocene as 2.5 Ma) and used the estimates of the δ¹³C of atmospheric CO₂ (based on the δ¹³C of benthic foraminifera) in those 0.5 Myr bins to correct the measured δ¹³C values to an equivalent value assuming a “modern” δ¹³C value of −8.0% (following Passey et al., 2002, 2009). The modern mean values for C₃ and C₄ plants, the apparent enrichment factors for each photosynthetic pathway, and the uncertainties for the mean values and enrichment factors served as inputs to the linear mixing model IsoError 1.04 (Phillips and Gregg, 2001), which were used to calculate the mean percentage C₃ and C₄ biomass for each corrected δ¹³C value, as well as the standard deviation for each mean value and the 95% confidence interval for each estimate. The enrichment in ¹³C by C₃ photosynthesis during fixation of CO₂ decreases with increasing irradiance and water deficit (Cernusak et al., 2013), reducing the spacing between C₃ and C₄ end-members, but given the temporal and geographic scope of the analyses here and for the sake of simplicity, we did not consider the potential influence of light availability and water deficits on δ¹³C value of C₃ plants in our interpretations of the carbon isotopic data.

Paleoclimate Simulations

We forced the vegetation models with simulated monthly climates from the HadCM3 family of coupled OAGCMs (Valdes et al., 2017) for Miocene (HADCM3BL-M2.1; Bradshaw et al., 2012) and Pliocene (HADCM3B-M1; Bragg et al., 2012; Conn, 2012, unpublished Master's thesis) paleogeography. Details of these models are given in Valdes et al. (2017). Grid cells are 2.5° latitude by 3.75° longitude and outputs taken from the OAGCM for each grid cell include monthly temperature, precipitation, and humidity. We used two Miocene simulations: a low pCO₂ case (280 ppmV), which is close to the best fit for proxy data during the late Miocene (Foster et al., 2017) and a high pCO₂ case (401 ppmV), which could represent conditions just after the Miocene climatic optimum (ca. 16.5–15 Ma) (e.g., Kürschner et al., 2008). We used three Pliocene simulations: a low pCO₂ case (280 ppmV), a moderate pCO₂ case (405 ppmV, considered here equivalent to the Miocene high pCO₂ case), and a high pCO₂ case (560 ppmV) that represents high pCO₂ forcing of the Mid Pliocene Warm Period (ca. 3.3–3.0 Ma) (e.g., Haywood et al., 2016). The SDGVM averages the climates for neighboring grid cells and restricts humidity to the range 30–95%, and for consistency these averaged climates were used for all three approaches.

Vegetation Models

We used the simulated climates to predict dominance or abundance of C₄ grasses in each grid cell employing three increasingly complex vegetation models. The simplest model is the crossover temperature model. A hallmark of C₄ plants is their dominance in high-light and high-temperature environments such as grasslands and savannas (Long, 1999; Sage et al., 1999). The dominant process-based explanation for environmental controls on C₃ and C₄ grass distributions is the crossover temperature hypothesis, which is based on the different quantum yields of grass species using each pathway (Ehleringer, 1978; Ehleringer et al., 1997). The quantum yield describes how much CO₂ is absorbed by a leaf compared to the amount of APAR it receives during light limited photosynthesis (Ehleringer and Björkman, 1977; Ehleringer, 1978; Collatz et al., 1998). The quantum yield in C₃ plants decreases with increasing leaf temperature, and increases with CO₂ at a given leaf temperature, essentially reflecting the influence of these factors on photorespiration (Ehleringer and Björkman, 1977; Pearcy and Ehleringer, 1983; Collatz et al., 1998; Sage, 2004). By comparison, the quantum yield of C₄ plants is relatively constant across a range of temperatures and CO₂ levels due to the C₄ carbon-concentrating mechanism. The point at which the quantum yield of C₃ grasses equals the quantum yield of C₄ grasses was defined as the “crossover temperature” (Ehleringer et al., 1997; Collatz et al., 1998). At temperatures below this cutoff, C₃ grasses should have a higher capacity to fix carbon, while C₄ grasses should have higher capacities at temperatures above this cutoff. The crossover temperature model simplifies physiological differences between C₃ and C₄ plants but has nevertheless been fairly successful in predicting large-scale distributions (Collatz et al., 1998; Still et al., 2003; Griffith et al., 2015).

For each of the five paleoclimate simulations, we estimated crossover temperatures based on the assumed values of atmospheric pCO₂. Representative photosynthetic fluxes were predicted using the coupled C₃ and C₄ leaf photosynthesis and stomatal conductance models of Collatz et al. (1991, 1992). Parameter values, such as maximum carboxylation rates (Vmax) and temperature response functions, were taken from Sellers et al. (1996). Vmax for C₃ grasses was assumed to be 90 μmol m⁻² s⁻¹ at 298 K, and 30 μmol m⁻² s⁻¹ for C₄ grasses at 298 K. These models, which are described in detail elsewhere, estimate gross photosynthetic rates as a function of temperature, relative humidity, insolation, and the partial pressures of atmospheric CO₂ and O₂. Crossover temperature values were estimated for light-limited conditions and represent the leaf temperature where C₃ and C₄ photosynthetic rates were equal. Crossover temperatures calculated for the 280, 401/405, and 560 ppmV cases are 14, 20, and 24°C, respectively.

For each grid cell in each paleoclimate simulation, the crossover temperature model as used here is simply the number of months for which the temperature is equal to or greater than the calculated crossover temperature (based on the assumed pCO₂ in each OAGCM simulation), which is an estimate of the number of months that production of C₄ grasses should be favored, all else being equal. Our use of the crossover temperature concept is intentionally simplistic as a first step in estimating the dominance or abundance of C₄ grasses.

Our second vegetation model, the Collatz model (Collatz et al., 1998), is slightly more complex, and includes a more reasonable modern climate threshold for C₄ grass dominance that accounts for both crossover temperature (22°C for contemporaneous pCO₂) and the moisture necessary for plant growth (≥25 mm of precipitation in a month that meets the crossover temperature criterion). Thus, for each grid cell in each paleoclimate simulation, the Collatz model determines the number of months for which the modeled temperature is equal to or greater than the calculated crossover temperature and for which precipitation is also equal to or >25 mm. The inclusion of a precipitation screen effectively limits the months considered to the growing season.

Finally, we used the Sheffield Dynamic Global Vegetation Model (SDGVM; Woodward and Lomas, 2004) to estimate the proportion of C₄ biomass in each OAGCM grid cell. Inputs for the SDGVM are monthly temperature, precipitation, and humidity (taken from each OAGCM simulation) and soil texture and dynamics. The SDGVM, which has been used extensively to model both modern and past vegetation (e.g., Beerling and Woodward, 2001; Bond et al., 2005; Scheiter et al., 2012), simulates photosynthesis, respiration, and transpiration for six natural PFTs, and soil carbon and nitrogen dynamics. PFTs include broadleaf deciduous trees, needleleaf deciduous trees, broadleaf evergreen trees, needleleaf evergreen trees, C₃ grasses, and C₄ grasses. Annual C₃ and C₄ grasses are the primary successional PFTs. If the climate permits, the four woody PFTs may acquire land cover at the expense of grasses over decadal and longer timescales, but trees are killed by several processes, including fire. At the beginning of each model year, the probability of fire determines the fraction of land which will be available for grass encroachment at any gridpoint. The SDGVM then assigns cover to C₃ and C₄ grasses based on their relative NPP, without explicit reference to a crossover temperature (Nemani and Running, 1996). With regard to this study, limitations of the SDGVM include the lack of savanna-specific PFTs and a detailed treatment of fire and herbivory. Although NPP and leaf area index are calculated daily and outputted monthly in the model, biomass is updated only at the end of the year, restricting our ability to examine seasonal root and stem growth. Here, we use average annual leaf area index for each PFT, which scales with biomass, to estimate C₄ grass biomass as a percentage of total simulated plant biomass in each grid cell. This third vegetation model represents a considerable increase in complexity relative to the crossover temperature and Collatz models as used here.

To compare the δ¹³C data to the model estimates of C₄ dominance or abundance, we used statistical summaries of all Miocene and Pliocene δ¹³C values for fossil biominerals and for paleosol carbonates separately in each grid cell. Given the global distribution of Miocene and Pliocene continental rock units and associated terrestrial vertebrate faunas, most continental grid cells in the model domains did not have carbon isotope records from Miocene or Pliocene sites, so by necessity only a subset of model grid cells are analyzed in each time interval. When aggregated at the geographic scale of grid cells, the mean and median δ¹³C values are virtually indistinguishable and yield similar results, so below we present only results for mean values. Because the observed and modeled results are not necessarily normally distributed, we used the non-parametric Kendall rank correlation coefficient (Kendall's tau, τ) to test for statistically significant correlations between % C₄ biomass based on mean δ¹³C values in grid cells and the modeled estimates of C₄ dominance (number of months that C₄ is favored) and biomass (SDGVM output). Kendall's τ ranges from +1 for perfect positive correlation of rank orders of the data to −1 for perfect inverse correlation of ranks, and a τ of 0 indicates complete lack of correlation of ranks. Kendall's τ was calculated using the Stats package in R (R Development Core Team, 2013) and in JMP Pro 13.

Results

Comparison of the percentage of C₄ biomass estimated from paleosol carbonate δ¹³C values with the number of months favoring C₄ photosynthesis based on the crossover temperature model is presented in Figure 4. Kendall's rank correlation coefficients (τ) for the two Miocene models are positive and low, but the correlation for the 401 ppmV case is statistically significant (p = 0.016). Kendall's τ for the Pliocene models are all statistically significant, and the correlations are stronger than for the Miocene models, but they are all negative, indicating that grid cells with more months above the modeled crossover temperature have generally lower δ¹³C values for paleosol carbonates, suggesting less C₄ biomass. Comparison of the results for the biomineral data with the crossover temperature model results are presented in Figure 5. Kendall's τ values are positive for all models, and most are statistically significant. The rank correlations are generally stronger than for the paleosol carbonate results, with τ values up to 0.395.

FIGURE 4

Figure 4. Comparisons of results of the crossover temperature model to mean estimated percent C₄ biomass from the paleosol carbonate compilation. (A) Pliocene, 280 ppmV CO₂. (B) Pliocene, 405 ppvM CO₂. (C) Pliocene, 560 ppmV CO₂. (D) Miocene, 280 ppmv CO₂. (E) Miocene, 401 ppmV CO₂.

FIGURE 5

Figure 5. Comparisons of results of the crossover temperature model to mean estimated percent C₄ biomass from the biomineral compilation. (A) Pliocene, 280 ppmV CO₂. (B) Pliocene, 405 ppvM CO₂. (C) Pliocene, 560 ppmV CO₂. (D) Miocene, 280 ppmv CO₂. (E) Miocene, 401 ppmV CO₂.

Comparison of the paleosol carbonate data with the Collatz model are presented in Figure 6. The results are similar to those for the crossover temperature model. The Miocene cases have positive τ values that are statistically significant, but the Pliocene cases have statistically significant, negative correlations between percentage of C₄ biomass estimated from the paleosol carbonate δ¹³C values and the number of months that favor C₄ photosynthesis. Comparison of the results for the biomineral data with the Collatz model results are presented in Figure 7 and are generally similar to the comparison of the biomineral data with the crossover temperature model (Figure 5): all τ values are positive and several are statistically significant, but the rank correlations are generally weak with most having τ values of 0.062 to 0.200.

FIGURE 6

Figure 6. Comparisons of results of the Collatz model to mean estimated percent C₄ biomass from the paleosol carbonate compilation. (A) Pliocene, 280 ppmV CO₂. (B) Pliocene, 405 ppvM CO₂. (C) Pliocene, 560 ppmV CO₂. (D) Miocene, 280 ppmv CO₂. (E) Miocene, 401 ppmV CO₂.

FIGURE 7

Figure 7. Comparisons of results of the Collatz model to mean estimated percent C₄ biomass from the biomineral compilation. (A) Pliocene, 280 ppmV CO₂. (B) Pliocene, 405 ppvM CO₂. (C) Pliocene, 560 ppmV CO₂. (D) Miocene, 280 ppmv CO₂. (E) Miocene, 401 ppmV CO₂.

Comparisons of the paleosol carbonate data with the SDGVM output for the two Miocene and three Pliocene cases are presented in Tables 1, 2, respectively. For both Miocene cases, all correlation coefficients are positive and many are statistically significant, although generally the correlations are weak (τ ranges from 0.052 to 0.322 for the 280 ppmV case and from 0.053 to 0.264 for the 401 ppmV case). For both sets of results, the strongest correlations are for the percentage of total annual biomass that is C₄ grasses (highest τ for both models) and percent C₄ cover. For the Pliocene cases, many correlations are statistically significant, but negative, and the correlations are generally weak.

TABLE 1

Table 1. Kendall's rank correlation coefficients (τ) and p-values for comparisons of estimates of mean percent C₄ biomass from δ¹³C-values of Miocene paleosol carbonate with results of two cases of the SDGVM.

TABLE 2

Table 2. Kendall's rank correlation coefficients (τ) and p-values for comparisons of estimates of mean percent C₄ biomass from δ¹³C-values of Pliocene paleosol carbonate with results of three cases of the SDGVM.

Comparisons of the biomineral data with the SDGVM output for the two Miocene and three Pliocene cases are presented in Tables 3, 4, respectively. For both Miocene cases, all correlations are positive and many are statistically significant, particularly for the 401 ppmV case. The correlations are generally weak, as for the paleosol carbonate data, and the maximum correlation coefficients are slightly lower than for the Miocene paleosol carbonates (τ of ca. 0.22–0.27 for percentage C₄ grass biomass and cover). The rank correlations of the biomineral estimates of C₄ biomass for the three Pliocene cases are the strongest of any of our comparisons, and most are statistically significant. The correlations are negative for three of the SDGVM output parameters for the 280 ppmV case, and otherwise all correlations are positive. The strongest rank correlations are moderately strong (τ ca. 0.4 for all three cases for percentage C₄ grass biomass and cover), particularly compared to most of the other results. It is also noteworthy that these significant, positive τ values are only for C₄-specific model outputs, whereas summary values for all PFTs such as total annual biomass do not show significant positive correlations with proxy data.

TABLE 3

Table 3. Kendall's rank correlation coefficients (τ) and p-values for comparisons of estimates of mean percent C₄ biomass from δ¹³C-values of Miocene biominerals with results of two cases of the SDGVM.

TABLE 4

Table 4. Kendall's rank correlation coefficients (τ) and p-values for comparisons of estimates of mean percent C₄ biomass from δ¹³C-values of Pliocene biominerals with results of three cases of the SDGVM.

Discussion

Our analyses indicate that none of the models have strong explanatory power for the paleo-proxy data at a global scale, with rank correlations reaching a maximum of ca. 0.4. Thus, the degree of process model complexity does not appear to matter substantially for better hindcasting of C₄ grass dominance, at least not when the data from different regions are analyzed together. In other words, it is not clear which of temperature, length of growing season, and ecological interactions (as variously simulated in the three models) mattered most for Neogene C₄ grass distribution and which other factors may have been equally or more important.

Nevertheless, some general conclusions can be drawn. First, many of the comparisons for which ranks are positively correlated are indeed statistically significant, and these results necessarily imply non-random associations between the measures of C₄ dominance and abundance and the predictions based on the three models. In other words, they provide evidence that climate, as the driver of our vegetation models, played a role in the distribution and abundance of C₄ plants (primarily grasses) during the Neogene. However, most of the positive correlations, including those that are statistically significant, are weak to moderate (all τ-values are <0.44 and most for the Miocene are <0.20). These results indicate that factors other than climate and atmospheric pCO₂ were as or more important controls on the dominance and abundance of C₄ plants during the Neogene, at least to the extent that our various models capture climate and atmospheric CO₂ impacts on vegetation physiology and dynamics. That the controls on the distribution of C₄ plants are multifaceted is well-understood. However, climate and atmospheric pCO₂ are the two parameters for which we have the most reliable methods of quantitative estimation in the geological record. Other factors, such as fire frequency and intensity and herbivore disturbance, cannot yet be estimated quantitatively in meaningful ways for most terrestrial sites. Thus, our results highlight the current limits on our ability to make sense of the isotopic proxy record of C₄ abundance in the past using our understanding of modern processes.

Recent studies in modern ecosystems have revealed a broad, but imperfect fit between modeled vegetation output and (proxy) vegetation data in terms of C₄ biomass and dominance (e.g., Ehleringer et al., 1997; Collatz et al., 1998; Cramer et al., 2001; Sitch et al., 2003; Still et al., 2003; Lehmann et al., 2011; Ardö, 2015). Although climate undoubtedly plays a major role in driving the balance between (C₄) grasses and trees in low- to mid-latitude ecosystems (e.g., Sankaran et al., 2005; Hirota et al., 2011; Staver et al., 2011), differential responses to climate in the distribution of savannas across continents point to the role of both disturbance (fire, herbivory) and ecosystem history (e.g., biogeography, soil type) in shaping vegetation dynamics (e.g., Bond, 2008; Lehmann et al., 2011). In light of these at least generally consistent results, the agreement between the carbon isotope data compilations and the vegetation model estimates of C₄ dominance and abundance is poor despite the frequency of positive and statistically significant results. For many of the cases, all three models predict strong dominance by, or high abundance of, C₄ grasses in many grid cells for which the isotopic data suggest only modest C₄ biomass. Additionally, in many cases the models predict essentially no C₄ biomass for many grid cells, for which estimates based on isotopic data range from 0% to almost 80% mean C₄ biomass. We used the average end member δ¹³C values for C₃ and C₄ plants to estimate the percentage of C₄ biomass, which could lead to overestimates of C₄ biomass from sites which experienced high irradiance and/or water deficits. Those conditions lead to a net decrease in carbon isotope fractionation during fixation of CO₂ by C₃ plants (Cernusak et al., 2013), which decreases the isotopic offset between C₃ and C₄ plants. However, even if we adjusted some or all values to account for this possible effect, the decrease in the abundance of C₄ biomass would be no >20% and many grid cells modeled as having no C₄ grasses would still have mean observed δ¹³C values consistent with 50–60% C₄ biomass.

Another important observation from our study is that the results for the biomineral data (predominantly from fossil mammalian herbivore tooth enamel) are stronger than for the paleosol data in terms of both frequency of statistical significance, and the strength and direction of the rank correlations. For example, many of the results for all three vegetation models yield negative correlations with the estimates of C₄ biomass from the paleosol carbon isotope data, and many of these cases are statistically significant, which is the opposite of expectations. Most striking is the uniformity of negative correlations for all of the comparisons with the Pliocene paleosol carbonate data, including for the most sophisticated of our vegetation models, the SDGVM. One explanation for the stronger correlations with biomineral isotopic data is the different biases inherent in the two types of proxy data we analyzed. The isotopic signal in vertebrate herbivores integrates vegetation growing across the landscape, whereas paleosols are essentially time-integrated point samples (Griffith et al., 2017). It may also be that consumption of aboveground biomass by herbivores better mirrors the variation in total C₄ biomass than paleosol carbonates, the carbon for which derives mostly from soil- and root-respired CO₂, thus more closely reflecting underground biomass, as in modern soils (Angelo and Pau, 2016). However, the fact that only the Pliocene data show contrasting correlations for the different types of isotopic data suggests that broad formational factors are not sufficient as an explanation.

Sources of Uncertainty in Model-Proxy Data Comparisons

Despite these weak and partly contradictory results, we see no reason to question the utility of carbon isotopic data for inferring the abundance of C₄ plants in soil biomass and aspects of the diet of primary consumers (see discussion in Cerling and Quade, 1993; Kohn and Cerling, 2002; Cerling et al., 2010). However, interpreting either type of isotopic data in terms of landscape scale patterns of C₄ abundance is perhaps not as straightforward as is often assumed. As discussed above, isotopic data from paleosol carbonates are likely to be much less spatially averaged (and more temporally averaged) than isotopic data from vertebrate teeth. Similarly, isotopic data from herbivore teeth suffer from a number of unique biases. For example, carbon isotope data for consumers for a single locality are subject to many filters that blur or obscure how the distribution of diets in terms of C₄ consumption by individuals and averages for species relate to the abundance of C₄ grasses locally or regionally. First, herbivores may not sample the vegetation evenly or completely, so that even if the fauna could be analyzed in its entirety, it would not provide an unbiased isotopic picture. Second, taphonomic factors can alter the taxonomic composition of faunas, both in terms of presence or absence of species in the fossil record and also relative abundance and evenness of species relative to the original living faunal community (Behrensmeyer et al., 2000). Sampling of fossil taxa rarely accounts for the relative abundance of species in the fossil assemblage, which could be an important consideration for interpreting distributions of δ¹³C values or species or faunal mean values in relation to the inferred nature of the ancient habitat. Third, in modern ecosystems in which essentially all grasses are C₄ and all shrubs and trees are C₃, such as much of low elevation east African savannas today, using δ¹³C values to classify consumers into distinct dietary categories (C₄ grazer, mixed feeder, C₃ browser, closed canopy C₃ browser) and then to classify faunas into proportions of species in each category, as Cerling et al. (2015) did, is a possible way to convert consumer isotopic data into quantitative models of biome type. However, in regions that have ecologically meaningful abundance of C₃ grasses, this approach may not work as well because some grazers will have a C₃ signal. These filters and biases mean that, for application to fossil faunas, careful consideration of feeding ecology, taphonomy, and sampling are necessary. Nevertheless, the isotopic proxies for C₄ are reliable in principle, if somewhat more nuanced than they are commonly treated, so we do not think that the isotopic data themselves are necessarily problematic.

Similarly, we do not question the basic concepts of the crossover temperature or the modern climatic threshold proposed by Collatz et al. (1998) as predictors of the physiological advantage of C₄ plants under the right environmental conditions. The fact that neither our simple crossover temperature model nor the only slightly more complex Collatz model performs well suggests that pCO₂, temperature, and moisture availability alone are not adequate to explain the history of C₄ grasses, despite the fact that the Collatz model does reasonably well at predicting the modern distribution of C₃, C₄, and mixed grasslands today (Collatz et al., 1998). Other local to regional biotic and abiotic factors on each continent must also have been important—and continue to be important—as has been suggested before (Fox and Koch, 2003; Edwards et al., 2010; Lehmann et al., 2014). Even the SDGVM results are not particularly compelling, despite the attempt in this model to account for a broader range of environmental variables and a more nuanced set of plant life strategies, and a demonstrated ability to reasonably recreate modern vegetation patterns (Cramer et al., 2001; Woodward et al., 2001).

Instead, we argue that an important reason for the relatively weak correlations between model output and proxy-data could be the difference in temporal and spatial scale and resolution between the paleosol and consumer records on the one hand and the OAGCM grid cells on the other. Individual paleosol samples are essentially point samples and at the scale of the OAGCM grid cells even multiple stratigraphic series of samples from a field area are still essentially point samples. The few cases of extensive lateral sampling at the outcrop and field area scale have shown considerable heterogeneity in the abundance of C₄ grasses on length scales of 10s to 1,000s of meters that reflect the influence of landscape processes, even in the context of regional transitions from C₃ to C₄ dominated ecosystems (e.g., Behrensmeyer et al., 2007 for the classic Siwaliks record of Pakistan). Landscape position, local topography, and edaphic factors such as grain size and clay content are critical controls on soil drainage and the potential for carbonate precipitation in addition to climate (Birkeland, 1999), so the paleosol record may not be representative of the average conditions across entire OAGCM grid cells. Vertebrate consumers collectively sample the landscape much more broadly than do soils, but only some of the largest species might have had individual home ranges that approach the scale of the model grid cells and we presume that relatively few species had long distance migration routes, as is the case for most modern mammals (Harris et al., 2009). In addition, as mentioned above, even these large herbivores are unlikely to have sampled the landscape in a perfectly representative way. Given that the consumer data reflect a mix of the feeding preferences of individuals and the typical behavior of populations and species, it is therefore possible that the biomineral dataset also does not completely reflect average conditions at the scale of every OAGCM cell (although our strongest results are for the Pliocene biomineral data in relation to the SDGVM output). Thus, the proxies may be sampling small parts of the area of a grid cell relatively accurately, but those may not be representative of the entire area of the grid cell, while the SDGVM results may also be relatively accurately predicting average conditions, but at the scale of the grid cell. Spatial averaging of non-linear physiological responses to temperature and CO₂ (“Jensen's inequality”; Denny, 2017), would result in over- or under-estimation of C₄ dominance and abundance in the large grid cells of the SDGVM.

The temporal scales of our comparisons may also not be commensurate. The isotopic datasets both span millions of years and are unevenly distributed in time and space. Moreover, in some areas (e.g., south Asia), the isotopic record in either the Miocene or the Pliocene or both has known patterns of temporal variability. The OAGCM outputs for the late Miocene and the Pliocene are single time slices that necessarily represent climate over millions of years for each interval and use geographic boundary conditions in the Miocene and Pliocene with their own uncertainties. Furthermore, global and, more importantly, regional climates were likely more dynamic during both periods in ways that could have been important for the evolution of C₄ grasslands. For the Pliocene, the multiple models at different pCO₂ may capture some aspects of the global scale dynamic of the transition into and out of the Mid Pliocene Warm Period (see Haywood et al., 2016). Also, given the relative brevity of the Pliocene compared to the Miocene, the geographic boundary conditions during the Pliocene may not have changed enough to matter. However, that is probably not the case for the Miocene given the histories of active tectonism in western North America, east Africa, and south Asia during this time, and the known influence of changes in orography on climate. These considerations highlight the fact that the paleoclimate simulations that underlie our vegetation models may be inaccurate in specific grid cells or even over regions in ways that contribute to the relatively low correlations between data and models that we observe.

Similarly, despite the comparative sophistication of the SDGVM, various factors may complicate its ability to model ancient ecosystems in a ways that limit the comparisons we have made here, particularly with regard to how grass-dominated ecosystems are modeled. The SDGVM does not have savanna-specific FTs such as fire-adapted, shade-intolerant trees, or a distinct shrub FT. The treatment of fire in this implementation of SDGVM is likely another important factor, given the prominence of fire as a control on modern grass-dominated ecosystems (e.g., Bond et al., 2005). Fire probability is calculated at the beginning of the year using climate envelopes and fuel loads depend on biomass, which is updated annually, so the model is relatively insensitive to vegetation dynamics; in addition, fire intensity is not modeled. The PFTs do not differ in flammability and have no fire-specific adaptations such as thick bark. The treatment of fire may be a reason in some cases for the SDGVM under-predicting C₄ biomass in areas that the isotopic data indicate abundant C₄ biomass. The model does not include herbivory, which like fire is a vital factor in maintaining grassy biomes, and has been implicated in the origins of savannas in Africa (Charles-Dominique et al., 2016). In the implementation we used, climate envelopes directly control tree distribution but only indirectly control grasses via their net primary production. The SDGVM does not have competition between the PFTs, except for crude competition for water resources, but without hydraulic redistribution. These points are not intended as criticisms of the SDGVM, since all models make simplifications. However, these factors may limit how well the SDGVM can model the types of biomes that were dominated by C₄ grasses today and in the past. This is particularly true if a non-analog combination of factors (biotic and abiotic) promoted the spread of C₄-dominated vegetation during the Neogene.

Finally, we note that all three vegetation models are underlain by the simulated climates, which have uncertainties that may be equally high or even higher than those of the SDGVM. Indeed, the climate simulated can be highly dependent on the geographical boundary conditions assumed (e.g., Zhang et al., 2012; Brierley and Fedorov, 2016; Ahlström et al., 2017). Validation of the climates simulated by the OAGCM are complicated by the same limitations of the paleoclimate proxy record that preclude quantitative analysis of the abiotic factors that controlled the history of C₄ grasslands globally over the Neogene. However, numerous model-proxy data comparisons for the Cenozoic attest to the inability of OAGCMs to reproduce climatic events in the past (e.g., Herold et al., 2010; Lunt et al., 2012), with the Paleocene-Eocene Thermal Maximum as a particularly well-studied example (McInerney and Wing, 2011). These model-data mismatches indicate fundamental gaps in our understanding of Earth's climate system that have to be filled before more meaningful comparisons can be made between simulations of past climates and proxy data for past vegetation.

Conclusions

In the end, given all of the possible complications in our data-model comparisons, as well as the complexity of the regional and local processes involved in the evolution of C₄ grasslands over the Neogene (and Quaternary), perhaps the statistically significant comparisons with rank correlation coefficients in the 0.2–0.4 range should actually be viewed as good agreement with the models, particularly the SDGVM. The models are intentional simplifications and abstractions of complex interaction properties and processes, and the historical patterns in the carbon isotopic (and other paleontological and geological) records are complicated in both space and time. Our results are consistent with prior findings that global changes in pCO₂ and global to regional changes in temperature and precipitation alone are not sufficient to explain the histories of C₄-dominated ecosystems since the beginning of the Neogene. Our results also highlight the need for means to estimate in the geological record other factors that are known from modern studies to be important controls on C₄ abundance, such as fire intensity and herbivore disturbance. Fire and herbivory should lead to model under-predictions of C₄ abundance, whereas relying only on climate and CO₂ might lead to over-predictions of C₄ abundance. Since abundance of C₄ grasses depends on disturbance factors that limit tree growth, then model omission or over simplistic representation of herbivory and fire can lead to an underestimation of C₄ grass dominance. Similarly, if C₄ grass dominance is simulated solely on the basis of physiological differences with C₃ plants (the crossover and Collatz models) without considering the ecological interactions between grasses and trees, then models will tend to overestimate C₄ grass abundance. Those aspects of the SDGVM that may complicate its ability to predict grassy biomes suggest that additional model complexity may also be necessary, and regional scale models with higher spatial resolution for multiple time steps through the Neogene and Quaternary are needed to address the scale issues we identified here. However, increasing model complexity and resolution may not yield better model-data comparisons if local ecological interactions and evolutionary history have equally strong effects as abiotic conditions, so that the physiological responses to environment are only one factor among several. Either way, our results should add caution to the reliance on current models for predicting the future of C₄ grasslands, and point to the importance of considering scale when comparing different sources of both ecological and paleoecological data.

Author Contributions

DF, SP, CS, CO, and CJS designed the project. DF compiled the isotopic database, performed the statistical analyses, and wrote the initial draft of the paper. CS calculated the crossover temperatures and Collatz model values. CB and SC performed the climate simulations. LT and DB performed the SDGVM simulations. DF, CS, SP, CO, LT, and CS edited and revised the paper.

Funding

A grant from the National Evolutionary Synthesis Center (NESCent; funded by NSF EF0905606) to CO, CAES, and CJS funded the Origins of C₄ grasslands: a new synthesis of phylogeny, ecology and paleobiology working group that led to this work.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Ben Passey for sharing his initial work on the biomineral isotope data compilation. We also thank the other members of the NESCent C₄ Grasslands Working Group for valuable discussions that helped shaped this project: T. Michael Anderson, William Bond, Erika Edwards, Elisabeth Forrestel, William Hoffmann, Ben Passey, Nicolas Salamin, and Melinda Smith.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2018.00147/full#supplementary-material

Data Sheet 1. Palesol carbonate stable isotope data compilation.

Data Sheet 2. Biomineral stable isotope data compilation.

References

Ahlström, A., Schurgers, G., and Smith, B. (2017). The large influence of climate model bias on terrestrial carbon cycle simulations. Environ. Res. Lett. 12:014004. doi: 10.1088/1748-9326/12/1/014004