Eocene emergence of highly calcifying coccolithophores despite declining atmospheric CO2

Coccolithophores, a group of unicellular calcifying phytoplankton, have been major contributors to marine carbonate production since the calcite plates that they produce (coccoliths) first appeared in the fossil record over 200 million years ago (Ma). The response of this process to changes in environment on evolutionary timescales remains poorly understood, particularly in warm climates. Here we integrate a dataset consisting of carbon isotope ratios of size-separated coccolith calcite from marine sediments with a cell-scale model to interrogate cellular carbon fluxes and pCO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{{\mathrm{CO}}_2}$$\end{document} through the Eocene (~55–34 Ma), Earth’s hottest interval of the past 100 million years. We show that the large coccolithophores that rose to dominate the oceans through the Eocene have higher calcification-to-carbon fixation ratios than their predecessors while the opposite is true for smaller coccolithophores. These changes, which occurred in the context of increasing ocean alkalization, may have played a role in an apparent positive carbon cycle feedback to decreasing pCO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{{\mathrm{CO}}_2}$$\end{document}. Our approach also provides independent support of multiproxy-based evidence for general pCO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{{\mathrm{CO}}_2}$$\end{document} decline through the Eocene in step with temperature. Together, this challenges the emerging view that a general decline in pCO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{{\mathrm{CO}}_2}$$\end{document} reduces calcification on evolutionary timescales. Highly calcifying, larger coccolithophores emerged as CO2 generally declined through the Eocene, despite cooling leading to lower organic-matter fixation rates, according to size-dependent coccolith carbon isotope analyses and cell-scale modelling

C occolith production and export removes alkalinity from the surface ocean, reducing its capacity to buffer changes in p CO 2 (ref. 1 ). On geologic timescales, coccolith burial in deep-sea sediments removes carbon from the ocean-atmosphere system 2 . From a biogeochemical perspective, carbon fixation in the surface ocean is limited by nutrient availability 3,4 , temperature 5 and light levels 6 ; the amount of calcite produced is then dictated by the particulate inorganic carbon (PIC) to particulate organic carbon (POC) ratio (PIC/POC) of biogenic material. PIC and POC are, respectively, the time-integrated rates of calcite precipitation (R calc ) and carbon fixation (R fix ). R calc /R fix (conceptually distinct from PIC/POC but equivalent at steady state) varies between species of coccolithophore, so changes in species composition through time can alter the amount of calcite transported to the deep ocean [7][8][9][10][11] . At the cellular level, R calc /R fix is a physiologically important parameter because calcification and carbon fixation fluxes have opposing effects on intracellular chemistry. These fluxes also drive the carbon isotope ratio of dissolved inorganic carbon (DIC) in the intracellular environment in opposite directions, which lays the foundation for reconstructing this ratio in ancient coccolithophores through the carbon isotope compositions of their intracellularly formed coccoliths 12 .
The prevailing view, based largely on the morphometry and weight of fossil coccoliths, is that coccolithophore calcification decreased during periods of lowered p CO 2 in the icehouse Neogene and glacial periods 10,[13][14][15][16][17] (although there are counter-examples 11,18,19 ; Supplementary Table 1). There is a degree of ambiguity in the meaning of 'calcification' , whether this refers to the amount of calcite per cell, production rates or the rate of calcite precipitated per unit carbon fixed. In this Article, we target the last of these explicitly by constraining R calc /R fix . We take an approach that integrates carbon isotopes and morphometric data, with a cellular carbon flux model, to simultaneously constrain R calc /R fix across different taxonomic groups and CO 2 levels. The fossil record captures how communities change on evolutionary timescales, and the Eocene in particular, characterized by high but declining temperature 20 and p CO 2 (refs. 21,22 ), provides a unique opportunity to gain this insight in a hothouse climate.
For decades, carbon isotope ratios of coccolith calcite were ignored in the field of palaeoclimatology because they exhibited large and variable deviations from calcite precipitated at equilibrium in the same environment, making them an unreliable target to reconstruct the carbon isotopic composition of DIC. More recently, these deviations, or carbon isotope vital effects (CIVEs), have themselves been shown to contain information regarding rates of calcification 12 . CIVEs are a function of coccolithophore physiology and ambient carbonate chemistry 12,13,23 . Coccolith calcite enriched (depleted) in 13 C compared with equilibrium calcite possesses a positive (negative) CIVE [24][25][26][27] . CIVE expression is driven by the molar ratio of carbon demand to carbon supply and the partitioning of carbon between calcification and carbon fixation 12 . The direction of the CIVE is largely set by R calc /R fix , which differs by a factor of around four in modern species (R calc /R fix is approximately 0.5 in Emiliania huxleyi and 2.0 in Calcidiscus leptoporus 8 ). A high R calc /R fix (typically greater than 1) produces a negative CIVE, while a low R calc /R fix (typically less than 1) produces a positive CIVE 12 . This difference is due to the effect of photosynthesis, which discriminates against 13 C, leaving the residual intracellular pool enriched in 13 C (ref. 12 ). The magnitude of the CIVE is set by the molar ratio of carbon demand to supply (carbon utilization), which increases with cell size and growth rate and decreases with ambient CO 2 (aq). Intracellular allocation of HCO 3 -(ref. 28 ) and variations associated with changing light conditions 23 impose secondary controls.
Cells of different sizes, and characterized by different values of R calc /R fix , express distinct CIVEs under the same CO 2 (aq) regime.

Eocene emergence of highly calcifying coccolithophores despite declining atmospheric CO 2
Coccolith size scales with cell size, so size separation of coccoliths enables the study of different-sized cells [29][30][31] . Analysis of multiple size fractions from the same sample also allows us to define a new parameter, CIVE mean , which is the isotopic composition of coccolith calcite of a specific size fraction minus the mean isotopic composition of all coccolith size fractions from that sample of sediment. This parameter is independent of equilibrium compositions and so can place constraints on R calc /R fix and CO 2 (aq), two key drivers of CIVEs, without independent knowledge of the carbon isotope composition of equilibrium calcite. We present a new coccolith calcite isotope dataset comprising isotope measurements from five size fractions of coccoliths across 37 depths from the Eocene section of Ocean Drilling Program (ODP) leg 208, site 1263, 1,300 km off the coast of southern Namibia (see Methods for site description). The carbon isotope ratios of the individual size fractions reveal rich variability compared with the bulk carbonate record (Fig. 1a). To explore the data, we used a cell-scale carbon isotopic flux model to find the optimal combination of R calc /R fix and CO 2 (aq) that provides the best fit to measured isotopic CIVE mean . Our isotopic data for the Eocene exhibit large deviations from the bulk record at approximately 46 million years ago (Ma) (Fig. 1a). After optimizing for values of R calc /R fix and CO 2 (aq), we were able to successfully model the distribution of measured CIVE mean (Fig. 1b).

Calcification across the Eocene
Eocene coccolithophores were some of the largest of the entire Cenozoic era 32 . Compared with modern, glacial and Miocene periods, Eocene coccolithophores grew under higher temperatures and higher p CO 2 (refs. 21,22 ). In our study, we identified four dominant genera across the Eocene: Coccolithus, Chiasmolithus, Discoaster and Reticulofenestra (Fig. 2). We found that R calc /R fix ratios across size ranges and genera were similar to those of modern-day coccolithophores 8 ; however, the existence of substantially larger cells 29,30 means that more calcite was produced per cell 7,29,30,33 (Fig. 2). On the basis of our modelled values of R calc /R fix and coccolith size measurements from microscopic analysis, we infer that a cell belonging to the Eocene Coccolithus genus (radius ~10 μm) would have consisted of ~1.3 times more calcite per cell than an equivalently sized large modern-day Coccolithus braarudii 34 . We have also been able to place constraints on the R calc /R fix of extinct genus Discoaster, which we conclude was surprisingly low. Nevertheless, owing to the apparently large size of these cells, they may have produced 2.8 times more calcite per cell than modern-day Coccolithus braarudii (Supplementary Information). Until now, attempts to reconstruct R calc /R fix for Discoaster with morphometric approaches have been hampered by a lack of intact coccospheres in the fossil record and the uniqueness of their star-shaped coccoliths, which have no modern analogue ( Supplementary Fig. 1). Isotopic compositions are determined by cellular chemical and isotopic disequilibrium that depends on relative rates of calcification to photosynthesis and thus circumvents these difficulties.
The 3-5 μm size fraction shows a continuous negative CIVE from 56 to 50 Ma (Fig. 1b). Previous studies also noted the presence of the largest negative CIVE in the smallest size fraction during the Eocene 28 ; however, this behaviour has previously proved difficult to interpret. In laboratory cultures, smaller species tend to show more positive CIVEs compared with larger species owing to their lower carbon demand and generally lower R calc /R fix (refs. 12,25,26 ). Our microscopic analysis shows that from 56 to 50 Ma, the smallest size fraction consisted mostly of Coccolithus while the largest size fraction was dominated by the genus Discoaster (Fig. 2). Given our modelled R calc /R fix values for the measured genera distribution and CIVE mean pattern during the early Eocene, we suggest that this negative CIVE mean in the 3-5 μm size fraction was probably the result of the smallest size fraction having a larger R calc /R fix than other size fractions (Fig. 2) as Coccolithus was heavily calcified relative to cell volume.
Our results show that there was a large shift in community composition from Discoaster to dominance by Coccolithus, Chiasmolithus and Reticulofenestra, and that the timing of this coincides with a decline in both CO 2 and temperature 35 . Our inverse modelling efforts suggest that the community shift corresponded to a change from coccolithophores with a low apparent R calc /R fix to those with a higher apparent R calc /R fix in the largest size fraction. Yet in the smallest size fraction, the community shift corresponds to a change from the genus Coccolithus to Reticulofenestra and a shift from high to low apparent R calc /R fix . There is a relative lack of studies that have focused on reconstructing R calc /R fix across the Eocene epoch, although a morphometric study focused on the PETM suggests that species such as Coccolithus pelagicus possessed ratios of ~2 (ref. 7 ), similar to our findings for the Coccolithus genus R calc /R fix . Coeval foraminifera are composed of the genus Acarinina, which probably had symbionts, hence leaving the foraminifera calcite enriched in 13 C; we investigate vital effects associated with foraminifera in the Supplementary Information. Also shown are contemporaneous bulk sediment δ 13 C values, which track the total sum of coccolith and foraminifera calcite and changes depending on their proportion, vital effects and mixed-layer whole-ocean signals. Covariation between the fine fractions and bulk is interpreted to imply that the bulk is largely composed of coccoliths and so their calcite dominates the geochemistry, although we cannot rule out some overgrowth (Supplementary Information). b, δ 13 C values of each sediment size fraction relative to their age interval mean (CIVE mean ). Points are measurements; lines are fitted model values. the error bars on the points represent uncertainty introduced through contamination of neighbouring size fractions (Methods). the coloured envelope surrounding each line represents 95% of model runs over 10,000 Monte Carlo simulations. the coloured vertical regions represent important periods during the Eocene, including the Palaeocene-Eocene thermal Maximum (PEtM), EtM2, early Eocene Climatic Optimum (EECO) and MECO.
It is additionally important to understand how our results fit into the wider context of the global carbon cycle. Evidence suggests that across the Eocene, the calcite compensation depth (CCD) deepened 36 , ocean pH rose 22 and an increase in global weathering [37][38][39][40] probably drove an increased flux of alkalinity to the oceans. Carbonate mass accumulation rates at site 1263 and at other sites across the ocean appear to have declined 41,42 (Supplementary Table 1). The putative deepening of the CCD implies an increase in the saturation state of the deep ocean, which suggests that there was either a whole-ocean decrease in calcite production rates or an increase in delivery of alkalinity to the deep ocean reservoir, possibly via increased respiration-driven dissolution of exported carbonate. Resolving the intensity of respiration dissolution during the Eocene could be addressed in future work by assessing the history of decoupling between the lysocline and CCD.

implications of temperature, light and nutrients on growth rate and calcification
Growth rate is a key parameter in all phytoplankton δ 13 C-based investigations of palaeo-CO 2 but is notoriously difficult to constr ain 12,23,[25][26][27] . In our model, growth rate was assumed to scale allometrically. On million-year timescales, the effects of variations in nutrient conditions, light and temperature on cellular growth rate, photosynthesis and calcification 43 are difficult to predict because of evolutionary change and species succession. We therefore assume that growth rate is constant for individual cell sizes across genera and incorporate uncertainties within our model (Methods). A constant growth rate is qualitatively supported by strontium-to-calcium ratios (Supplementary Information). Lab-based studies have shown that coccolithophores probably exhibit a species-specific optimum-like response to variations in temperature 5 , light and carbonate chemistry 44 . For an observed CIVE, a net decrease (increase) in the growth rate of a size fraction would require a correspondingly lower (higher) predicted ambient CO 2 concentration. Future improvements may be made with an allometric formulation of growth-rate relationship that includes coefficients for temperature, light and nutrient conditions, which would allow for such variations to be quantitatively investigated within our model. Within our sensitivity analysis ( Supplementary Fig. 2), we found that, although variations to assumed growth rate altered absolute values of R calc /R fix and relative CO 2 , given here as the ratio of modelled output CO 2 (aq) to modelled input CO 2 (aq) (see Methods for further derivation of relative CO 2 and reasoning for its use over absolute values), trends remained the same.
In our study, we assume that R calc /R fix is constant for each genus and size fraction throughout our record. Predicted R calc /R fix values therefore reflect the best fit for the entire Eocene. This is a limitation to our approach: the long-term environmental change across the Eocene may be associated with variations to R calc /R fix within a genus and size fraction, which, due to our assumption that R calc /R fix is constant across our interval, cannot be captured. Previous culture studies on multiple species have found complex responses for variations in R calc /R fix given changes to environmental parameters 45 , but the translations of these results to evolutionary timescales is not clear. Changes in R calc /R fix within a genus and size fraction could potentially be constrained if this modelling approach were paired with an independent high-quality CO 2 (aq) record as an input and high-resolution record of δ 13 C DIC .

Eocene CO 2 levels
Our model implies that relative CO 2 (aq) was higher in the early Eocene and decreased across the epoch as CIVE mean emerge and increase in magnitude, in line with previous data 21 (Fig. 3). Given site 1263 was probably close to equilibrium with the atmosphere (Methods), modelled CO 2 (aq) can be assumed to reflect atmospheric p CO 2 . Unlike the use of other carbon isotope-based p CO 2  proxies [46][47][48] , the use of mean-normalized model predictions fitted to mean-normalized measurements makes our approach internally dependent on the isotopic composition of coccolith calcite and independent of seawater carbon isotopic values. Independence comes at the cost of one degree of freedom per depth interval (total of 37) and results in a solution that is non-unique in absolute CO 2 (aq) but robust with respect to relative change in CO 2 (aq) and R calc /R fix (Supplementary Fig. 3). Loose constraints on the absolute value and change in p CO 2 are provided by estimates of seawater carbon isotopic values from foraminiferal calcite ( Supplementary  Fig. 4). When carbon utilization is low (particularly at high external CO 2 (aq)) the effect of CO 2 (aq) on CIVE expression is small. During Eocene hyperthermals (for example ~40 Ma), the model therefore results in a greater spread of possible CO 2 (aq) values. Qualitatively, the convergence of CIVE mean at 53.5-54.0 Ma indicates higher p CO 2 , given reduced CIVEs are observed at higher [CO 2 (aq)] in culture experiments 12 , which is supported by our modelled predictions (Figs. 1 and 3). Long-term temperature change across the Eocene was considered when calculating p CO 2 atm across the interval via Henry's law (Supplementary Fig. 4). Previous estimates of Eocene CO 2 vary between ~500 and 2,000 ppm (refs. [20][21][22][49][50][51] ) but agree on the presence of a long-term decline, with substantial departures at 40 Ma and 54 Ma (refs. 21,22 ), in line with our results (Fig. 3). Our predicted peak in relative CO 2 (aq) during the early Eocene approximately 53.5-54.0 Ma may be attributed to Eocene Thermal Maximum 2 (ETM2), a rapid period of warming [52][53][54] . Further, high relative CO 2 (aq) at ~40 Ma is probably associated with the middle Eocene Climatic Optimum (MECO 55 ). The MECO is believed to represent a period of elevated volcanic outgassing, decreased silicate weathering due to a decreased weatherability after the extended warm Eocene 56 and rapid shallowing of the CCD 36 . Before the MECO, our relative CO 2 (aq) curve shows a highly significant (P < 0.001) negative correlation with global δ 18 O benthic data (Fig. 3b), a widely accepted proxy for deep-sea temperature (in the absence of ice sheets) in the Eocene that reflects longer-term climatic change. The observed covariance between reconstructed CO 2 (aq) and global δ 18 O benthic data supports the robustness of our findings. From 41.5 Ma to 38.5 Ma, we find that the relative CO 2 (aq) curve shows no correlation to global δ 18 O benthic . Data from 41.5 Ma to 38.5 Ma (4 of 37 data points) were excluded from the regression analysis. The exclusion is justified given the potential decoupling between δ 18 O benthic , atmospheric CO 2 and carbonate storage within the ocean during the MECO owing to an unstable nascent cryosphere 50 . Variations in ice volume at high latitudes and altitudes from 41.5 to 38.5 Ma rather than temperature could cause a deviation in the relationship between δ 18 O benthic and CO 2 .
The question thus remains: what parameters drive long-term adaptive changes in calcification? In addition to declining temperature and CO 2 (aq), the Eocene was probably characterized by an increase in ocean alkalinity 37,38,40 , nutrients and pH 22 . These conditions appear to have selected for large cells with higher R calc /R fix , with the inverse being true for the smallest cells. For the largest cells, whose small surface-area-to-volume ratios may make them particularly vulnerable to self-acidification of the cytosol by calcification, increasing surface ocean pH 22 may have resulted in a more favourable gradient for the efflux of H + ions and thus enabled the emergence of large, more heavily calcified coccolithophores 57,58 . This effect may have been less pronounced in the small cells, which may therefore have been limited by other environmental parameters. To fully elucidate the environmental drivers of calcification change in coccolithophores, further work to understand why coccolithophores calcify will be necessary 58 .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of , defined as model output CO 2 (aq) relative to model input CO 2 , is derived from our optimized model with size-differentiated isotopic constraints (see Methods for detailed explanation of relative CO 2 ). the coloured markers represent previously published pCO 2 data 21,22 across the Eocene (right-hand y axis). the black circles represent the median relative CO 2 (aq) over 10,000 iterations of our Monte Carlo optimization (left-hand y axis), with error bars representing plus and minus 2-sigma of all realisations for each age point (95% of realizations). the black line represents a three-point moving average of our output, with the grey shaded region representing the moving average of plus 2 sigma and minus 2 sigma. Our modelled low relative CO 2 (aq) during the middle Eocene coincides with a step-wise increase in proxies linked to an increase in the rate of chemical weathering and therefore CO 2 draw-down [37][38][39][40] . the dashed black line during the MECO joins points without a moving average. Using a moving average at this time greatly skews the trend line given the large, and uncertain, deviation to higher relative CO 2 . Also shown is the approximate interval of the EECO. b, the five-point moving average of relative CO 2 (aq) is plotted against the five-point moving average of δ 18 O benthic . the five-point moving average was used to avoid noise in the data and evaluate the longer-term trend. Error bars represent one sigma uncertainty on reconstructed CO 2 (aq) values. the equation of the line has been overlaid on the figure. the four points in red represent the points from the MECO and younger, which have been excluded from the P-value calculation (P < 0.001). Without the MECO data included, our relative CO 2 (aq) and δ 18 O benthic show a robust negative correlation. the warmer early Eocene, characterized by more negative values of δ 18 O benthic , has on average higher values of relative CO 2 (aq). c, the slope of relative CO 2 (aq) with δ 18 O benthic for each of the 10,000 Monte Carlo simulations was calculated and plotted as a histogram. Only one value was greater than zero.

Methods
Sample selection and preparation. The sediment samples were from ODP leg 208, site 1263. Site 1263 has remained at a relatively constant latitude for most of the Cenozoic 59 . Ocean general circulation models suggest reduced Atlantic upwelling during the Eocene, with the Walvis ridge probably outside of any significant zones of upwelling 60,61 . We thus posit that site 1263 was in equilibrium with the atmosphere. Site 1263 was drilled from a current water depth of 2,717 m, with the palaeodepth during the Eocene believed to be around 1,500 m (ref. 59 ). Thirty-seven samples were selected at approximately 0.5 Myr intervals through the early and middle Eocene.
Sediment samples were placed in 100 ml of pH-neutral de-ionized water and placed on an orbital shaker overnight for disaggregation. Samples were then filtered through 64 μm, 38 μm and 20 μm sieves to remove large particles such as foraminifera and their fragments. The <20 μm fraction was separated into four smaller fractions with microfiltration 62 : 3-5 μm, 5-8 μm, 8-10 μm and 10-20 μm. The 10-20 μm fraction was further separated into 10-15 μm and 15-20 μm fractions via the settling method 63 . Through size separation and microscopic analysis, we constrained the average coccolith length for each genus across all size fractions. The coccolith length was subsequently used to determine cell size through established morphometric relationships ( Supplementary Information) 7,29 . Nannofossil counts were conducted over a gridded 400 μm 2 area. Nannofossil relative abundance is presented as percentage contribution of coccoliths to total counted and averaged across two counts. The microfiltration technique was assumed to be imperfect; we therefore ascribed an uncertainty to the measured coccolith carbon isotopic values relative to an international standard (δ 13 C) assuming contamination from neighbouring size fractions. The 15-20 μm fraction and the 3-5 μm fraction were assumed to have a 25% contamination from the 10-15 μm fraction and the 5-8 μm fraction, respectively, in line with previous studies 64 . The remaining size fractions were assumed to have a 12.5% contamination from the bounding size fractions. Scanning electron microscope images show good preservation of coccolith material; however, we address the possibility of slight overgrowth on individual coccoliths/nannoliths in the Supplementary Information. Age model. The age of each sample has been calculated on the basis of three astronomically tuned age models for site 1263 that best represent the middle and late Eocene 65-67 . The age model from 56 Ma to 50 Ma is based on a regression model of the biostratigraphic datums to best fit the data around the PETM, with our oldest point believed to be during peak PETM, given its isotopic signature.
Isotopic measurements. The δ 13 C lith was determined from ~150 μg of size-separated sedimentary material. In total, 178 samples were processed (5 size fractions × 37 depths = 185, minus 7 samples lost due to machine error). Material was analysed via a Delta V Advantage isotope mass spectrometer fitted with a Gas Bench II in the Department of Earth Sciences at the University of Oxford. An isotopic measurement was made for each of the five size fractions from each time interval. Samples were calibrated to PeeDee Belemnite via the international NBS-19 standard. All measurements are expressed relative to the Vienna PeeDee Belemnite standard, and the analytical precision was better 0.04‰ across repeat measurements for δ 13 C (1σ). The final error on isotopic measurements encompasses both analytical error and error introduced from the microfiltration technique (Fig. 1).

Development of the forward model.
We developed a forward model that used coccolithophore taxonomic abundance (resolved to genus level) and a previously published cellular isotopic flux model 12 to generate abundance-weighted predictions of size-fraction-specific CIVE mean , with individual taxon-specific CIVE mean calculated as a function of R calc /R fix and CO 2 (aq). Our approach separates the dependence of the CIVE from the isotopic composition of DIC by subtracting the mean CIVE for each depth (age point) from both the modelled predictions and the data. We use an iterative optimization approach to search parameter space for 20 size-and genus-specific values of R calc /R fix (5 size fractions within the 4 dominant genera) and 37 (age-or depth-specific) values of CO 2 (aq) to maximize the fit of model predictions to the mean-normalized CIVE record (See Supplementary Fig. 5 for full work flow). Our approach allowed changes in relative CO 2 , modelled output CO 2 (aq) relative to model input CO 2 (aq), to be constrained independent of δ 13 C DIC . We also used a foraminiferal record of δ 13 C DIC to place additional constraints on absolute p CO2 (Supplementary Information). We used a Monte Carlo approach to explore the sensitivity of constrained values to model inputs ( Supplementary Information).
The forward model predicts the δ 13 C of each size fraction as a weighted sum of the predicted δ 13 C for taxa that contribute to each size fraction and their observed relative abundances. Predicted δ 13 C values for each size and taxon-specific unit (STU; 4 genera across 5 size fractions giving 20 in total) were calculated from a previously published model as a function of ambient carbonate chemistry (all STUs are assumed to experience the same ambient environment at each time point) and STU-specific parameters: R calc /R fix , cell size and growth rate (assumed constant for each STU throughout the studied interval). Growth rate was assumed to follow an allometric scaling relationship with cell size 8 , cell size was determined from coccolith size 7,29,33 and abundances were counted directly (model inputs are summarised in Supplementary Table 2). The assumption of constant growth rate for each STU is supported by Sr/Ca data (Supplementary Fig. 6). Parameters fitted to culture data in the original paper remain unchanged 12 . Carbon speciation in ambient seawater was calculated via the Seacarb package for R 68 with an assumed constant calcite saturation state. The remaining parameters, CO 2 (aq) and R calc /R fix (referred to as RR in the following mathematical representation), defined the parameter space to be explored. Initial values for CO 2 (aq) and R calc /R fix were initialized from a random uniform distribution (Selecting initial CO 2 and R calc /R fix ).
For each STU-specific RR i (i ∈ 1 : I, I = 20 STUs) and each age-specific CO 2 (aq) j (k ∈ 1 : K, K = 37 depths), the isotopic flux model was used to calculate a δ 13 C value, D ij (I × K = 740 in total). Modelled estimates of size-fraction-specific δ 13 C, E jk ( j ∈ 1 : J, J = 5 size fractions) then involved weighting each age and STU combination by its relative contribution to the analysed calcite for each size fraction. To do this, D was combined with a matrix of the fractional abundance of each genus in each size fraction (B, with elements B ik ) on the basis of data from microscope counts, a matrix of genus coccolith weight (C, with elements C ik ) calculated from previously published morphometric relationships 7,33 and a sorting matrix to produce a J × K (5 × 37) matrix of age-and size-fraction-specific δ 13 C values (E): A, B, C, D and E are all matrices; ⊙ represents the mathematical operation of element-wise multiplication and ⊘ represents element-wise division, and A (5 × 20) is given by: The overall result is to transform D (STU-specific δ 13 C I × K matrix) to E (a J × K matrix that represents the weighted δ 13 C values for each size fraction that would have been measured at each time interval).
The modelled δ 13 C mean for each time interval was then subtracted from each size fraction δ 13 C at the corresponding time interval to remove the dependence on absolute isotopic compositions and generate a mean-centred value of CIVE (F) for each size fraction: where E 1:J,k is the mean for modelled δ 13 C across all five size fractions at each age point; F represents a J × K matrix of mean-centred model-predicted δ 13 C values subsequently used for the fitting of R calc /R fix and CO 2 .
Optimization. The measured data were treated similarly to remove the dependence on δ 13 C DIC and generate a mean-corrected measured isotope matrix, H: where G jk is the absolute δ 13 C calcite measured for each size fraction at each depth and G 1:J,k is the mean δ 13 C calcite across all five size fractions at each age point. The misfit function between the mean-corrected data (H) and the mean-corrected model-predicted values (F) is given by: where n is the total number of measured and modelled pairs (J × K minus any missing values). There were 7 missing values across 185 size and depth combinations, giving a total of 178 data points to constrain 57 parameter values (20 R calc /R fix + 37 CO 2 ). An iterative optimization approach was taken to find the best CO 2 (aq) and R calc /R fix combination to generate model predictions (F) that fit the measured mean-centred CIVE data (H). We used a two-step iterative optimization to minimize the offset: step one consisted of optimizing for R calc /R fix given an initial constant CO 2 (aq); step two consisted of optimizing for CO 2 (aq) using the values of R calc /R fix determined in step one. Taking the fitted values of CO 2 (aq) for the new initial value, the optimization process was repeated until R calc /R fix and CO 2 (aq) were stable. An unconstrained Nelder-Mead optimizer (fminsearch, Matlab) was used to determine fitted values of CO 2 (aq) and R calc /R fix that reduced the misfit between the measured and modelled CIVE mean (according to equation (4)) for each iteration. Both CO 2 (aq) and R calc /R fix were allowed to be only positive values.
Selecting initial CO 2 and R calc /R fix . The model outputs of R calc /R fix and CO 2 (aq) using CIVE mean are dependent on their initial values. However, relative change is robust across depths and across STUs (Supplementary Fig. 3). These parameters described a valley of stability in parameter space whereby the average output CO 2 (aq) trades off against the average output R calc /R fix . We therefore present CO 2 as a relative change where relative CO 2 is given by the optimized CO 2 (aq) output divided by the initial prescribed CO 2 (aq) value. Differences in the trends in CO 2 between different initial values were minor; however, to non-arbitrarily choose initial values of CO 2 and R calc /R fix , we introduced an additional optimization step that used an independent estimate of δ 13 C DIC derived from coeval planktic foraminifera ( Supplementary Fig. 4).F From a log-uniform distribution between 0.5 and 10, 10,000 initial values of R calc /R fix were selected, while an initial 10,000 CO 2 (aq) values were selected from a uniform distribution between 20 and 150 μmols kg −1 , with external carbon conditions allowed to reach concentrations as low as 5 μmols kg −1 . The optimization outlined in equations (1)-(3) was undertaken for each of the 10,000 combinations, and a fitted array of CO 2 (aq) and R calc /R fix values was generated for each combination. The fitted values of CO 2 and R calc /R fix were then passed to the forward model to generate 10,000 E matrices of modelled absolute CIVEs.
The δ 13 C of planktic foraminifera were analysed for ten depths with the associated δ 13 C DIC determined assuming a constant vital effect (P; Supplementary  Information). The δ 13 C DIC was added onto E j:J for the corresponding depths, and the misfit between the mean raw size fraction data (G) and the predicted mean of the E matrices + δ 13 C DIC were calculated: where foram depths is a set where each value is the depth at which planktic foraminifera were picked and measured, P is the vector of δ 13 C DIC values and N is the total number of size-fraction-specific δ 13 C values from the ten samples where we analysed coeval foraminifera.
Comparison of the modelled CIVEs with δ 13 C DIC provides constraints on the absolute magnitude of R calc /R fix and CO 2 (aq) (and p CO2 ; Supplementary Fig. 4). However, estimating δ 13 C DIC from planktic foraminifera relies on poorly constrained vital effects, which compromises the elegance of our approach. We thus present relative CO 2 (aq) changes based on an internally consistent analysis of coccolith calcite (independent of δ 13 C DIC ) and consider estimates of absolute CO 2 (aq) to be only loosely constrained by the additional consideration of the δ 13 C of foraminiferal calcite.
Monte Carlo exploration of model output. A Monte Carlo approach was employed to further explore how uncertainties in growth rate, coccolith abundance data, coccolith weight and external carbonate saturation state impact the results over 10,000 runs. For each iteration, values were resampled from their uncertainties (Supplementary Table 2). In addition, the initial target value for H was selected from the uncertainty ascribed due to contamination during the filtration process (Sample selection and preparation). The distribution was assumed normal around its mean. Initial R calc /R fix was fixed at 1.45, and the initial [CO 2 (aq)] at 29 μmol kg −1 , the values that provided the best fit of measured CIVEs with modelled CIVEs relative to δ 13 C DIC . Assuming an average p CO2 of 1,000 ppm, pH of 7.7 and temperature of 25 °C across our time interval 21,22 , the corresponding average [CO 2 (aq)] calculated via Seacarb in R 68 would be ~29 μmol kg −1 , which further validates our choice. Final values of R calc /R fix and relative CO 2 (aq) presented in Figs. 2 and 3 represent the output of R calc /R fix from the Monte Carlo simulation over 10,000 iterations.
We also performed a sensitivity analysis on our Monte Carlo simulation to determine the degree of uncertainty that each of our input variables had on our final predicted values of CIVE, relative CO 2 (aq) and R calc /R fix over 100 iterations ( Supplementary Fig. 2). For each of the sensitivity analyses, uncertainty on all but the target variable was assumed to be 0, while the uncertainty on the target variable was ascribed as listed in Supplementary Table 2. We found the model similarly sensitive to all parameters within the ranges investigated. The robustness of our outputs to variations in our model assumptions such as the validity of an allometrically scaled growth rate and an assumed constant R calc /R fix was also investigated. Our results show that the trends in CO 2 (aq) and R calc /R fix were robust across growth rates and that the model is unable to fit measured CIVE mean given a constant R calc /R fix (Supplementary Fig. 7). When growth rate was fixed at one for all sizes, much larger values of relative CO 2 (aq) were required to fit the modelled data. Finally, we investigated the effect of changes in carbonate chemistry parameter space given variations in Ca 2+ concentration across the Eocene ( Supplementary  Fig. 7). While marginal differences were apparent, the trends in CO 2 (aq) and R calc /R fix were robust. We opted not to include the change in Ca 2+ over time to avoid imposing trends into our data that were not internally derived from the measured values.
The model produced a narrow range of fitted values of R calc /R fix when the genus contributed a large proportion to the total sediment size fraction (Fig. 2). If a genus is underrepresented in a specific size fraction, the possible range of R calc /R fix is much wider since it contributes minimal carbonate within the size-fraction δ 13 C signal. The R calc /R fix values of Chiasmolithus in the 3-5 μm and the 5-8 μm size fractions and Reticulofenestra in the 15-20 μm size fraction are poorly constrained and thus not presented.