Atmospheric CO2 during the Mid-Piacenzian Warm Period and the M2 glaciation

The Piacenzian stage of the Pliocene (2.6 to 3.6 Ma) is the most recent past interval of sustained global warmth with mean global temperatures markedly higher (by ~2–3 °C) than today. Quantifying CO2 levels during the mid-Piacenzian Warm Period (mPWP) provides a means, therefore, to deepen our understanding of Earth System behaviour in a warm climate state. Here we present a new high-resolution record of atmospheric CO2 using the δ11B-pH proxy from 3.35 to 3.15 million years ago (Ma) at a temporal resolution of 1 sample per 3–6 thousand years (kyrs). Our study interval covers both the coolest marine isotope stage of the mPWP, M2 (~3.3 Ma) and the transition into its warmest phase including interglacial KM5c (centered on ~3.205 Ma) which has a similar orbital configuration to present. We find that CO2 ranged from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{394}}}_{{\boldsymbol{-}}{\bf{9}}}^{{\boldsymbol{+}}{\bf{34}}}$$\end{document}394−9+34 ppm to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{330}}}_{{\boldsymbol{-}}{\bf{21}}}^{{\boldsymbol{+}}{\bf{14}}}$$\end{document}330−21+14 ppm: with CO2 during the KM5c interglacial being \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\bf{391}}}_{{\boldsymbol{-}}{\bf{28}}}^{{\boldsymbol{+}}{\bf{30}}}$$\end{document}391−28+30 ppm (at 95% confidence). Our findings corroborate the idea that changes in atmospheric CO2 levels played a distinct role in climate variability during the mPWP. They also facilitate ongoing data-model comparisons and suggest that, at present rates of human emissions, there will be more CO2 in Earth’s atmosphere by 2025 than at any time in at least the last 3.3 million years.

www.nature.com/scientificreports www.nature.com/scientificreports/ To address this data deficiency, we developed δ 11 B-based CO 2 estimates from Ocean Drilling Program (ODP) Site 999 in the Caribbean ( Supplementary Fig. 1) at a resolution of 1 sample per 3-6 kyr for the time interval 3.15 to 3.35 Ma, encompassing the M2 glaciation and the KM5 interglacial (including KM5c). Although focusing predominantly on the mixed layer dwelling planktic species Globigerinoides ruber (45 new data points, 63 in total), we also present new measurements of Trilobus sacculifer (2 new data points, 5 total) on the same samples to provide a check on the consistency of the δ 11 B-pH calibration for G. ruber which has been recently called into question 17 .

Results
Our new high-resolution CO 2 record is shown in Fig. 2 (and Supplementary Fig. 2 Table 2) and is consistent with earlier studies [12][13][14] in showing that CO 2 was higher than the pre-industrial during the mPWP (mean=360 ppm). Our more detailed record reveals that CO 2 variations ranged from − + 331 11 13 ppm to − + 389 8 38 ppm (based on the mean and distribution of CO 2 in the <25% and >75% interquartile range during the whole studied period; ref. 18 ), with a peak-to-trough range of 47-68 ppm determined by a Welch T-test of the data within the upper and lower quartiles (at 95% confidence; p < 0.01). From a comparison with benthic δ 18 O data from ODP 999 (Supplementary Table 3) and the δ 18 O stack 19 (LR04) we observe that cold marine isotope stages (e.g. KM2; Fig. 2) are typically closely associated with low CO 2 and warm stages with high CO 2 levels. However, during the prominent M2 cold stage and through the warming out of M2, CO 2 appears to lag benthic δ 18 O by ~10 kyr (Supplementary Fig. 5). This lag is not attributable to age model uncertainty because it is present when comparing δ 11 B-derived CO 2 and benthic δ 18 O from the same samples (Fig. 2). CO 2 during the interglacial KM5c, determined using the mean of the δ 11 B of the five samples in this interglacial, is estimated at − + 371 29 32 ppm (at 95% confidence). Using a broader time window for the KM5 interglacial does not significantly influence this estimate (see Supplementary Table 1).

and Supplementary
Our estimates of borate ion, pH and CO 2 derived from the δ 11 B of T. sacculifer (without sac-like final chamber) from ODP 999 (our new data, and ref. 14 ) and ODP 926 20 overlap well with those based on G. ruber from ODP 999 in the mPWP ( Fig. 2 and Supplementary Fig. 3). Calculated CO 2 does not exhibit substantive inter-species offset with a mean difference of 10 ± 29 ppm with no consistent bias towards higher or lower pH/CO 2 .
Accuracy of δ 11 B-co 2 from G. ruber in the Plio-Pleistocene. Ref. 17 suggested that the δ 11 B-pH calibration of G. ruber may have evolved through time and that G. ruber may suffer from morphotype-differences in δ 11 B-derived pH estimates for the surface Pliocene ocean resulting in underestimates of the true pH and corresponding overestimates of true CO 2 . The principal evidence presented for this assertion was the disagreement between Pliocene CO 2 calculated using the T. sacculifer data of Bartoli et al. 12 and the G. ruber data from Martinez-Boti et al. 13 (Fig. 1). As shown here, when T. sacculifer and G. ruber are compared from the same samples and measured with the same analytical technique (in this case MC-ICPMS) there is no significant offset between the methods in terms of reconstructed borate ion δ 11 B, pH or CO 2 ( Fig. 2 and Supplementary Fig. 3).   53 ) compared to the benthic isotope stack of ref. 19 . M2 glacial and early Pleistocene strong glacials are highlighted in blue for context. Interglacial KM5 and KM5c are highlighted in yellow and orange, respectively. Note that there are no estimates for the M2 glacial and very few across the mid-Piacenzian warm period (mPWP), the low resolution of previous studies makes pin-pointing individual interglacials such as the KM5c future analogue difficult. These studies also differ in their estimates of Mid-Piacenzian CO 2 [12][13][14]20 . This finding suggests that the G. ruber δ 11 B-pH calibration has not evolved through time and therefore that the pH (and hence CO 2 ) we reconstruct here is an accurate record of the surface water carbonate system parameter. This finding indicates that disagreements between published Pliocene δ 11 B-based datasets 12,13,20,21 are likely attributable to either: (i) sampling driven aliasing due to the relatively large short-term CO 2 variability in the mPWP (e.g. Figure 2), (ii) differences in core site location between studies (with possible different local CO 2 disequilibrium), (iii) the lack of a comparison between species on exactly the same sample, or (iv) the well-documented analytical offset between MC-ICPMS and negative-ion thermal ionization mass spectrometry (NTIMS; see refs. 21,22 ). We note that, if the offset is attributable to analytical issues, the agreement between data generated by both methodologies for the younger Pleistocene time slices examined here (Fig. 2) confirms the suggestion of Sosdian et al. 20 that the T. sacculifer dataset of Bartoli et al. 12 (measured with NTIMS) requires an additional analytical correction (see ref. 16 for details), beyond that currently used for the NTIMS δ 11 B datasets of refs. 17,23 .  13 , Chalk et al. 18 ), red squares are Trilobatus sacculifer at ODP 999 (this study and Seki et al. 14 ), purple squares are T. sacculifer from ODP 668 (Honisch et al. 23 ) and blue squares are T. sacculifer from ODP 926 (Sosdian et al. 20 ). Black solid line shows ice core-derived CO 2 from ref. 58 . Left; Late Pleistocene CO 2 from boron isotopes 14,18,23,52 and ice core data. Also shown are CO 2 projections in line with RCP8.5 at current emission rates to the year 2040 (black broken line). Middle column; MPT CO 2 estimates 18,23 including disturbed ice estimates 24,25 (Note: age adjusted for scale). Right; mPWP estimates of CO 2 (this study combined with Martinez-Boti et al. 13 ), new data from T. sacculifer is shown in red squares and shows no offset from G. ruber estimates. Second panel: Time periods as above, LR04 and ODP 999 δ 18 O from C. wuellerstorfi 18,19,53 . Third panel: Iron mass accumulation rate from the Southern Atlantic ODP Site 1090 28 . Fourth panel: % Northern Component Water (NCW) estimated from δ 13 C in benthic foraminifera (grey) and ɛ Nd from fish debris (dark green) in the deep North Atlantic (core U1313 31 ). Note the lag of ocean circulation and CO 2 relative to the M2 glaciation 31 . www.nature.com/scientificreports www.nature.com/scientificreports/ mpWp co 2 cycles -variability and causes. Our new high-resolution data set clearly demonstrates that the mPWP is an interval of relatively high CO 2 , with a mean of 360 ppm whereas lower values are observed during the mid-and late-Pleistocene (262 ppm and 241 ppm 18 , respectively Fig. 2). Our dense sampling permits, for the first time, an assessment of CO 2 variability during mPWP on orbital timescales, although the length of our record still precludes reliable time series analysis. The CO 2 cycles (and derived climate CO 2 forcing) we observe in the mPWP (including the M2 glaciation) are similar, yet smaller in amplitude, than those evident in δ 11 B-CO 2 data from the Late Pleistocene (LP = 0-250 kyr; ref. 18 ) and early mid-Pleistocene Transition (eMPT = 1050-1250 kyr; Figs. 2 and 3). The variation in climate forcing between these different intervals correlates with the magnitude of δ 18 O variability, although there is an apparent increased response in δ 18 O during the late Pleistocene (δ 18 O range ~ 2‰) relative to the mPWP (δ 18 O range ~ 0.5-1‰ Fig. 3; ref. 13 ) likely reflecting the increased influence of ice-volume change on δ 18 O following the northern hemisphere glaciation 13 .
We also compare the δ 11 B-derived CO 2 data with ice-core based CO 2 for the Pleistocene intervals 18,23 ( Supplementary Fig. 4). We expect the range of CO 2 variability from the various Antarctic ice cores to be narrower than our marine-based records because of the greater precision associated with the ice core records (± 6 ppm vs. ±20-30 ppm) and the CO 2 data from blue ice from the mid-Pleistocene may not capture a full climate cycle 24,25 . Despite these caveats, and as discussed elsewhere 18,25 , there is good agreement between the ice-core and δ 11 B-derived CO 2 , providing confidence in the accuracy of the distribution (and absolute values) we determine here for the mPWP.
Three phenomena are suggested to exert a major control over glacial-interglacial CO 2 change in the Late Pleistocene record: (i) changes in stratification south of the Antarctic Polar Front influencing the ventilation of CO 2 -rich circumpolar deep water; 26 (ii) variations in the magnitude of dust-borne Fe fertilization in the sub-Antarctic that influences the strength of the biological pump in this region; 27,28 (iii) changes in the extent to which southern component deep water fills the North Atlantic increasing deep ocean carbon storage 29 . Ref. 30 proposed that polar waters in the North Pacific and the entire Southern Ocean (SO) were destratified prior to the onset of Northern Hemisphere Glaciation at 2.7 Ma. This mechanism may contribute to elevated CO 2 during the mPWP 30 , but it would also rule out water mass ventilation in the high-latitude SO as an important driver of the CO 2 cycles we observe in the mPWP. Both the accumulation and variability of dust-borne Fe at ODP Site 1090 (Fig. 2) were reduced during the mPWP to a fraction of that observed in the other time intervals examined (~ factor of four reduction during the mPWP, compared to the late Pleistocene), suggesting that dust supply did not play a major role in generating the observed CO 2 cycles in the mPWP.
Chemical stratification of the North Atlantic Ocean, due to incursions of southern component water (SCW), during the mPWP is suggested by gradients in ɛ Nd and δ 13 C from ref. 31 , albeit at reduced magnitude compared to the latest Pleistocene and eMPT. Ref. 31 showed that the SCW-signal was particularly marked after the M2 glacial maximum as shown by reduced % NCW (northern component water, Fig. 2), and also observed that ɛ Nd lagged δ 18 O by 10 kyr during the onset of the M2 glaciation, as is also evident in our new CO 2 data (Fig. 2,  and supplementary Fig. 5). In general, those intervals of the mPWP where SCW contributed significantly to the waters of the deep North Atlantic were characterized by low CO 2 (and vice versa; Fig. 2), yet with variable leads and lags (supplementary Fig. 5). This association is consistent with the suggestion 31,32 that glacial expansion of a CO 2 -charged SCW reservoir plays at least a first order role in driving orbital CO 2 cycles, even prior to the onset of northern hemisphere glaciation at 2.7 Ma.  19 and many other temperature records, such as arctic air temperature 33 and sea surface temperature 34,35 (including site ODP 999, supplementary Fig. 6) It is also often described as a failed attempt at www.nature.com/scientificreports www.nature.com/scientificreports/ Northern Hemisphere Glaciation 34,36,37 , that eventually initiated ~600 kyr later at ~2.6 Ma. Ref. 13,14 used δ 11 B-CO 2 to suggest that CO 2 dropped below 280 ppm for the first time during the intensification of Northern Hemisphere Glaciation (iNHG), consistent with the suggestion that 280 ppm CO 2 is an important threshold below which extensive glaciation of continents in the northern hemisphere develop 38 . Our new data reveal that while the lower bound of the error envelop in our smoothed CO 2 record approaches this value, atmospheric CO 2 during M2 is unlikely to have crossed this threshold. Furthermore, in our records, CO 2 variability associated with M2 lags δ 18 O by 10 kyr, which also means that minimum CO 2 is not coincident with minimum northern hemisphere orbital forcing ( Supplementary Fig. 5). Additional support for CO 2 playing a secondary role in M2 is that other periods of low CO 2 are evident in the mPWP (Fig. 2). For instance, the KM2 glaciation is clearly evident in the benthic δ 18 O data, Mg/Ca-SST at site 999 ( Fig. 2 and Supplementary Fig. 6) and our CO 2 record, but is not considered to be a major glacial interval 19 (Fig. 2). This therefore suggests that other boundary conditions, such as orbital configuration 39 (Supplementary Fig. 5) were perhaps dominant in triggering the M2 cold stage.
M2 Glaciation -CO 2 lags δ 18 o. Minor dephasing between ɛ Nd , δ 13 C and benthic δ 18 O has been observed in the late Pleistocene 40 , with carbon budget change lagging δ 18 O and preceding the change in ocean circulation, but with smaller lags than we observe during M2 (~2.5 ky vs. ~10 ky).
Similarly, while variations in atmospheric CO 2 are not entirely in phase with ice volume changes over the late Pleistocene, the leads-lags are small 41 (<3 ky) and are not readily observed when comparing δ 11 B-derived CO 2 and benthic foraminiferal δ 18 O time series (e.g. Figure 2 left panels).
This therefore either requires the operation of subtly distinct carbon-cycle dynamics during the Pliocene, and M2 in particular, or implies some other driver operated to explain the observed ~10 kyr lag of CO 2 during M2.
One possible non-carbon cycle driver for the observed lag at M2 could be a preservation bias in our data because partial dissolution of planktic foraminiferal tests drives δ 11 B towards more negative isotopic composition (lower pH, higher CO 2 ) in some species 14,42 . However, while our chosen species for pH/CO 2 reconstruction, G. ruber, is known to be relatively susceptible to partial dissolution on the seafloor 43 , its δ 11 B signal has been observed to be robust 14,44 . Furthermore, a recent study 45 of test weight and fragmentation at ODP 999 showed that tests become better preserved during M2 ( Supplementary Fig. 7), as observed during glacial periods of the late Pleistocene 46 . This is inconsistent with the high CO 2 observed during the descent into M2 maximum (especially given fragmentation and δ 18 O are in phase, Supplementary Fig. 7), ruling out an effect of partial dissolution on our CO 2 reconstruction and observed lag.
An alternative explanation for the delayed pH change observed at ODP 999 during M2 lies in local water mass changes during this glacial episode. Although the Central American Seaway was closed to deep water by this time 36 , it is possible that limited exchange of surface water was still occurring around M2 34 and the early Pleistocene 47 . Temperature data from Mg/Ca in T. sacculifer 34 and G. ruber at site 999 (this study, Supplementary  Fig. 6, data from the East Equatorial Pacific site 1241 48 is also shown for comparison) show a cooling prior to the M2 maximum. A connection between the East Equatorial Pacific and the Caribbean, or local intensification in upwelling could have brought cold, nutrient-and carbon-rich waters (low Mg/Ca, low pH) to Site 999, potentially explaining the apparent high CO 2 at the inception of M2. However, because no noted decline in the abundance of G. ruber occurred, given this is a species known to favour oligotrophic conditions 49 , it seems unlikely this was accompanied by a significant influx of such water. Also, a mechanism for increased influx of EEP water during sea level regression is lacking. Importantly, during KM5c (3205ky), the CAS likely remained closed impeding the influx of Pacific waters to site 999, as shown by elevated temperature (relatively to the M2 interval) in the Caribbean ( Supplementary Fig. 6), hence local changes in hydrography are unlikely to have unduly impacted our CO 2 estimates for this central interval.
The cause of the apparent lag of CO 2 compared to δ 18 O during the inception of M2 therefore remains uncertain, but because this lag is also present in North Atlantic ɛ Nd records 31 we favour a carbon cycle-based interpretation. A Southern Ocean-driven explanation for the CO 2 lag during M2 is perhaps indicated by the observation that the tail of the CO 2 decline during M2 is out of phase with decline in northern hemisphere insolation, but apparently in-phase with insolation decline at 65°S (i.e. offset by one half precession cycle, Supplementary Fig. 5).
Past to future -when will we exceed mPWP CO 2 levels?. Atmospheric CO 2 has been increasing since the industrial revolution from a background of ~280 ppm, reaching 411 ppm in 2019. The mPWP is often cited as being the last time CO 2 levels were as high as today 13 , although we note methane and other greenhouse gases remain poorly constrained and may contribute to extra radiative forcing 50 .
Our refined view of CO 2 during this interval however reveals that in terms of the mean (360 ppm), mPWP values were exceeded in the mid-1990s. Our upper quartile range ( − + 389 8 38 ppm) suggests that CO 2 during the mPWP is very likely to have been ≤ 427 ppm (using the distinctions of the IPCC). Atmospheric CO 2 rose by 2.5 ppm from 2017 to 2018, if this rate is sustained, our new data indicate that CO 2 will surpass even the highest mPWP values within the next 4 to 5 years (i.e. by 2024-2025).

Conclusions
Our new δ 11 B-data for the mPWP provide a tightly constrained and robust estimate of CO 2 during KM5c of − + 371 29 32 ppm (at 95% confidence), and documents CO 2 variability during the mPWP from − + 331 11 13 ppm to − + 389 8 38 ppm. This extended view suggests that changes in ocean carbon storage may play an important role in natural variability in CO 2 on orbital timescales in warmer than present climate states. By better constraining the upper levels of CO 2 during the mPWP we conclude that, given current rates of CO 2 increase, we will very likely surpass mPWP values within 4 to 5 years, meaning that, by 2024/2025 levels of atmospheric CO 2 will be higher than any point in the last 3.3 million years.
www.nature.com/scientificreports www.nature.com/scientificreports/ Methods ODP Site 999 is located in the Caribbean Sea and has been reliably used to reconstruct past atmospheric CO 2 in a number of previous studies 18,51,52 . Air-sea disequilibrium for CO 2 in the surface waters in the region of ODP Site 999 are close to 0 (+20 ppm) today and ref. 13 suggest this remained the case for at least the last 3.3 million years. An age model for the interval 3.15 to 3.35 Ma studied here was constructed by analyzing the epibenthic foraminifera Cibicidoides wuellerstorfi for δ 18 O in each sample, combining these data with similar data from ref. 53 , and tuning the combined record to the LR04 benthic δ 18 O stack 19 (Fig. 1). From each sample, ~120 individual tests of the mixed layer dwelling species Globigerinoides ruber white sensu stricto (300-355 μm) were hand separated, clay removed, oxidatively cleaned and analyzed for boron isotopic (δ 11 B) and trace element composition (e.g. Mg/Ca, Al/Ca) at the University of Southampton using well-established procedures 21,52 . These data were converted to pH using the G. ruber δ 11 B-pH calibration of ref. 51 and a modern δ 11 B of seawater 54 (39.6‰). Sea surface temperatures (SST) were derived from each samples Mg/Ca ratio using the calibrations of ref. 55 corrected for the changing Mg/Ca for seawater following ref. 20 based on ref. 56 .
In order to calculate atmospheric CO 2 from the reconstructed surface water pH we use dissolved inorganic carbon (DIC) from the reconstructions of ref. 20 . Uncertainty in this term is also propagated into CO 2 uncertainty using a Monte Carlo approach (n = 10,000), but with a uniform distribution encompassing the range of the reconstructed DIC for our study interval (1765 to 1840 μmol/kg). CO 2 was then determined by the maximum probability of all 10,000 realisations.
To obtain CO 2 during the KM5c interglacial the mean and associated uncertainty of the reconstructed borate ion δ 11 B and SST of the data lying within +/− 7 kyr of the peak of KM5c (n = 5) was determined. In our age model we defined that the peak of KM5c was at 3212 kyr, although expanding the window to ±15 kyr, or uncertainty in the identification of the peak of KM5c to within 10 kyr, did not have a significant effect on our calculated mean (Supplementary Table 1).
This average was then used to calculate pH and CO 2 , but only considering the uncertainty in δ 11 B sw and DIC detailed above. This approach was chosen because it allows the random uncertainties arising from SST and δ 11 B measurement error to be reduced through replication, whilst still realistically accounting for the systematic uncertainties in δ 11 B sw and DIC.
In order to place our new data into a Plio-Pleistocene context, we recalculated the δ 11 B data of ref. 13 from 3.0 to 3.3 Ma using the same methodology outlined here and combined it with our new data. The difference in the method is the choice of the second carbonate parameter where constant modern alkalinity is used in ref. 13 whereas we use DIC 20 in this study. These new and published data are then combined into a single record with an average temporal resolution of 1 sample every 4 kyr. This was then interpolated to 1 kyr resolution and a 6-point running mean was used to reduce the influence of analytical uncertainty on our reconstructed CO 2 record. Uncertainty in the smoothed record was determined using the output of the original Monte Carlo simulation. CO 2 forcing was calculated relatively to preindustrial CO 2 (278 ppm) and is defined as: www.nature.com/scientificreports www.nature.com/scientificreports/ www.nature.com/scientificreports www.nature.com/scientificreports/