Introduction

Greenhouse gas emissions resulting from human activity are considered to be an important driver of the current global warming. Accurate understanding of this cause-effect relationship is essential to propose realistic scenarios of the future evolution of Earth’s climate, but the human monitoring history of these variables may be very tenuous to properly take into account the relative importance of their large time-scale dynamics into the climate system. Geological records offer an opportunity to address this time-scale issue as they provide track of the long-term dynamics of atmospheric greenhouse gas contents and surface temperature throughout Earth’s history. An analysis of this relationship throughout the whole Phanerozoic pointed out atmospheric pCO2 as the main multi-scale driver of global warming1. However, recent improvement and refining of methods to reconstruct pCO2 evidence that this conclusion was based on outdated imprecise estimates2, 3. Additionally, some factors involved in long-term climate dynamics such as O2 have been put forward to challenge this paradigm4. This raises questions about the extent to which we are correctly interpreting the relative influence of atmospheric pCO2 forcing on climate5.

The analysis of the relationship between pCO2 and temperature for the Phanerozoic as a whole may not be pertinent because it comprised shifts between drastically different climate states of the Earth in which the relative impact and the dynamics of atmospheric CO2 may have strongly differed. Of these, Greenhouse is considered to be the default planetary climate state of the Earth, prevailing in more than 70% of the Phanerozoic6. The Cretaceous is one of the longest and most studied Greenhouse periods of Earth’s history, with an extensive background of pCO2 7,8,9,10 and temperature11,12,13 reconstructions. This makes the Cretaceous an exceptional ‘laboratory’ of the past to explore the relationship between the two factors. Nevertheless, pCO2 reconstructions made so far for the Cretaceous defined temporal trends that are not consistent with those described for temperature11,12,13 or the carbon cycle dynamics evidenced by carbon isotopes14,15,16. This raises questions about how accurate are pCO2 reconstructions made so far and how significant was the relative contribution of pCO2 to temperature evolution during the Cretaceous.

Carbon isotopes are relevant tools for describing past carbon cycle and CO2 dynamics14,15,16,17,18,19,20. Carbon isotope compositions of terrestrial plants are particularly interesting as they are strongly linked to atmospheric CO2, which is the main source of carbon assimilated by plants via photosynthesis. The recent construction of models relating carbon isotope fractionation by plants and pCO2 19,20,21 coupled with the good preservation of the pristine carbon isotope signal within plant fossil materials22,23,24,25, make carbon isotope compositions of plant fossils valuable proxies to analyze the evolution of atmospheric pCO2 during past periods.

Herein the long-term evolution of atmospheric pCO2 during the Cretaceous is reconstructed based on the carbon isotope compositions of leaves of the fossil conifer Frenelopsis13Cleaf). pCO2 is estimated from twelve stages spanning the Barremian–Santonian interval (129.4–83.6 Ma) during which major disturbance events of the carbon cycle related to CO2 release to the atmosphere14, 15, 18 and major variations in temperature11,12,13 described for the Cretaceous took place. We explore the correspondence between estimated pCO2 trends and those previously published for temperature to analyze the relative influence of pCO2 forcing on climate during the Cretaceous.

Results and Discussion

Carbon isotope compositions of the fossil conifer Frenelopsis

δ13Cleaf values from the twelve studied stages range from −27.7 to −22.1‰ in average (Fig. 1; Table 1). In most of these stages, δ13Cleaf values show similar standard deviations (ca. 0.9‰), thus reflecting a similar range of local environmental variability. Santonian stages are the only exceptions, which may reflect higher climatic instability related to the development of Oceanic Anoxic Event (OAE) 3. Carbon isotope compositions of Frenelopsis in the Barremian, upper Albian and lower Santonian stages are notably 13C-depleted compared to others, even though it is not significant in the latter case (Fig. 1). These lower δ13Cleaf values are consistent with decreases in the carbon isotope composition of atmospheric CO213CCO2) recently described by Barral et al.16 (Fig. 1). Overall, there is a good agreement between δ13Cleaf and δ13CCO2 curves, though trends in the former are generally magnified. This magnification is more pronounced during the Cenomanian-Turonian interval, coinciding with a critical increase in global temperatures11, which may have slightly reduced carbon isotope fractionation by plants (Δ13Cleaf) by reducing the discrimination effect of RuBisCO against 13C26, 27.

Figure 1
figure 1

Comparison of the carbon isotope compositions of leaves of the fossil conifer Frenelopsis13Cleaf) obtained for each studied episode with the evolution of the carbon isotope composition of atmospheric CO213CCO2) throughout the Cretaceous estimated by Barral et al.16. Envelopes and bar errors correspond to ±1σ.

Table 1 δ13Cleaf, δ13CCO2, Δ13Cleaf and pCO2 values for the twelve studied Cretaceous stages. μ and σ correspond to mean and standard deviation, respectively.

Atmospheric CO2 concentration during the Cretaceous

Δ13Cleaf values inferred from δ13Cleaf and concomitant δ13CCO2 values (Table 1) were used to estimate pCO2 values. pCO2 estimates are in the range of ca. 150–650 ppm for the Barremian–Santonian interval, and the resolution is better than 55 ppm (i.e. 1σ, standard deviation) for most of the stages (Fig. 2; Table 1). This range of values is consistent with the conclusions by Franks et al.20 from a method based on universal equations of leaf gas exchange using stomatal anatomy and carbon isotope ratios of fossil leaves suggesting pCO2 values below 1000 ppm during most of the Phanerozoic, although their conclusions remain controversial28,29,30. Our pCO2 estimates average ca. 300 ppm for the studied interval, which is in accordance with the most recent revision of GEOCARBSULF model3 and the range of values obtained in most recent works dealing with proxy-based reconstructions7,8,9,10 (Fig. 2). All these reconstructions suggest that atmospheric CO2 concentrations were comparable to those of present-day atmosphere, which is a conclusion in deep contrast with the classical view that the Cretaceous period was characterized by atmospheric CO2 concentrations considerably higher than today31. This is also consistent with several other approaches. Boron isotope compositions (δ11B) of marine biogenic carbonates are closely related to pCO2 and pH32,33,34. Higher atmospheric concentrations would imply more CO2 dissolution into the ocean, causing a lowering of seawater pH and δ11B of marine biogenic carbonates33, 34. However, δ11B of Cretaceous and present-day brachiopods as well as ophiolitic serpentinites are close to each other35, thus suggesting that seawater pH and atmospheric pCO2 were of similar magnitude. Analysis of reef abundance and diversity throughout the Phanerozoic also agrees with our results, indicating that during the Cretaceous no significant long-term rise in ocean acidification, reflecting a long-term rise in pCO2, took place36. Only discrete short-term events of ocean acidification have been described during the Cretaceous based on marine calcifiers, associated with major Cretaceous OAEs and Hothouse episodes37, 38.

Figure 2
figure 2

Comparison of the δ13Cleaf-based pCO2 estimates with most recent proxy-based and model-based pCO2 reconstructions for the Cretaceous3, 7,8,9,10. Envelopes and bar errors represent ±1σ. GCSM: GEOCARBSULF Model reconstructions; PB: Paleosol Barometer reconstructions; SI: Stomatal Index reconstructions.

Was pCO2 a long-term driver of temperature during the Cretaceous?

Several events of volcanic CO2 release into the atmosphere linked to the emplacement of large igneous provinces have been inferred from pronounced negative shifts of both carbon and strontium isotope compositions of marine carbonates14, 38 associated with the main Cretaceous OAEs14, 18, 37. These CO2 release events are also evident regarding the evolution of δ13CCO2 throughout the Cretaceous16 (Fig. 1) and are most likely responsible for temperature increases at the kyr scale12, 13. The biggest of these CO2 release events in terms of the magnitude of carbon isotope excursions was related to OAE1b and lasted ca. 25 kyr14, 39. This event was responsible for a fast warming of the atmosphere of ca. 0.3 °C and followed by cooler climatic conditions during a period of ca. 45 kyr39. Greater kyr-scale increases in temperature of ca. 6–7 °C have been described associated with the CO2 release events related to OAE1a12, 13 and OAE211 (Fig. 3). However, our results indicate that drastic Myr-scale pCO2 drawdowns took place during the warmest time intervals described for the Cretaceous11,12,13: a drawdown of ca. 310 ppm (from ca. 500 to ca. 190 ppm) during the upper Barremian-lower Aptian interval (ca. 7.5 Myr), and ca. 225 ppm (from ca. 410 ppm to ca. 190 ppm) during the upper Albian-lower Cenomanian interval (ca. 5 Myr) (Fig. 3). In fact, at the Myr-scale, trends in pCO2 are even inverse compared to temperature during significant time intervals such as the upper Barremian–lower Aptian (considering the higher resolution temperature curves based on TEX86 proxy12, 13) and the upper Albian–Santonian in which pCO2 minima are coeval with thermal maxima (Fig. 3). If pCO2 was rather low and usually decoupled from temperature at the Myr-scale, it may not have been the main long-term driver of global warming during the Cretaceous.

Figure 3
figure 3

Reconstructed pCO2 compared with available sea surface temperature (SST) records11,12,13 for the Cretaceous. Vertical grey bars represent major Oceanic Anoxic Events (OAE) occurring during the Cretaceous18. Envelopes represent ±1σ.

Time scale dependency of the relationship between pCO2 and temperature

Episodes of enhanced primary production have been associated with the kyr-scale CO2 release events that took place during major OAEs and the warmest intervals of the Cretaceous40,41,42,43,44. Barclay et al.7 showed that a 20% increase in pCO2 during OAE2 triggered extensive primary productivity leading to a global increase in rates of organic carbon burial, which eventually resulted in a pCO2 decrease down to 26% in ca. 100 kyr. This indicates that, as a consequence of a critical CO2 release event to the atmosphere, primary productivity can account for pCO2 drawdowns of at least up to the 105 yr scale. The time resolution of our record is not high enough to fully discuss pCO2 trends on the kyr scale; however, when the time uncertainty of our estimates allows clear perusal, pCO2 levels are consistently low well after the events of CO2 release to the atmosphere described in the literature18, 45 (Fig. 3). This indicates that the aforementioned primary productivity effect may exceed the kyr scales resulting in a negative balance or a low-level stasis of pCO2 up to the Myr scales. If primary production were responsible for a long-term decrease in atmospheric pCO2 during the Cretaceous, we should find a co-varying increase in atmospheric pO2 as a product of long-term enhanced photosynthesis. This is in accordance with the long-term trends obtained by the recent revision of the GEOCARBSULF model, which describes both a drawdown of pCO2 and an increase in pO2 levels from the earliest stages of the Cretaceous up to the mid Cretaceous3. Poulsen et al.4 emphasized the importance of O2 in long-term climate dynamics as under low pO2 conditions shortwave scattering by clouds and air molecules is reduced; this implies a significant increase in surface shortwave forcing, leading to increase of atmospheric water vapour and consequent enhanced greenhouse forcing which increases global surface temperatures. Thus, their simulations predict that temperature increases under low pO2 and high pCO2 conditions. However, during the Cretaceous temperature was maintained high during the late Albian–Santonian interval coinciding with relatively high pO2 3 and low pCO2 conditions3, 7 (Fig. 3). This observation reflects that both pO2 and pCO2 parameters may not be enough to explain long-term evolution of temperature.

In addition to the conclusions of previous work on atmospheric CO2 dynamics at kyr time scales7 our results indicate that the relationship between CO2 and temperature was probably scale-dependent during the Cretaceous. In the framework of a Greenhouse climate state functioning, CO2 appears to have been a main driver of global warming and primary production at timescales up to the kyr7. However, our results indicate that it was decoupled from temperature and most likely a long-term consequence of the climate-biological system at Myr scales. Direct extrapolation of these conclusions to time intervals corresponding to different climate states of the Earth is not pertinent because C cycle dynamics, and thus the evolution of pCO2, might be strongly modulated by different global scale functioning of climate factors. For instance, temperature strongly affects carbon sequestration by living organisms26, 27, and global oceanic and atmospheric circulation patterns are involved in heat transport46 and distribution of climate zones47 that strongly influence primary productivity. Different dynamics and relative effects of those factors within different climate states of the Earth6 may result in different functioning of atmospheric CO2 dynamics. However, nowadays anthropogenic CO2 release into the atmosphere also constitutes a main driver of climate change at a human time scale. At larger scales, primary productivity appears to have been an important driver of atmospheric CO2 concentration during both interglacial48 and glacial intervals49. Further work on the Myr-scale interrelationships between primary productivity, CO2 and temperature within these different climate state intervals of the Earth is needed to explore how universal are the conclusions obtained herein. Primary production constitutes a major factor involved in climate dynamics and that regulates the balance between pO2 and pCO2. Thus, productivity must be taken into consideration in the evaluation of the relative effect of O2 and CO2 on climate dynamics. Present-days anthropogenic CO2 release event can be compared in magnitude with those of the Cretaceous, reflected in a ca. −1.8‰ shift in δ13CCO2 since the industrial revolution50 against the ca. −2‰ shift associated with OAE1b16 (Fig. 1). However, its critically shorter time of development makes it a more intense disturbance event, and so its forthcoming consequences are difficult to predict. Evaluating the response time of primary production to cushion CO2 forcing will be critical to allow the prediction of climate evolution and its impact in the near future of life on Earth.

Methods

Twelve stages corresponding to eleven sites from the Cretaceous of western Europe were used for pCO2 reconstruction (see Supplementary Fig. S1). The sites represent similar terrestrial environments in which the fossil conifer genus Frenelopsis was present, and they are distributed over a narrow palaeolatitude band (30–40° N) of the same climate zone throughout the time range analyzed51 (see Supplementary Table S1). These stages and sites were selected to reduce as much as possible environmental biases on δ13Cleaf.

Frenelopsis was chosen as a pCO2 proxy because it was frequent and abundant in Cretaceous ecosystems52,53,54,55. Frenelopsis is a valuable material for performing stable carbon isotope analysis due to its resistance to decay and diagenesis. It provides accurate δ13Cleaf values, which are comparable to those of living plant leaves25, 56, 57.

Frenelopsis leaves were treated with HCl 32% following the protocol recommended by Barral et al.25. Thirty Frenelopsis leaves per locality were separately reduced to powder using an agate mortar. Three replicates of ca. 100 μg per leaf were randomly selected and weighted using a precision balance Sartorius ME36S and loaded into tin capsules. Carbon isotope analyses were performed using a continuous flow IRMS configuration, being He the carrier gas, using a EuroEA3028™–HT elemental analyzer working in combustion mode and interfaced with an IP60 isotope ratio mass spectrometer. Carbon isotope ratios were calibrated against the international standard IAEA-C4 (δ13C = −23.96‰) and casein (δ13C = −22.67‰) as an internal reference. Based on repeated analysis of standards, analytical reproducibility was better than ±0.2‰.

pCO2 values were estimated for each stage using the equation relating pCO2 to δ13Cleaf described by Schubert and Jahren19, 21 (r = 0.94):

$${{\rm{\Delta }}}^{13}{{\rm{C}}}_{{\rm{leaf}}}=[(28.26)(0.22)(p{{\rm{CO}}}_{2}+23.9)]/[28.26+(0.22)(p{{\rm{CO}}}_{2}+23.9)]$$

Δ13Cleaf values were estimated from δ13Cleaf and concomitant δ13CCO2 values by using the equation proposed by Farquhar et al.58:

$${{\rm{\Delta }}}^{13}{{\rm{C}}}_{{\rm{leaf}}}={{\rm{\delta }}}^{13}{{\rm{C}}}_{{\rm{Co2}}}-{{\rm{\delta }}}^{13}{{\rm{C}}}_{{\rm{leaf}}}/1+({{\rm{\delta }}}^{13}{{\rm{C}}}_{{\rm{co2}}}/1000)$$

δ13CCO2 values were taken from the δ13Ccarb-based δ13CCO2 estimates provided by Barral et al.16. We avoid discrimination bias due to different physiological functioning related to taxonomic specificity by using always the same taxon, and thus δ13Cleaf faithfully represents environmental variability. In order to avoid misrepresented pCO2 estimates due to subtle uncorrelations between the beds from which δ13CCO2 and δ13Cleaf were obtained, we accounted for δ13CCO2 variability in the uncertainty timespan of the Frenelopsis localities by averaging and estimating cumulated error for each stage. This was made by randomly drawing 10,000-repetition Gaussian distributions for each single case within the interval of uncertainty of a given stage, defined by the pre-determined means and standard deviations of each case, and calculating the mean and the standard deviations of the 10,000 randomly generated sets. This procedure was also implemented to include the associated errors to each δ13Cleaf, δ13CCO2 and Δ13Cleaf values into the aforementioned equations to account for all sources of error in the resolution of pCO2 estimates.