First absolute seasonal temperature estimates for greenhouse climate from clumped isotopes in bivalve shells

ABSTRACT


Introduction
Seasonal extremes were of vital importance for the evolution and distribution of life over geological history 1 .The effects of greenhouse warming on seasonal variability in temperature and the hydrological cycle are still poorly constrained, while being of considerable interest for projecting future climate and its impact on the ongoing biodiversity crisis 2,3 .Reconstructions of deep time (pre-Quaternary) greenhouse periods yield valuable insights into the dynamics of warm climates and the ecological response to forcing mechanisms such as rising atmospheric CO2 levels 4,5 .Accurate reconstructions are imperative to evaluate climate model predictions under dissimilar climate states 6 , in particular seasonal range for which there is little quantitative evidence to constrain climate model behavior.The warm, ice free Late Cretaceous period presents a valuable reference to assess seasonal variability under greenhouse conditions 7,8 .
However, the reliability of past seasonal reconstructions is undermined by poorly constrained assumptions.This hampers our understanding of past warm climates and hinders accurate evaluation of climate models 13,14 .Two potentially significant biases resulting from assumptions underlying SST reconstructions are seasonal bias and bias resulting from assumed seawater composition (δ 18 Osw bias).
Seasonal bias occurs if a proxy is interpreted as representing annual mean conditions but is in fact biased to a particular season.Since fossil species producing the material that constitutes SST archives may not have a close modern relative for proxy calibration 16 , uncertainties about their growth seasons may unpredictably bias reconstructions.This bias limits our understanding of the behavior of greenhouse climates, which further leads to misinterpretation of model-data comparisons of past warm climates 15 and hinders the use of paleoclimate data for informing future climate predictions.Seawater oxygen isotope composition (δ 18 Osw) forms an important input parameter into the widely used carbonate δ 18 Oc temperature proxy 17 , but remains poorly constrained across geological timescales 18,19 .Biases in assumed δ 18 Osw composition thus undermine SST reconstructions.
The advent of carbonate clumped isotope (Δ47) SST reconstructions on a seasonal scale promises to eliminate these two biases 20 .The clumped isotope thermometer yields accurate SST reconstructions independent of δ 18 Osw assumptions 21,22 .This technique also allows the reconstruction of δ 18 Osw, yielding information about the (local) hydrological cycle, an important aspect of climate rarely constrained in deep time, and rectifying bias in the popular carbonate δ 18 Oc temperature proxy.Recent advances in clumped isotope instrumentation and standardization have reconciled previous inter-lab disagreements and shown that many carbonate paleoarchives (e.g.foraminifera, bivalves and eggshells) conform to the theoretical Δ47 temperature calibration 23,24 (see Supplementary Methods and Discussion).The large sample sizes required for individual Δ47-based temperature estimates (>2 mg) have complicated paleoseasonality reconstructions using this accurate method 25 , but a recently developed statistical approach enables its use for seasonality reconstructions 20 .
Here we use clumped isotope analyses on microsampled (~100 µg) profiles through fossil bivalve shells to obtain, for the first time, absolute SST and δ 18 Osw seasonality reconstructions of a greenhouse climate.We apply this new method on well-preserved oyster (Rastellum diluvianum and Acutostrea incurva) and rudist (Biradiolites suecicus) shells from Campanian (78.1±0.3Ma 26 ) coastal localities of the Kristianstad Basin in southern Sweden (46±3°N paleolatitude 27 ; see Fig. 1 and METHODS).We further supplement these reconstructions with fully coupled climate model simulations of the Campanian greenhouse (see METHODS) to explore their implications for the "equable climate" hypothesis.

Results
All specimens showed clear seasonal δ 18 Oc fluctuations of −2.0-0.0‰ in R. diluvianum, -2.0-0.0‰ in A. incurva and -2.7‰--1.0‰ in B. suecicus on which shell chronologies were based (see Methods).The assumption that periodic δ 18 Oc fluctuations reflect seasonality is demonstrated to be a valid basis for constructing intra-shell chronologies in nearly all modern environments 20 .
Detailed step-by-step results of the data processing routine are shown in Supplementary Methods and 20 .Fig. 2 and Table 1 show monthly Δ47, SST and δ 18 Osw reconstructions for each specimen.Uncertainties at the 95% confidence level on monthly SST vary between 1.8 and 4.2°C owing to variable monthly sampling density related to intra-shell growth rate variability (Fig. 2).
While growth rate varied through the year (Fig. 2A), all monthly time bins in each specimen contained enough datapoints to allow separate monthly SST and δ 18 Osw reconstructions.
We compare our reconstructed SST with local and global Campanian SSTs modelled using the HadCM3L model, one of the most developed paleoclimate models to date and part of the IPCC intercomparison assessment reports 6,14 .We present global Campanian latitudinal gradients in summer, winter and MAT (Fig. 3A) as well as monthly SST in the Boreal Chalk Sea (Fig. 3B) for both 2× and 4× preindustrial atmospheric pCO2 simulations (see Methods).Model results are summarized in Supplementary Data 5.The modelled Campanian latitudinal SST gradient (difference between tropics and high-latitude MAT; 26°C in both simulations) resembles the modern (25°C gradient).Modelled global mean Campanian SST seasonality (difference between warmest and coldest month) is lower (6.6°C)than that of the modern ocean (8.6°C) under 2× preindustrial pCO2 conditions and similar to the present (8.2°C) in the 4× preindustrial pCO2 simulation, in disagreement with the hypothesis of reduced seasonality during greenhouse conditions.Campanian modelled MAT is ~18°C and ~22°C under 2× and 4× preindustrial atmospheric pCO2, respectively, compared to ~14°C in the modern ocean (NOAA, 2020), yielding an equilibrium climate sensitivity, or global warming per doubling of atmospheric CO2 concentration, of ~4°C 15 .Specifically, simulated seasonal SST ranges in the Campanian Kristianstad Basin of 7±3-20±2°C and 12±2-26±2°C for 2× and 4× preindustrial atmospheric pCO2 forcing, respectively, are significantly warmer than in present day southern Sweden (3±0.8-17±0.4°C 28).

Discussion
Our novel Δ47-based monthly SST and δ 18 Osw reconstructions from A. incurva and B. suecicus are statistically indistinguishable from 4× preindustrial pCO2 simulations (p>0.05) and significantly warmer than the 2× preindustrial pCO2 simulations (>4°C higher MAT, p<0.05) of local SST seasonality (p>0.05;Fig. 3).Higher (p<0.05)SST (+4-5°C) and δ 18 Osw (+1.0-1.5‰) in R. diluvianum are likely caused by local differences in its shallower, inter-tidal (< 5 m) environment 27 .Temporary areal exposure during low tides could have elevated temperatures and δ 18 Osw recorded in R. diluvianum year-round by direct sunlight and evaporation, as in modern inter-tidal oyster species 29 .By comparison, the deeper (5-15m) subtidal environments of A. incurva and B. suecicus were unaffected by these processes and may have received more water with an open marine δ 18 Osw signature (closer to the -1‰VSMOW assumed for the ice-free Cretaceous 30 ), especially in winter.These local environmental differences are not resolved in the climate model simulations but show the unprecedented detail of local SST and δ 18 Osw reconstructions from clumped isotope analyses on bivalve shells (see Supplementary Discussion).The ~1‰ δ 18 Osw seasonality shows that summers in the Campanian Kristianstad Basin either experienced excess evaporation, which increases δ 18 Osw by preferentially removing isotopically light seawater, or reduced precipitation, which supplies isotopically light meteoric water, reducing δ 18 Osw.Both processes lead to comparatively dry summers and wet winters.
Strong seasonal fluctuations in δ 18 Osw (up to 1.3‰ in B. suecicus) and regular deviations from the commonly assumed -1‰VSMOW δ 18 Osw value 26 lead to large differences (up to 8.9°C in R. diluvianum) between SST estimates based on Δ47 and the widely used δ 18 Oc proxy (Fig. 2).
The risk of assuming constant δ 18 Osw is even more clearly illustrated by significantly (+3.5-6.0°C)higher δ 18 Oc-based seasonal temperature reconstructions for B. suecicus compared to A. incurva, while both specimens grew under similar SST seasonality conditions (Fig. 2B).Similarly, δ 18 Ocbased temperature reconstructions of A. incurva and R. diluvianum are indistinguishable, while the paleoenvironment of R. diluvianum was 4-5°C warmer year-round (Fig. 2B), illustrating that the constant δ 18 Osw assumption is only valid in settings with negligible δ 18 Osw seasonality and where δ 18 Osw is known.Low-latitude Tethyan SST seasonality reconstructions based on rudist δ 18 Oc 10 agree with model simulations, which may indicate that δ 18 Osw seasonality is less important in open marine settings, although data-model agreement is by no means solid evidence for correct δ 18 Osw assumptions, which should always be independently verified (Fig. 3A).Our findings corroborate previous Δ47-based and proxy comparison studies which also report a significant cold bias (~-8°C) in δ 18 Oc-based SST reconstructions due to inaccurate δ 18 Osw assumptions 8,31 .However, these studies did not account for seasonal biases.
Seasonal variability in growth rates in all specimens (Fig 2A ) illustrates how bulk sampling of bio-archives can lead to significant biases in MAT reconstructions compared to our more accurate estimates of MAT as an average of Δ47-based monthly SST.In this case, slower summer and autumn growth (months 4-7 in Fig. 2A), especially in A. incurva and B. suecicus, would cause bulk analysis of shell material to underestimate MAT.Indeed, our Campanian mid-latitude SST ranges (~15-27°C, MAT of 20°C) are significantly higher than previous reconstructions of the same paleolatitude based on fish tooth δ 18 Oc (15-20°C 11 ), chalk δ 18 Oc (12-15°C 12 ), bulk mollusk Δ47 (5-12°C 32 ), TEX86 (15-20°C 8 ) and sub-annual mollusk δ 18 Oc (15-22°C 10,26 ; Fig. 3A).All these reconstructions are potentially affected by either seasonal or δ 18 Osw bias, or both.Our more accurate reconstructions of δ 18 Osw and SST on a seasonal scale aid in evaluating these biases and correct for them by combining long-term MAT reconstructions with snapshots of climate on the seasonal scale.
Given the increase in frequency and duration of growth stops in modern mollusks with increasing latitude 33 , seasonal biases are likely more common in higher latitude environments.The accuracy of our new method for SST reconstruction and the remarkable agreement between Δ47based SST ranges and our climate model strongly suggest that the average seasonal range reconstructed from our three specimens ( Robust agreement between our reconstructions and the 4× preindustrial pCO2 model simulation down to the monthly scale provides strong evidence for considerably warmer (~8°C) higher latitudes during the Late Cretaceous greenhouse compared to the present day.Significant disagreement of summer, winter and annual SST reconstructions from every specimen in this study with the 2× preindustrial pCO2 simulation strongly favor warmer (4× preindustrial pCO2) climate conditions (see Supplementary Data 4).Bio-archives from mid to high latitudes are likely much more sensitive to δ 18 Osw and seasonality bias than low-latitude records, contributing to the flawed paradigm of shallow latitudinal temperature gradients during greenhouse climates.Instead, our results concur with the recent trend of converging data and model reconstructions yielding modernscale Late Cretaceous latitudinal temperature gradients 13 , thereby challenging the hypothesis of "equable climate" during greenhouse periods 9 .Moreover, our unique absolute monthly SST reconstructions and model simulations corroborate growing evidence against the hypothesis of reduced temperature seasonality in greenhouse climates 30 .Future work should aim to further test these hypotheses by applying the clumped isotope seasonality method on bio-archives from a range of latitudes in greenhouse climate periods.Results from B. suecicus represent the first Δ47-based SST reconstructions from rudist bivalves, introducing an abundant archive for accurate Mesozoic SST seasonality reconstructions with which these new insights can be evaluated.

Conclusions
Our new absolute temperature seasonality reconstructions merit critical evaluation of classical paleoclimate records that risk bias, such as those based on δ 18 Oc (assuming constant δ 18 Osw 7,10 ), bulk analyses of fossil material with growth seasonality (e.g.mollusks and brachiopods 31 ) or a fixed growth season (e.g.planktic foraminifera 36 ) and organic proxies that may be seasonally biased (e.g.TEX86 and U k' 37 5,37 ).In addition, our monthly δ 18 Osw reconstructions for the first time allow evaluation of local seasonality in the hydrological cycle from accretionary bio-archives, revealing dry summers and wet winters in the Campanian Kristianstad Basin.This unique advantage of Δ47-based seasonality reconstructions enables the reconstruction of previously unknown high-resolution variability in salinity, local rainfall and evaporation in past climates.
Combined with longer-term, global-scale paleoclimate records and models, our new method for absolute monthly SST and δ 18 Osw reconstructions has the potential to resolve critical disagreements between SST proxies, reduce biases of deep-time paleoclimate reconstructions, shed light on new aspects of past climate seasonality and reconcile proxy reconstructions and model simulations of greenhouse climate.
The three distinct localities contain a rich (> 200 species), well preserved Campanian rocky shore fauna 27,38 and were all deposited at the peak of transgression of the latest early Campanian, as supported by the restriction of these deposits to the Belemnellocamax mammillatus belemnite biozone and Sr-isotope chemostratigraphy 26,38,42 .The tectonic quiescence of the region since the Late Cretaceous limited burial and promoted excellent shell preservation 27 .Burial of loosely compacted sediments of the studied localities was limited to a maximum of 40 meters 27 .We can therefore conclude that burial temperatures never exceeded 80°C and that solid-state reordering did not affect clumped isotope results from these specimens 43,44 (see Supplementary Methods).
The Kristianstad Basin represents the highest latitude location for the occurrence of rudist bivalves known to date 38 .

Materials
Fossil R. diluvianum oysters were found in situ clinging to the sides of large boulders at a paleodepth of <5 m 27 , while the B. suecis rudist and A. incurva oysters were found in life position in a deeper setting (5-15m) among skeletal fragments on the paleo-seafloor 45 (see Supplementary Methods).The preservation of multiple specimens from this site (including the ones used here) was demonstrated through electron and visible light microscopy, trace element (e.g.Sr/Ca and Mn/Ca) analyses and ultrastructure preservation, results of which are reported in detail in 26 and 45 and summarized in Supplementary Methods.

Sampling
Powdered samples (~300 µg) were drilled in growth direction from polished cross sections through the shell's axis of maximum growth using a Dremel ® 3000 rotary drill (Robert Bosch Ltd., Uxbridge, UK) operated at slow rotation equipped with a 300 µm wide tungsten carbide drill bit.
High (~100 µm) uniform sampling resolution was achieved by careful abrasive drilling using the side of the drill parallel to the growth lines in the shell.In oyster shells, the well-preserved dense foliated calcite was targeted, while in the rudist the dense outer calcite was sampled, avoiding the honeycomb structure in the inner part of the outer shell layer which is more susceptible to diagenetic alteration 46 .A total of 145 samples was obtained, from which 338 aliquots of ~100 µg were analyzed (see Table 1).

Clumped isotope analyses
Clumped isotope (Δ47) analyses were carried out on Thermo Fisher Scientific MAT253 and 253 Plus mass spectrometers coupled to a Kiel IV carbonate preparation devices.Calcite samples (individual replicates of ~90 μg for MAT253 Plus and ~150 μg for MAT253) were reacted at 70 °C with nominally anhydrous (103%) phosphoric acid.The resulting CO2 gas was cleaned from water and organic compounds with two cryogenic LN2 traps and a PoraPak Q trap kept at −40 °C.
The purified sample gases were analyzed in micro-volume LIDI mode with 400 s integration time against a clean CO2 working gas (δ 13 C = -2.82‰;δ 18 O = -4.67‰),corrected for the pressure baseline 41,42 and converted into the absolute reference frame by creating an empirical transfer function from the daily analyzed ETH calcite standards (ETH-1, -2, -3) and their accepted values 49 .
All isotope data were calculated using the new IUPAC parameters following 49 and Δ47 values were projected to a 25 °C acid reaction temperature with a correction factor of 0.071 ‰ 22 .Longterm Δ47 reproducibility standard deviation was 0.04‰ (0.039‰ for MAT253 Plus and 0.045‰ for MAT253) based on repeated measurements of ~100 µg aliquots of our check standard IAEA C2 (Δ47 of 0.719‰; measured over a 20-month period; see Supplementary Data 8 for details).
No statistical difference was found between results from both instruments (see Supplementary Data 8).For the δ 18 Oc compositions we applied an acid correction factor of 1.00871 17 and reported versus VPDB with a typical reproducibility below 0.13‰ (95% confidence level).Results were combined with δ 18 Oc data previously measured in the same shells 26,45 (Supplementary Data 2) to improve the confidence of seasonal age models and the temporal resolution of SST and δ 18 Osw reconstructions.

Absolute paleoseasonality reconstructions
We reconstructed absolute SST seasonality by aligning Δ47 data relative to the seasonal cycle observed in δ 18 Oc using an age modelling routine 51 (Supplementary Data 1 and 9).Note that while chronologies were based on seasonal oscillations in δ 18 Oc records, the resulting age model is not compromised by unconstrained seasonal variability in δ 18 Osw (see discussion in [REF20] and Supplementary Methods).Since only the shape of the seasonally oscillations in δ 18 Oc is used for age modelling, age model results are independent on the absolute SST and δ 18 Osw seasonality and yield accurate results as long as the shape of the δ 18 Oc curve exhibits annual cyclicity (see 50 and 20 ; Supplementary Methods).We used a statistical approach to combine aliquots for Δ47based seasonality reconstructions.A step-by-step explanation of our Δ47-δ 18 Oc seasonality routine as well as a detailed evaluation of its precision and accuracy on a range of Δ47-δ 18 Oc datasets is provided in 20 and in Supplementary Methods.The number of 100 μg Δ47 aliquots to combine into monthly SST estimates is optimized by grouping aliquots from the same month in different growth years.Analytical uncertainties are propagated through this optimization procedure using Monte Carlo simulations (details in Supplementary Methods and Supplementary Data 10).
SST's are calculated from Δ47 values in monthly time bins (1/12 th of the seasonal cycle) using the temperature calibration by 51 recalculated in 23 , and δ 18 Osw is reconstructed from Δ47-SST and δ 18 Oc following 17 (Supplementary Methods and Supplementary Data 3).The accuracy of this statistical approach for combining Δ47 aliquots for seasonal SST and δ 18 Osw reconstructions is tested on a diverse group of modern datasets and evaluated in 20 .It is demonstrated that this method achieves the ideal compromise between eliminating bias and retaining high reproducibility while keeping SST and δ 18 Osw reconstructions independent of the δ 18 Oc values on which the age model is based 20 (see also Supplementary Methods).The clumped isotope temperature calibration by 23 is statistically indistinguishable from the temperature relationship based on theoretical principles within the temperature range discussed in 24 and is the culmination of recent convergence of measurement results between labs across the world and inter-lab standardization efforts 22,49 .
Seasonality is defined as the difference between the average temperatures in the warmest and coldest month, while mean annual temperature (MAT) is expressed as the average of all monthly temperatures, following USGS definitions 52 .Statistical analyses of seasonality, differences between specimens and differences between data and model are summarized in Supplementary Data 4.

Climate model
We utilize a fully equilibrated (>11,000 model years) paleoclimate model (HadCM3L) Campanian (78 Ma) simulation.Model boundary conditions (topography, bathymetry, solar luminosity) for the Campanian are the same as in 14 .We evaluate model simulations with radiative forcing (pCO2) set to 560 ppmV (2× preindustrial concentration) and 1120 ppmV (4× preindustrial concentration), within the range of pCO2 reconstructions for the Campanian as compiled by 53   The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries.This map has been provided by the authors.

Supplementary Files
This is a list of supplementary les associated with this preprint.Click to download. SupplementaryTextNatCommEarthEnv.pdf 15-27°C range, MAT of 20°C) represents the most accurate SST seasonality reconstructions for the Campanian Boreal Chalk Sea to date.Since shallow marine bio-archives can record local climate conditions at higher spatial and temporal resolution than conventional (open ocean) archives, our monthly resolved Δ47 records showcase a tool for eliminating widespread biases related to seasonal variability and δ 18 Osw assumptions on SST reconstructions across time and space.

Fig. 1 :Fig. 2 :
Fig. 1: Campanian (78 Ma) paleogeography and geological setting A) Global paleogeography used in climate model 14 B) Northern Europe, black star indicates the Kristianstad Basin C) Southern Sweden with Kristianstad Basin (in dark grey, submerged in the Campanian).Colored dots indicate the three sampled localities on the paleoshoreline with schematic representations of the species.In all maps, blue color indicates sea surface and light grey indicates emergent land.Maps B) and C) are adapted from 38 .

Fig. 3 :
Fig. 3: Comparison between model and reconstructions.A) Campanian latitudinal SST gradient with vertical orange, purple and green bars showing seasonality reconstructions and dashed black lines indicating modeled mean annual temperatures (560 ppmV = 2× preindustrial and 1120 ppmV = 4× preindstrial pCO2) with grey envelopes representing seasonality.Black symbols and bars show previous SST reconstructions 7,8,10,11,12,32,39,40 .The shaded yellow envelope indicates modern seasonal SST range 28 .Horizontal dashed lines mark modern and Campanian global MAT.B) Monthly SST reconstructions (in orange, purple and green) and local model outputs (in grey) in the Boreal Chalk Sea.Diamonds indicate monthly SST means, with red and blue diamonds showing monthly summer and winter extremes, respectively.Shaded envelopes show 95% confidence levels and color coding follows Fig. 2.

Figures Figure 1
Figures

Figure 3 Comparison
Figure 3

Table 1 : Overview of analytical results (δ 18 Oc and Δ47) and reconstructions Measurement results Monthly reconstructions
N = Number of measurements and the age is estimated from the age modelling results (Supplementary Data 1), CM = Coldest Month, MA = Mean Annual, WM = Warmest Month.All uncertainties are given as 95% confidence levels.
25,55 a modern astronomical configuration with dynamic vegetation.Details on the HadCM3L model are provided in Supplementary Methods and in14.Local seasonal SSTs are calculated for the paleorotated Kristianstad Basin 41 (42.5-50°N, 7.5-15°E; Supplementary Data 5) from averages of the upper ocean grid boxes in the model simulation.The model has a spatial resolution of 3.75° × 2.5° and uses 20 layers in ocean depth, of which the upper ocean box averages the top 10 meters of the ocean.Hence the average SST of the Kristianstad Basin is biased against the shallowest coastal regions of the basin, such as the locality of R. diluvianum54,55.For comparison, modern SST data come from the National Oceanic and Atmospheric Administration25(Supplementary Data 6 and Supplementary Methods).