Oligocene-early Miocene paradox of pCO₂ inferred from alkenone carbon isotopic fractionation and sea surface temperature trends

Atmospheric carbon dioxide decline is hypothesized to drive the progressive cooling over the Cenozoic. However, the decline in the phytoplankton carbon isotopic fractionation (ε_p) from the early Oligocene to Miocene time interval, interpreted as a long term CO₂ decline, differs from the apparent long term stability in climate indicators like benthic oxygen isotopes. Here, we produce two new long-term records of ε_p over the Oligocene to early Miocene time interval from widely separated locations at IODP Site 1406 and ODP 1168 and increase the resolution of determinations at the equatorial Atlantic ODP 925. These new results confirm a global trend of ε_p decline occurring during this interval. Rapid 3 ‰ declines are found from 27 to 24.5 million years ago (Ma) and 24 to 22.5 Ma, and minimum ε_p is attained at 19 Ma. Between 29.7 and 28.7 Ma at IODP 1406, a 20–30 ky sampling resolution reveals orbital scale 100 kyr cyclicity in ε_p. Making use of alkenone-based sea surface temperature (SST) estimates and benthic δ¹⁸O estimates extracted from the same samples, we perform a direct comparison with ε_p to evaluate the relationship with climate. We observe that across the long Oligocene to early Miocene interval, ε_p is positively correlated to SST only at the Southern Ocean Site 1168, but not with SST at the North Atlantic Site 1406. Accounting for the temperature-driven growth rate or cell size effects on ε_p does not lead to stronger correlations between ε_p and benthic δ¹⁸O nor stronger correlations between ε_p and SST at Site 1406. Moreover, at orbital timescale, the relationship between ε_p and benthic δ¹⁸O, albeit weak, implies greater ice volume or colder deep ocean at higher CO₂. Despite the persistence of climate paradox, the reproducible trends in three widely separated sites, which experienced contrasting temperature evolution and likely experienced different variations in nutrient availability, suggest that a common CO₂ forcing is likely the dominant control on the long term trends in ε_p. Changing ocean heat transport to the North Atlantic may contribute to the observed decoupling of long term ε_p and SST in this location.

1 Introduction

1.1 Oligocene to Miocene long term CO₂ and climate trends

Geological records provide key context to current assessments of the consequences of rising atmospheric CO₂ on ice sheet stability and oceanic temperatures (Foster et al., 2017; Golledge, 2020; Zachos et al., 2008). The Oligocene to Miocene time interval has been proposed to represent a nonlinear transition between the “greenhouse” and “icehouse” stages of Earth history (Miller et al., 1991; Zachos et al., 2001) useful to evaluate the Earth system climate sensitivity to the hypothesized progressive CO₂ drawdown of the Cenozoic (Deconto et al., 2008; Zhang et al., 2013). However, the long term decline in CO₂ estimated by existing proxy records contrasts with the rather stable climatic state with multimillion year warming (e.g. Late Oligocene Warming) and cooling (e.g. Mi-1 glaciation) trends interpreted from deep ocean (Cramer et al., 2011; Lear et al., 2000), and surface ocean records (Guitián et al., 2019; Liu et al., 2009; O'Brien et al., 2020) and variable Antarctic Ice sheet volume and sea level (Lear et al., 2004; Liebrand et al., 2017; Miller et al., 2020).

The long term pCO₂ trends from the Oligocene to early Miocene are derived from the sensitivity of marine algae to pCO₂ based on the carbon isotopic fractionation in organic matter during photosynthesis (ε_p) of marine phytoplankton (Henderiks and Pagani, 2008; Pagani et al., 1999, 2000, 2005; Super et al., 2018; Zhang et al., 2013) while published δ¹¹B based CO₂ estimates cover the latest Oligocene into early Miocene (younger than 24 Ma). The algae isotopic fractionation can be reconstructed in the past from sediments by the analysis of δ¹³C of the organic lipids and reconstruction of the δ¹³C of the DIC in the seawater from which biomass was produced. Fractionation (ε_p) is predicted to be higher when CO₂ availability is high relative to cellular carbon demand. A decrease in atmospheric CO₂ and consequently in CO₂ of the surface ocean should therefore lead to a decrease in ε_p at any given site. However, in addition to CO₂, the ε_p in phytoplankton is affected by physiological factors such as the rate of carbon fixation, which may vary over time in a given location due to variations in temperature or the supply of light (Rau et al., 1996; Stoll et al., 2019).

One approach to evaluate the relative contribution of physiological factors vs CO₂ is to produce ε_p records from sites of widely contrasting oceanographic setting, where the CO₂ signal may be expected to be common to both locations but the environmental factors affecting the fractionation such as nutrient availability might not be expected to change in unison. In this study, we produce a new long-term record of ε_p over the Oligocene to Miocene time interval at two new, widely separated locations: IODP Site1406 in the subtropical North Atlantic off the Newfoundland coast, and ODP 1168 in the Southern Ocean off Tasmania. We also increase the resolution of determinations at the equatorial Atlantic ODP 925. The existing ε_p-based CO₂ estimations for the Oligocene are derived from ∼ 1 million year resolution measurements from two sites (Site 925 and 516) on the South American margin of the equatorial and South Atlantic; in the early Miocene an additional North Atlantic record (Site 608) provides data (CenCO2PIP Consortium, 2023).

Our new < 1 million year resolution ε_p records from these two mid- latitude locations allow us to directly compare ε_p with estimates of SST derived from alkenones extracted from the same samples, since unlike very warm tropical locations, the $U_{\mathrm{37}}^{k^{\prime }}$ index still retains sensitivity to temperature in the mid-latitudes during the Oligocene and early Miocene. Additionally, we compare ε_p with benthic δ¹⁸O available from the same sediments. Variations in benthic δ¹⁸O are controlled by changes in both deep-water temperature and deep ocean δ¹⁸O_sw which reflects ice volume. We further measure ε_p and benthic δ¹⁸O at approximately 20–30 ky resolution over a series of eccentricity cycles in the early Oligocene at IODP 1406. The dataset allows a robust evaluation of the relationship between ε_p and climate for this time interval. We further discuss the significance of the observed ε_p record with the implications for the phytoplankton sensitivity over multimillion year timescales over the Cenozoic.

1.2 An overview of alkenone ε_ppCO₂ proxy

The carbon isotopic fractionation in phytoplankton during photosynthesis is affected not only by the aqueous carbon dioxide (CO_2[aq]) but also by physiological factors related to the cellular uptake of carbon. Physiological factors were initially modelled from the assumption of diffusive carbon acquisition in phytoplankton cells (Rau et al., 1996), where higher ε_p could be induced by higher CO_2(aq), lower instantaneous growth rates, or a higher cellular surface area to volume ratio. Both cellular permeability and the carbon isotopic fractionation by the Rubisco enzyme have been assumed to be constant, with Rubisco fractionation typically estimated between 25 ‰ and 29 ‰ for alkenone producers (Pagani, 2014). Traditional attempts to reconstruct pCO₂ from ε_p have simplified this original diffusive model by relating ε_p and CO₂ with a single factor b defined to include all physiological parameters affecting the fractionation, and ε_f representing the fractionation of the Rubisco enzyme (Jasper et al., 1994).

$$\begin{array}{}\text{(1)}& {\mathit{\epsilon }}_{\mathrm{p}}={\mathit{\epsilon }}_{\mathrm{f}}-{\displaystyle \frac{b}{{{\mathrm{CO}}_{\mathrm{2}}}_{\left[\mathrm{aq}\right]}}}\end{array}$$

The b-value has been estimated from modern photic zone and culture samples, for which CO_2(aq) is independently known. For sedimentary alkenones, previous pCO₂ calculates have either (1) assumed the modern b-value for that oceanographic setting remained constant in the past (e.g. Zhang et al., 2013), (2) applied modern relationships between b-value and phosphate and a simulated paleo-surface ocean phosphate concentration at the site (Pagani et al., 2011) (3) estimated the difference between the modern b-value at the site and the paleo-setting b-value from productivity proxies (Bolton et al., 2016) or (4) applied variation in the b-value at the site based on proxies for coccolithophore size (Henderiks and Pagani, 2007). Despite the appeal of this approach, a recent re-evaluation of cultures and field observations suggest the b term is not well predicted by growth rate, light or cell size alone in a diffusive model but that additional effects occur from carbon concentration mechanisms (CCM) on carbon uptake at lower CO₂ concentrations, which cause a deviation in the CO₂ dependence from the theoretical hyperbolic relationship (Stoll et al., 2019; Hernández-Almeida et al., 2020). A further challenge to the physical diffusive model is that the Rubisco fractionation in coccolithophores has been measured in-vitro as 11 ‰ rather than 25 ‰ (Boller et al., 2011), suggesting that fractionations larger than 11 ‰ may reflect the operation of additional enzymatic fractionations (Wilkes et al., 2018). The lower Rubisco fractionation implies a lower sensitivity of ε_p to CO₂ (e.g. as explored in González-Lanchas et al., 2021).

A meta-analysis of experimental culture data (Stoll et al., 2019) suggests that ε_p features a logarithmic dependence on CO₂, rather than the hyperbolic dependence implied by Rau et al. (1996). This analysis does not resolve the mechanisms for the observed relationship between ε_p and CO₂, but over the range of CO_2(aq) from 5 to 30 µM, it provides an empirical relationship for interpreting the magnitude of CO_2(aq) change implied by a given ε_p change. The culture dataset illustrates more broadly how ε_p is the sum of its dependencies on ln(CO₂), ln(light), and growth rate μ_i and cell radius:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\mathit{\epsilon }}_{\mathrm{p}}& =\mathrm{2.66}\ln \left({\mathrm{CO}}_{\mathrm{2}}\right)+\mathrm{2.33}\ln \left(\mathrm{light}\right)-\mathrm{6.98}{\mathit{\mu }}_i\\ & -\mathrm{1.28}\mathrm{radius}+\mathrm{6.26}\end{array}\end{array}$$

where CO_2(aq) is in µM, light is in µE, growth rate μ_i is d⁻¹ and radius is in microns (see Stoll et al., 2019 for confidence intervals on the regression).

From this empirical culture calibration, two challenges remain for the estimation of past CO₂ from ε_p measurements derived from sedimentary alkenones. First, its use would require an estimation of the cell radius, light during the season and depth of alkenone production, and the growth rate. While cell size can be estimated from coccolith length (Henderiks and Pagani, 2007), determining the absolute light and growth rate is rarely possible. Since the equation is a linear sum of these influences, these non-CO₂ variables may be integrated into the intercept (e.g. as in González-Lanchas et al., 2021), where the intercept (I) would decrease with higher growth rates and larger cell sizes and increase with higher light.

$$\begin{array}{}\text{(3)}& {\mathit{\epsilon }}_{\mathrm{p}}=\mathrm{2.66}\ln \left({\mathrm{CO}}_{\mathrm{2}}\right)+I\end{array}$$

Yet, as with Eq. (2), there remains the challenge of determining which value should be used for the intercept for past conditions, and whether a constant or variable I is more appropriate for a given site since there are limited proxies for algal growth rate. Recent culture studies document a 0.5 ‰ decrease in ε_p per 1 °C warming (Torres Romero et al., 2024), and show that this magnitude is identical to the product of ε_p dependence on growth rate (Stoll et al., 2019) and the modeled temperature dependence of coccolithophore growth rates (Krumhardt et al., 2017) derived from diverse culture and field studies (Fielding, 2013; Sherman et al., 2016; Behrenfeld et al., 2005). This suggests that growth rate, and I, may vary over time at a given location if temperature is variable. Therefore, records of SST from alkenone unsaturation or other proxies provide the opportunity to deconvolve the effects of temperature-driven growth rate variations on ε_p even when the absolute growth rate is not known.

In this study, given the potential for oceanographic conditions at the studied sites to differ significantly from those in the modern ocean at these locations, and the concomitant high uncertainties in estimating an appropriate b-value for the traditional approach or I for Eq. (3), we do not calculate the past absolute CO₂ concentration from our ε_p measurements. Instead, we account for the potential influence of temperature-driven growth rate changes on our ε_p records using alkenone temperature estimates derived from the same samples. Similarly, we evaluate the potential impact of cell size variations on the ε_p changes. Then, we employ the sensitivity of ε_p to CO₂ in Eq. (3) to estimate possible relative changes in CO₂ in the case where the other nutrient-stimulated growth rate or light influences on ε_p were constant during the studied interval at each site evaluating evidence for this assumption.

2 Sites and sediments

We have selected two widely separated paleo locations for this study, from the mid latitude North Atlantic Ocean and the high to mid latitudes of the Southern Ocean, from IODP 1406 (40°21.0^′ N, 51°39.0^′ W; 3814 m b.s.l.), ODP 1168 (43°36.5^′ S, 144°24.7^′ E, and 2463 m b.s.l.), respectively (Fig. 1). A total of 43 and 34 sediment sampling spreading from 30 to 17 Ma at each site were selected. Additionally, 61 samples (at approximately 15 ky sampling interval for bulk carbonate isotopes) were obtained from IODP 1406 within the 30–29 Ma time window, of which 29 were processed for benthic foraminiferal isotopes and 22 yielded biomarkers sufficient for analysis. We also measured an additional six sample set from equatorial Site ODP 925 (4°12.25^′ N, 43°29.33^′ W, 3042 m b.s.l.) in order to increase the million-year scale resolution of the previous longest Oligocene record in alkenone carbon fractionation (Zhang et al., 2013).The age model for Sites 1406 and 1168 has been updated using new ${}^{\mathrm{87}}\mathrm{Sr}{/}^{\mathrm{86}}$Sr isotope stratigraphy (Stoll et al., 2024) and the age-modelling software Bacon (Blaauw and Christen, 2011). The new age model for IODP Site 1406 is comparable to previously published chronologies (e.g. as in Guitián et al., 2019; Van Peer et al., 2017) but clarifies the duration of the upper Miocene hiatus between 33.3 and 34.7 m core depths (CCSF-A) as extending from 18.5 to 21 Ma. The age model of ODP Site 1168 was revised with the Sr isotope stratigraphy in the interval from 562 to 278 m b.s.f. The Site 1168 chronology is significantly shifted for most of the Oligocene to early Miocene compared with previous biostratigraphically-based age models and shipboard magnetostratigraphy (Pfuhl and McCave, 2003). Strontium isotope stratigraphy identifies a condensed interval from 22.5 to 21.6 Ma, but suggests sustained sedimentation thereafter through 16 Ma. The age model is most uncertain between 27 and 25 Ma where the Sr isotopic curve has a low rate of change. For the two ODP 1168 samples deeper than the Sr isotope measurements (562 m b.s.f.), and those from Site ODP 925 we use previous age constraints as published previously by Guitián et al. (2020).

Figure 1 Location of the study sites. Reconstructed map of continental distribution over the 30 Ma (grey)–17 Ma (black) time interval. Modified after the plate tectonic reconstruction service from the Ocean Drilling Stratigraphic Network (https://www.odsn.de/, last access: 4 November 2025) using the data from Hay et al. (1999).

The paleolatitude reconstruction for the Oligocene to early Miocene barely changes the position of Atlantic sites, in contrast, ODP 1168 moved from 55 to 48° S between 30 and 15 Ma (Torsvik et al., 2012; van Hinsbergen et al., 2015). Paleoecological reconstructions from dinoflagellates confirm that the waters above Site 1168 were continually influenced by the Leeuwin Current and located well equatorward of the Subtropical Front (Hou et al., 2023b) during the Miocene. Paleodepth estimates for coastal site ODP 1168 suggest a gradual deepening from the Eocene onwards (Exon et al., 2001).

3 Methods

3.1 Alkenone purification, quantification and δ¹³C analysis

Biomarkers from sediments of IODP 1406 and ODP 1168 were extracted from 30 g of freeze-dried sediment using an Accelerated Solvent Extractor 350 with CH₂Cl${}_{\mathrm{2}}/\mathrm{MeOH}$ (9:1 $\mathrm{v}/\mathrm{v}$) solvent for four static cycles at 100 °C and further silica gel column chromatography protocols for purification of the ketone fraction containing the alkenones (see Guitián et al., 2019 for details).

Alkenone ratios were obtained with a Thermo Scientific Trace 1310 Gas Chromatograph (GC)-FID. Originally, for IODP 1406 and 1168 samples the GC-FID was equipped with a non-polar (60 m × 0.25 mm × 0.25 µm) capillary column (ZB-1ms, Zebron^™) at ETH Zurich and at Lamont-Doherty Earth Observatory from Columbia University (Guitián et al., 2019; Guitián and Stoll, 2021). However, ODP Site 1168 samples older than 22.4 Ma featured more complex chromatograms and a high diversity of compounds. To reduce the effects of coelution, samples were additionally analyzed on a 105 m column RTX-200ms at ETH Zurich, which improved separation of long chain ketones (Rama-Corredor et al., 2018). The following temperature program was used: 1 min at 50°C, temperature gradient of 40 °C min⁻¹ to 200°C and 5 °C min⁻¹ to 300°C, hold for 45 min, and increased to 320°C at 10 °C min⁻¹ and hold for 8 min. Carrier gas was Helium at a flow rate of 1.5 mL min⁻¹. In-house standards and replicates injected at every sequence ensured instrument precision. A subset of IODP Site 1406 and samples younger than 22.4 Ma from Site 1168 were re-measured with the RTX-200 to ensure replicability (Table S1 in the Supplement). Method used for each organic analysis is described in the supplementary material dataset.

Sea surface temperature was calculated from $U_{\mathrm{37}}^{k^{\prime }}$ ratio using the Bayspline calibration (Tierney and Tingley, 2018) for all samples in IODP 1406 and for the young set (< 23.1 Ma) of samples in ODP 1168. Because for those, we find $U_{\mathrm{37}}^{k^{\prime }}$ ratios within the analytical uncertainty using both columns, we report the original ZB-1 results for all Site 1406 samples and Site 1168 samples younger than 22.4 Ma/408.22 m depth. The RTX-200 column provided substantially improved resolution of C38 peaks, allowing quantification of C_38:2 and C_38:3 ME peaks. For samples between the ages of 23.1 and 29.1 Ma in ODP 1168 the RTX-200 column still did not sufficiently resolve coelutions on the C_37:3 peaks. Therefore, for this interval we provide temperatures estimated from the $U_{\mathrm{38}\mathrm{ME}}^{k^{\prime }}$ ratio applying the Novak et al. (2022) core top calibration (Table S1). For the ODP 925 equatorial site samples, the C_37:3 methyl ketone is under the detection limit, therefore we further purified and analyzed the extracted organics as in Guitián et al. (2019) to get the temperatures from the TEX₈₆ ratio using the BAYSPAR calibration by Tierney and Tingley (2015).

Compound-specific δ¹³C measurements were performed on a Thermo Scientific Trace 1310 Gas Chromatograph coupled to a Thermo Scientific GC Isolink II, a Conflo IV, and a Delta V Plus Mass Spectrometer at ETH Zurich. Oxygen was flushed through the combustion reactor for one hour at the beginning of each sequence and seed oxidized for one minute before each injection. Alkenones from ODP Site 1168 and IODP 1406 were analyzed on a GC equipped with a non-polar capillary column (60 m × 0.25 mm × 0.25 µm) (ZB-1ms, Zebron^™) and 5 m guard column. Helium was used as carrier gas flow with 2 mL min⁻¹. GC oven was set to 90 °C, ramped to 250 at 25 °C min⁻¹, to 313 °C at 1 °C min⁻¹ and finally to 320 °C at 10 °C min⁻¹. The GC oven was then maintained isothermally for 20 min. A subset of IODP 1406 samples were measured additionally at the Lamont-Doherty Earth Observatory on equivalent instrumentation but with some modifications for improved sensitivity (Baczynski et al., 2018) and following similar GC procedures, with similar results.

From ODP Site 1168, a subset of samples older than 24.5 Ma featuring more complex chromatograms were rerun on a GC-irMS equipped with a RTX-200 ms column. GC oven was set at 50 °C ramped to 275 °C at 40 °C min⁻¹, to 295 °C at 0.5 °C, hold for 22 min, and finally ramped to 320 °C at 10 °C min⁻¹ and hold for 5 min. Flow rate was 1.5 mL min⁻¹. Comparison of a subset of samples from Site 1406 and ODP 1168 younger than 22.4 Ma showed that δ¹³C C37:2 were similar on both ZB-1 and RTX-200 columns. We consequently report here δ¹³C C37:2 from the ZB-1 runs, with the exception of the samples from Site 1168 older than 22.4 Ma. All values are reported here in parts per mil (‰) relative to VPDB (Vienna Pee Dee Belemnite). Sample replicates, in-house alkenone standard (provided by G. O'Neil, Western Washington University, and C. M. Reddy, Woods Hole Oceanographic Institution), and known isotopic mixtures A5 and B4 (supplied by A. Schimmelmann, Univ. of Indiana) were simultaneously measured to determine the analytical accuracy of the measurement and an uncertainty of 0.5 ‰.

3.2 Estimation of aqueous carbon dioxide δ¹³C

Isotopic composition of CO_2[aq] is estimated from the temperature dependent fractionation between DIC and aqueous CO₂ during alkenone production of Rau et al. (1996) based on Mook et al. (1974) and Freeman and Hayes (1992):

$$\begin{array}{}\text{(4)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]\mathrm{aq}}={\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{DIC}}+\mathrm{23.644}-\left({\displaystyle \frac{\mathrm{9701.5}}{T}}\right)\end{array}$$

We calculate the δ¹³C DIC from the δ¹³C measured on the bulk carbonate, which is dominated by calcareous nannofossils, for which previous studies show Reticulofenestra to be the most abundant genera (Guitián et al., 2020). Because there is no divergence of vital effects between small and large coccoliths in the late Oligocene to early Miocene (Bolton and Stoll, 2013), we propose that the offset between coccolith δ¹³C and DIC is likely to remain constant. We subtract 0.5 ‰ from the δ¹³C bulk to calculate δ¹³C DIC, based on average alkenone-producing coccoliths cultured at DIC < 4 mM compiled in Stoll et al. (2019). Support for estimating photosynthetic fractionation from coccolith δ¹³C is provided by recent culture studies of G. oceanica (Torres Romero et al., 2024). In previous studies, the δ¹³C DIC has also been estimated from the δ¹³C of calcium carbonate of benthic foraminifera with the assumption of a constant and known offset between the δ¹³C DIC of the deep and surface ocean. Although the foraminifera content in Site 1406 and 925 is very low, sediments feature sufficient well preserved benthic foraminifera, mainly epifaunal Cibicidoides spp. in the size range larger than 200 µm. At ODP Site 1168 benthic foraminifera were scarce for picking for isotopes in many intervals and the progressive evolution of water depth at the site may change the δ¹³C offset between the benthic environment and the surface ocean over time (Exon et al., 2001). For an additional sensitivity test to evaluate the significance of the method of DIC estimation and facilitate comparison to other published ε_p records calculated from benthic δ¹³C, we also estimate surface ocean DIC by adding a constant offset of +2 ‰ to the δ¹³C benthic measurements, following previous Miocene and Oligocene studies (Guitián et al., 2019; Pagani et al., 2011; Zhang et al., 2013).

Bulk carbonate and benthic foraminifera were measured using analytical techniques described in Guitián et al. (2019) with the guidelines from Breitenbach and Bernasconi (2011) for small carbonate samples on a GAS BENCH II Delta V Plus irMS from Thermo Scientific with international (NBS-19 & 18) and in-house carbonate as standards achieving a precision of 0.07 ‰.

3.3 Calculation of ε_p 37.2  and CO₂ doublings

Carbon isotopic fractionation (ε_p), describes the fractionation occurring during photosynthesis when carbon is fixed into algal cellular biomass (δ¹³C_org) from the ambient aqueous CO₂ (δ¹³C${}_{{\mathrm{CO}}_{\mathrm{2}}\mathrm{aq}}$) (Freeman and Hayes, 1992):

$$\begin{array}{}\text{(5)}& {\mathit{\epsilon }}_{\mathrm{p}}=\left({\displaystyle \frac{\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\left[{\mathrm{CO}}_{\mathrm{2}}\right]\mathrm{aq}}+\mathrm{1000}\right)}{({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{org}}+\mathrm{1000})}}-\mathrm{1}\right)\cdot \mathrm{1000}\end{array}$$

Organic δ¹³C is obtained from the δ¹³C analysis of haptophyte specific alkenone di-unsaturated C_37.2. Culture experiments showed that the lipid organic matter is depleted in ¹³C relative to the whole cell isotopic composition by 4.2 ‰, a correction that needs to be applied (Wilkes et al., 2018; Popp et al., 1998):

$$\begin{array}{}\text{(6)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{org}}=\left[\right({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{37}:\mathrm{2}}+\mathrm{1000})\cdot ((\mathrm{4.2}/\mathrm{1000})+\mathrm{1})-\mathrm{1000}]\end{array}$$

Uncertainties were propagated by a full Monte Carlo (n= 10 000) simulation following Tanner et al. (2020).

To compare our new records with previous data spanning the same time interval, we discuss published ε_p datasets recently compiled by the paleo CO₂ community (CenCO2PIP Consortium, 2023), from DSDP 516 in the South Atlantic (Pagani et al., 2000, 2005, 2011), DSDP 608 in the North Atlantic (Super et al., 2018), and the equatorial site from ODP 925 (Zhang et al., 2013). For these, we ensure that ε_p for the published records is calculated from biomarker-based paleothermometers. The most recent publications from DSDP 608 and 925 used glycerol dialkyl glycerol tetraethers (GDGT) derived estimations from TEX₈₆. To better compare our results with DSDP 516, where originally temperatures were derived from δ¹⁸O of planktic foraminifera for the Miocene section and GDGTs for part of the Oligocene, we have updated the ε_p calculations using a running averaging of the recent higher resolution GDGT temperature reconstructions from Auderset et al. (2022) at the same site.

For our data of paired ε_p and alkenone SST, we calculate the shift in ε_p which is expected from temperature-stimulated growth rates. Using the relationship of −0.48 ‰ (95 % CI = −0.37 to −0.95 ‰) per 1 °C SST (Torres Romero et al., 2024), we adjusted each samples ε_p absolute value by using the difference between the SST estimated for that sample and the average SST during the studied interval at that site. We complete a similar exercise for cell radius, calculating the deviation in ε_p only relative to the median cell size, for each point using the culture dependence of ε_p on cell radius shown in Eq. (2) (Stoll et al., 2019). Biogenic silica (bioSi) was determined on 20 samples from ODP 1168 following methods described previously in Guitián et al. (2020).

From the ε_p time series we estimate the change in CO₂ relative to the maximum values at 29 Ma, using the adjustment in ε_p for temperature sensitive growth rate described in the previous paragraph, and Eq. (3) as applied in González-Lanchas et al. (2021), where I reflects the size and light influences on ε_p and is assumed constant across all time intervals, and the ε_p dependence on ln [CO₂[aq]] of 2.66 is the 50th percentile estimate of the modern cultures. We then estimate the doubling/halving of CO₂ relative to the CO₂ at the reference age (R) applying the solubility for the measured temperature (Zeebe and Wolf-Gladrow, 2001) which can be reduced to:

$$\begin{array}{}\text{(7)}& \mathrm{doubling}{\mathrm{CO}}_{\mathrm{2}}={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{p}}\left(t\right)-{\mathit{\epsilon }}_{\mathrm{p}}\left(R\right)}{\mathrm{2.66}\ln \mathrm{2}}}+{\log }_{\mathrm{2}}\left({\displaystyle \frac{\mathrm{sol}\left(R\right)}{\mathrm{sol}\left(t\right)}}\right)\end{array}$$

4 Results and Discussion

4.1 Trends in ε_p in the Oligocene to early Miocene

4.1.1 New ε_p records from sites 1406 and 1168

In both long-term records from Site 1406 and Site 1168, δ¹³C of C_37.2 alkenones range from very low values near −30 ‰ in the early Oligocene (30–28 Ma) increasing to −24 ‰ by 20–18 Ma (Fig. 2). The new calculated ε_p decrease from the Oligocene to the early Miocene at the highest resolution North Atlantic Site 1406, features abrupt 3 ‰ declines from 27 to 24.5 Ma. At ODP 1168 the lack of ε_p measurements prior to 25 Ma hinders identification of a clear transition. Newly obtained ODP 925 ε_p determinations within the interval 25 to 19 Ma are in agreement with previous determinations at this site (Zhang et al., 2013) showing also the broad interval of ε_p decrease from 27 to 24.5 Ma. The trends in ε_p calculated from benthic δ¹³C are similar to those calculated from the coccolith-dominated bulk δ¹³C. In the high resolution section from 29.6 to 28.8 Ma in Site 1406 there is no long-term trend, but orbital scale ε_p variations exceed 1.5 ‰ in amplitude (Figs. 3 and S1). Over several ∼ 100 ky orbital cycles, variations of 0.75 ‰ benthic δ¹⁸O and bulk δ¹⁸O are observed, consistent with previous findings of high 100 ky power in benthic δ¹⁸O in other sites during this time period (Liebrand et al., 2017).

Figure 2 Analytical results of this study. (a) Carbonate stable isotopes for benthic foraminifera (red crosses), data from this work for each site; black line shows results from Westerhold et al. (2020) and bulk sediment carbonate (green lines and cricles). (b) Alkenone C_37.2δ¹³C. (c) Calculated alkenone carbon fractionation, solid line DIC δ¹³C is derived from bulk carbonate, transparent line from picked benthic foraminifera at the same samples. (d) Temperature estimates. IODP 1406, including SST and benthic δ¹⁸O dataset from Guitián et al. (2019); ODP 1168, dark green δ¹³C bulk carbonate shows the 4-point moving average used to calculate ε_p at the site; ODP 925, filled circles are new measurements for this study, white symbols are published data (Zhang et al., 2013), being ε_p recalculated following method described in text, SST is derived from GDGT at this site. Alkenone SST-derived show 1σ error bars and GDGTs 925-TEX₈₆ refer to 2σ.

Figure 3 High resolution sampling from IODP 1406, showing δ¹⁸O of bulk carbonate and benthic foraminifera, alkenone SST estimates, and ε_p calculated from bulk carbonate (circles) and benthic foraminifera assuming a constant offset (crosses).

4.1.2 Comparison of 1168 and 1406 ε_p records with published Atlantic records

The overall decline in ε_p through the time interval of our records, is broadly comparable to the trend in published ε_p datasets recently compiled by the paleo CO₂ community (CenCO2PIP Consortium, 2023), which exhibits a long term decrease in the late Oligocene through the transition to the Miocene and to overall low and stable values in the early Miocene (Fig. 4). However, several factors complicate a detailed comparison of our new and the previously published records.

Figure 4 Oligocene to Miocene global long term ε_p trends. Comparison of the new obtained ε_p records with previously published alkenone measurements. All ε_p estimates have been recalculated following the methodology described in text with data source described on Table S2.

The inference of rapid declines is also affected by the age models and the correlation of rapid ε_p shifts among different sites might be hindered by uncertainties in chronology among the different sites. Although all records are presented here on the GTS 2012 (Gradstein et al., 2012), ODP 1168 and IODP 1406 age models rely on Sr isotope stratigraphy (Stoll et al., 2024), whereas ODP 925, DSDP 516 and ODP 608 are based exclusively on biostratigraphic and magnetostratigraphic reversals datums (Curry et al., 1995; Guitián et al., 2020; CenCO2PIP Consortium, 2023). As previous studies document for sites 1168 and 1406, Sr isotopic stratigraphy can adjust age determinations by 0.5 to 1 Myr or even up to 2 Myr in a few cases (Stoll et al., 2024).

Additionally, differences in the absolute value of ε_p among records may also reflect contrasting approaches to the reconstruction of DIC δ¹³C in the different studies. At DSDP 608, the DIC δ¹³C was reconstructed from surface-dwelling foraminifera G. quadrilobatus (Pagani et al., 1999), while at DSDP 516 the Miocene section was estimated from planktic foraminifera and most of the Oligocene samples DIC δ¹³C was determined from fine fraction (Pagani et al., 2000, 2005). Published ODP 925 ε_p has been recalculated here with DIC δ¹³C derived from bulk carbonates of nearby samples, to resolve the previous divergent estimates from planktic and benthic foraminifera (Zhang et al., 2013).

The longest record from DSDP Site 516 exhibits a general ε_p decline from the Oligocene to early Miocene. However, due to lower resolution at this site, we cannot evaluate if there is an abrupt 3 ‰ decline from 27 to 24.5 Ma as seen in sites 1406 and 925. A steep ε_p decline between 21 and 20 Ma in Site 516 may be within age uncertainty of the decrease observed between 20 and 19 Ma at ODP 1168 and ODP 925; additional Sr isotope stratigraphy at Site 516 in this time interval could help test the synchronicity. The late Oligocene at DSDP 516 features a transient 5 ‰ positive excursion in ε_p between 24.5 and 24.9 Ma, which is not reflected at 1406, or 1168 sites. With current information, we cannot assess if this difference reflects age model uncertainty, potentially unresolved coelutions from GC-IRMS chromatography, or aliasing. The early Miocene record at DSDP Site 608 shows a more variable ε_p with a much steeper decline through the early Miocene and higher amplitude variation compared to other sites. The characteristic minimum in ε_p from 18 to 17 Ma is potentially within the age model uncertainty of the 19 Ma minimum in ε_p at 1168 and the 18.5 Ma minimum identified at Site 1406. Updated stratigraphy could contrast more robustly the timing of these events. If there is high amplitude short term variability in ε_p in the early Miocene as in the Oligocene time interval (Fig. 2), there is also the potential for low resolution sampling to undersample high frequency temporal variability and generate aliasing artefacts.

4.2 Potential for size and nutrient effects on ε_p

In addition to CO₂, ε_p may be influenced by changes in cell size and cellular growth rate regulated by light and nutrients. There is no long-term trend in mean coccolith size in these records (Guitián et al., 2020) and estimating the impact of it on the ε_p records shows a negligible effect on long-term trends (Fig. 5). At discrete time intervals of IODP Site 1406, the correction for the size effect reduces most values older than 27 Ma and produces a steeper decrease on ε_p towards the Oligocene Miocene transition.

Figure 5 Timeseries of ε_p at IODP 1406 and ODP 1168 including (a) the measured ε_p record (solid line) and the estimated ε_p resulted once size and temperature effect is applied following Torres et al. (2024), and Stoll et al. (2019) (transparent shadows). (b) Coccolith size record from Guitián et al. (2020); (c) SST estimates (d) Biogenic silica measurements (Guitián et al., 2020 and this study); (e) reference CENOGRID benthic δ¹⁸O curve (Westerhold et al., 2020). Grey squares point ε_p decreasing intervals.

For the statistical model of Eq. (2), it is complex to identify proxy records for any possible effect of nutrient-stimulation of growth rate or changes in the mean light conditions at the depth of growth. In modern spatial gradients in the ocean, these factors are often coupled, so that settings characterized by deep mixing and high nutrient supply rates to stimulate growth, are also characterized by lower mean light levels due to the deeper mixing, both factors lowering ε_p (Hernández-Almeida et al., 2020).

As one possible nutrient indicator, a higher concentration of biogenic silica (bioSi) in sediments may reflect a higher rate of bioSi delivery to the seafloor due to higher export production produced by siliceous organisms (mainly diatoms) in the ocean (Ragueneau et al., 2000). In the modern ocean, regions with abundant dissolved Si in the photic zone are regions also characterized by higher concentrations of macronutrients such as P and N. However, bioSi is an imperfect indicator of past surface nutrient content because coccolithophores have a minimal Si requirement, and Si remineralization in the ocean does not occur at the same rate as soft-tissue nutrients such as N and P. At IODP 1406, bioSi concentrations generally increase from the Oligocene to earliest Miocene, potentially indicating a gradual increase in the concentration of dissolved Si in surface waters at the site (Fig. 5). If the increase in dissolved silica observed at the North Atlantic is correlated to an increase in dissolved P or N, it could contribute to increase in growth rate, and therefore likely increase in biomass and chlorophyll, which would reduce light in the water column, both factors potentially contributing to the observed long term decrease of ε_p. However, the actual correlation between bioSi and ε_p is not that strong (Fig. S2), suggesting that while increased nutrient concentrations could contribute to the long-term evolution of ε_p, the specific steps of ε_p decline are less likely to be driven by increased nutrients and growth rate.

The drivers for increasing bioSi burial rates at Site 1406 are not clear. They could reflect a global increase in nutrient delivery or local effects. Important changes in the rate of continental weathering within the Oligocene- early Miocene are often interpreted from the evolution of radiogenic isotopes of Sr, Li and Os (Misra and Froelich, 2012) including the steep rise in ${}^{\mathrm{87}}\mathrm{Sr}{/}^{\mathrm{86}}$Sr, although the precise origin of the late Eocene and Miocene increase in ${}^{\mathrm{87}}\mathrm{Sr}{/}^{\mathrm{86}}$Sr remains under discussion (Rugenstein et al., 2019). On a global scale, the nutrient delivery may be conditioned by the riverine supply of P from continental erosion and weathering of P containing minerals. Yet, on the time scales examined in our records, much longer than the residence time of P, the net effect on nutrient concentrations depends on the balance of the supply and the nutrient removal in sediments.

While a significant increase in erosion and weathering and nutrient inventory is one mechanism to contribute to the long term decline in ε_p via enhanced algal growth rates, an increase in erosion and weathering can itself contribute to a CO₂ drawdown by CO₂ consumption through silicate weathering and enhanced burial of organic carbon in delta regions (Raymo and Ruddiman, 1992). If the biogenic Si increase at 1406 were representative of a global trend, an increase in nutrient supply may have contributed to ε_p decline through both CO₂ decline and increased nutrient stimulation of phytoplankton growth. A global decline in ε_p solely from increased weathering and nutrient concentrations without a CO₂ decline would require that in the Oligocene, the nutrient release from silicate weathering was less coupled to carbon burial than in the late Neogene. If the periodically glaciated margin of Antarctica is a major locus for increased erosion and weathering in the Oligocene (Reilly et al., 2002), release of nutrients and radiogenic isotopes may have occurred in the continental margins, but with much less organic carbon burial than the modern Himalaya system due to limited terrestrial biomass on Antarctica and temperature and sea-ice limited oceanic biomass production in the marine regions.

On the other hand, the long term trend of increased bioSi is not observed in the Southern ocean Site 1168 (Fig. 5). The available Miocene bioSi at ODP 1168 is stable with no change across the steep ε_p drop from the latest Oligocene to early Miocene. The increasing distance of Site 1168 from the coastline with basin subsidence may have decreased the availability of Si from the early Oligocene through the early Miocene, imparting a local effect superimposed on any potential global trend. However, likely not only Si but also other nutrients would decrease with increasing distance from the coast. If the long term trend in ε_p at both 1406 and 1168 sites were conditioned by increased nutrient availability, faster growth rates, and lower light levels it would require bioSi accumulation rates at Site 1168 to be decoupled from the overall changes in nutrient availability, which we consider less likely. Consequently, we propose that the similarity in trend and magnitude of the long term ε_p decline in both sites (and in tropical Site 925), is more consistent with a global forcing of ε_p, which may be most plausibly driven by a significant decrease in atmospheric CO₂ and CO₂aq.

4.3 Relationship between ε_p and SST and benthic δ¹⁸O

The new ε_p data from Site 1406 and Site 1168 provide the first records of ε_p from the early Oligocene to early Miocene with alkenone unsaturation indices as independent estimations of SST for the precise time intervals of ε_p determination. Since they are biomarkers derived from the same organism, alkenone-derived SST estimates correspond to the same season and growth depth as the alkenone ε_p determinations. There are two processes which may influence the relationship between temperature and ε_p. First, higher temperatures lead to higher phytoplankton carbon fixation rates, decreasing ε_p. Secondly, higher CO₂ would increase ε_p and through radiative forcing lead to warmer global average temperature and SSTs. The expected relationship between ε_p and SST from either process could be obscured by a superposition of temperature effects on growth rate and a climatic correlation of ε_p with mean air temperature.

4.3.1 Million year scale relationships

Across the overall time interval, Site 1406 ε_p is weakly inversely correlated with SST, whereas Site 1168 ε_p is weakly positively correlated with SST (Figs. 6 and S3, Table S3). Our estimation of the growth rate effect due to warmer temperatures shows that it has a very limited impact on the long term ε_p trend, amplifying slightly the long-term excursion in ODP 1168 and imparting a minor increase in ε_p in the late Oligocene 25.5 to 24 Ma in IODP 1406 but otherwise not affecting the sign of the overall trend (Figs. 5, 6 and S3). At our studied sites across the 30 to 17 Ma time interval, the long term average warming of 2–3 °C is insufficient to account for the 7 ‰ decline in ε_p due to a temperature-driven growth rate effect. Unlike alkenone-based SST, the published TEX₈₆ SST record at Site 1168 (Hou et al., 2023a) does not indicate a transition to lower temperatures from the Oligocene to early Miocene, suggesting different temperature trends in the season or depth niches of the different proxy carriers.

Figure 6 (a–b) For paired samples, relationships between the measured ε_p and SST for Site 1406 and 1168. (c–d) Relationship between SST and ε_p, when the temperature-dependent growth rate effect on ε_p has been removed using the observed dependence of from culture experiments of 0.5 ‰ (0.37 ‰–0.95 ‰ range indicated by error bars) decrease in ε_p per 1 °C warming (Torres Romero et al., 2024). (e–f) For paired samples, measured ε_p and benthic δ¹⁸O. Dashed dark line in the plot background shows the regression for the overall dataset. See Fig. S3 and Table S3 for overall statistical relationships.

Over the studied time interval, SST and ε_p similarly show insignificant correlation for the previously published ε_p, records with updated age models and ε_p calculations (Fig. S4), although temperatures estimates are derived from GDGTs which might not reflect the same depth and/or season of coccolithophore growth. Negative covariance is observed at DSDP 608 from 19 to 16 Ma, at DSDP 516 before 21 Ma and for the few samples from 27 to 24 Ma at ODP 925. Some of these intervals feature significant temperature changes of 4 to 5 °C, and therefore the temperature-growth rate effect on ε_p may be significant, and the negative slopes observed in some intervals are consistent with this being the dominant effect. At 1406, during the older intervals of positive correlation of SST and ε_p, potentially the growth rate stimulation due to higher SST was balanced by a decrease of nutrient availability during warmer temperatures as suggested by the bioSi evolution (Fig. S2), whereas during younger time intervals, temperature exerted a dominant effect on growth rate.

Benthic δ¹⁸O was measured in multiple time intervals in Site 1406. Benthic δ¹⁸O has been proposed to reflect global surface temperature (Evans et al., 2024; Hansen et al., 2013) and as such may be less sensitive than SST to regional reorganizations of heat transport. Alternatively, benthic δ¹⁸O has been proposed to be highly sensitive to the areal extent of the Antarctic ice sheet due to its cooling effect on surface ocean temperatures in regions of deepwater formation (Lisiecki and Raymo, 2005; Shackleton, 1987; Bradshaw et al., 2021). If the global surface temperature change translated to changes in surface ocean temperatures at Site 1406 and Site 1168, we would expect the temperature-growth rate effect to generate a direct correlation between benthic δ¹⁸O and ε_p. If the radiative forcing effect on global temperature change were dominant, we would expect an inverse correlation between ε_p and benthic δ¹⁸O. As for SST, only the time intervals older than 25.4 Ma exhibit the inverse correlation expected from radiative forcing, whereas other intervals suggest neutral slope which may reflect the superposition of growth rate and CO₂-radiative effects on ε_p.

4.3.2 Relationships between ε_p, temperature and benthic δ¹⁸O at orbital timescales

In the high resolution sampling between 29.0 and 29.6 Ma, despite a significant 1 ‰ range in δ¹⁸O benthic and δ¹⁸O bulk, we likewise observe no inverse relationship between ε_p and δ¹⁸O benthic, or between ε_p and δ¹⁸O bulk (Fig. 7). We also observe no significant correlation between ε_p and alkenone SST. Because the magnitude of SST variation is small over this time interval, the impact of temperature-stimulated carbon fixation rates is not a significant impact on the relationship between ε_p and SST or δ¹⁸O benthic – a temperature-corrected ε_p record for the 29.6 to 29 Ma interval would still not exhibit an inverse relationship between ε_p and δ¹⁸O benthic (Fig. 7e) as observed in the late Pleistocene glacial cycles (Hernández-Almeida et al., 2023). If ε_p variations are dominantly responding to CO₂, our results suggest that low CO₂ is not contributing to greater ice volume and/or colder ocean temperatures on 100 ky cycles and that the relationship between Antarctic ice growth and CO₂ may be more complex at this time.

Figure 7 For the high resolution 29 to 29.6 Ma interval at Site 1406, relationship between calculated ε_p and temperature, δ¹⁸O from benthic foraminifera, and δ¹⁸O from bulk carbonates. Panels (a) through (c) illustrate the measured ε_p. In panels (d), (e) and (f) ε_p is corrected for temperature-driven growth rate variation, using the temperature effect on ε_p following findings from cultures of 0.5 ‰ per 1 °C (Torres Romero et al., 2024). Vertical error bars illustrate ε_p values when the slope of the correction ranges from 0.37 to 0.95 ‰ per 1 °C warming. Dashed line shows the linear regression for all plotted samples. ^* p value < 0.1, ${}^{**}$ p value < 0.05.

4.4 A Climate and CO₂ paradox from the Oligocene to early Miocene

The long term trends between 30 and 16 Ma based on new ε_p data at two sites and recalculation of previous ε_p studies with uniform methods cannot be attributed to a temperature effect on growth rate and ε_p, nor to changes in the cell size of the alkenone producing community. Both effects are small in magnitude according to the sensitivities observed in cultures and do not alter the long-term trend (Fig. 5). Therefore, the long-term ε_p decline must have a significant global driver, with the most obvious being a decline in pCO₂.

Although the calculation of absolute CO₂ concentrations from ε_p in the Oligocene and early Miocene remains challenging, the logarithmic dependence of ε_p on CO_2[aq] observed in cultures allows us to estimate the relative changes in CO₂ if the sensitivity of ε_p to CO₂ in the Oligocene were similar to modern cultured species using Eq. (7). If we incorporate a temperature correction and apply the 50th percentile estimate of the modern culture ε_p dependence on ln [CO_2[aq]] of 2.66, it implies major changes in CO₂ concentrations, with potentially 4 halvings of CO₂ concentration from 29 to 16 Ma (Fig. 8). Modern General Circulation Models (GCM) summarized by the IPCC estimate climate sensitivity as “very likely” in the range of 2 to 5 °C per doubling or halving of CO₂ (Forster et al., 2021), which if representative for the Oligocene to early Miocene, would imply 8 to 20 °C of cooling of earth's mean surface temperature (6 to 15 °C incorporating the lower confidence interval of modern culture ε_p dependence on ln [CO_2[aq]] of 3.5, which would imply 3 halvings of CO₂). Although ocean is 70 % of the globe and temperature changes are around 1.5-fold less than land temperature (Sutton et al., 2007), such a temperature change of at least 6 °C would be expected to be reflected in paleoceanographic proxies.

Figure 8 Implications of CO₂ as main climate driver. (a) pCO₂ doubling for the discussed sites from ε_p referenced at 29 Ma (Site 608 referenced to maximum at 23 Ma). Solid lines are calculated using the SST corrected ε_p (b) pCO₂ estimates from boron and leaf gas compiled in CenCO2PIP Consortium, (2023) (Sosdian et al., 2018; Greenop et al., 2019; Liang et al., 2022b; Londoño et al., 2018; Liang et al., 2022a; Tesfamichael et al., 2017; Reichgelt et al., 2020; Moraweck et al., 2019; Roth-Nebelsick et al., 2014; Sun et al., 2017; Erdei et al., 2012). (c) Benthic δ¹⁸O global compilation (Westerhold et al., 2020). (d) Schematic representation of main paleoceanographic and paleogeographic changes over the studied time interval for the Northern and Southern Hemisphere. ACC arrow refers to a shallow circulation not synonymous with the deep late Miocene ACC described by Evangelinos et al. (2024). Yellow and blue vertical bands show main warming and cold periods as Middle Miocene Climatic Optimum (MMCO), Mi-1 glaciation, Late Oligocene Warming (LOW) and Middle Oligocene Glacial Interval (MOGI). Antarctic Circumpolar Current. NCW: Northern Component Water.

Similar to phytoplankton proxy records, the available low resolution leaf gas CO₂ records suggest a decline in CO₂ from the mid to latest Oligocene. However, in contrast to phytoplankton proxy records indicating a significant long term decline in CO₂ from the early Oligocene through mid-Miocene, leaf gas CO₂ proxies suggest higher CO₂ in the early Miocene than the Oligocene due to a positive shift across the OMT. Boron isotope-based CO₂ records from 24 to 18 Ma show significant variability with no clear trend, although the higher density of data around the OMT suggests a CO₂ rise from 23 to 20 Ma which may be consistent with the trend observed in the ε_p record at Site 1406, which has the highest resolution for this time interval.

The late Oligocene climate and CO₂ paradox has been discussed based on previously published lower resolution ε_p record from Site 925 (O'Brien et al., 2020). Our new results from two additional sites confirm the steep CO₂ decline through the late Oligocene warming and underscore the paradox. On a global scale, biomarker SST estimates do not show evidence for systematic cooling during the CO₂ decline (Liu et al., 2018; O'Brien et al., 2020; Guitián et al., 2019). If the interpretation of ε_p as a CO₂ decline is correct, it suggests that climate sensitivity was either significantly weaker so that no appreciable change in global mean surface temperature occurred, or that available paleotemperature estimates reflect a significant misinterpretation of measured biomarker signals. During this time the inferred CO₂ decline also coincides with sequence stratigraphic evidence for ice margin retreat in Antarctica (Levy et al., 2019; Salabarnada et al., 2018), and sea level transgressions inferred from estimates of deep sea δ¹⁸O_sw and $\mathrm{Mg}/\mathrm{Ca}$ records (Miller et al., 2020), suggesting a different relationship between ice expansion and CO₂ than characterized the late Neogene.

For the late Oligocene to early Miocene, the Southern Ocean Site 1168 is the only surface ocean temperature record which exhibits a colder early Miocene coincident with the record of a large magnitude pCO₂ decline. This long term cooling is despite the equatorward drift of the site over this time interval (Guitián and Stoll, 2021). If the ODP Site 1168 temperature trend is more representative of global average temperature trends, whereas the long term alkenone temperature record at Newfoundland Ridge Site 1406 and Site 1404 (Liu et al. 2018) is dominated by variations in the heat transport from the Gulf Stream, then the 1168 temperature trend may reflect the signal of radiative greenhouse forcing. Yet, temperature trends at either site may be subject to both global factors as well as regional temperatures, and with only two sites with temperature records paired to ε_p proxy records it is difficult to ascertain which, if any, of the sites may better reflect global temperature forcing.

A decoupling was at one time proposed for the late Miocene based on apparent negligible pCO₂ change and substantial cooling of SST (LaRiviere et al., 2012). Revisions of the alkenone carbon fractionation to CO₂ calibration approaches for low pCO₂ periods have refined the record from the last 15 Ma, revealing clear pCO₂ -SST covariation (Stoll et al., 2019; Rae et al., 2021). However, the Oligocene paradox is not easily resolvable from updated calibration of the ε_p-CO₂ relationship. The late Oligocene paradox arises from an inverse correlation between ε_p and SST reconstructions in regions other than the Southern Ocean such as the North Atlantic, and a lack of correlation between ε_p and the global climate signal in benthic δ¹⁸O trends. The discrepancies between alkenone and published TEX₈₆ at ODP 1168 suggests continued reevaluation of SST proxy interpretation are needed, along with evaluation of the potential influence of changing surface ocean circulation on SST in some locations such as the North Atlantic. Additionally, the divergence of CO₂ trends among ε_p and boron isotopes suggest that further interrogation of ocean chemistry and biogeochemical cycles potentially affecting the growth and physiology of alkenone producers and the calculation of CO₂ from boron isotopes, are crucial to reconcile climate sensitivity to CO₂ in the Oligocene to early Miocene.

5 Conclusions

The new long term alkenone ε_p records from the Oligocene to early Miocene at North Atlantic Site IODP 1406 and Southern Ocean Site ODP 1168 reveal a significant 8 to 10 ‰ shift. The highest resolution Site 1406 record resolve an abrupt 3 ‰ decline from 27 to 24.5 Ma. The long term trend is comparable with previous lower resolution analysis when they are recalculated with the same methodology.

In addition to CO₂, ε_p may be modified by changes in cellular surface area to volume ratio and growth rate regulated by light, temperature and nutrients. However, our assessment of these effects using records of coccolith size and alkenone temperature estimates for exact time intervals of ε_p determination, shows that size and temperature-stimulated growth rate effects have a negligible impact in the long term declining trend. The similarity of ε_p in widely separated sites experiencing contrasting temperature histories strongly suggests a global CO₂ decline as the most likely cause of the declining ε_p. At the same time, our high-resolution sampling reveals significant orbital scale variability in ε_p and underscores the potential for aliasing in low resolution records. Higher resolution ε_p time series, and more precise age models on legacy ε_p records to facilitate more confident comparisons of trends among sites, will provide a better characterization of the key long term trends.

Our results highlight the paradox of complex relationships between CO₂ indicators and SST at both the orbital and multi-million year timescales. The higher resolution sampling between 29.7 to 28.7 Ma shows that orbital ε_p maxima do not coincide with orbital minima in ice volume and/or warmer deep ocean temperature. Similarly, through the late Oligocene warming, CO₂ decline contrasts with evidence for Antarctic ice retreat and evidence of stable or warming SST. The transition from late Oligocene to early Miocene, reaching minimum CO₂ around 19 Ma, is coincident with significant cooling only in the Southern Ocean Site 1168, but not the North Atlantic site which may be more affected by changes in ocean heat transport.