Absolute seasonal temperature estimates from clumped isotopes in bivalve shells suggest warm and variable greenhouse climate

Seasonal variability in sea surface temperatures plays a fundamental role in climate dynamics and species distribution. Seasonal bias can also severely compromise the accuracy of mean annual temperature reconstructions. It is therefore essential to better understand seasonal variability in climates of the past. Many reconstructions of climate in deep time neglect this issue and rely on controversial assumptions, such as estimates of sea water oxygen isotope composition. Here we present absolute seasonal temperature reconstructions based on clumped isotope measurements in bivalve shells which, critically, do not rely on these assumptions. We reconstruct highly precise monthly sea surface temperatures at around 50 °N latitude from individual oyster and rudist shells of the Campanian greenhouse period about 78 million years ago, when the seasonal range at 50 °N comprised 15 to 27 °C. In agreement with fully coupled climate model simulations, we find that greenhouse climates outside the tropics were warmer and more seasonal than previously thought. We conclude that seasonal bias and assumptions about seawater composition can distort temperature reconstructions and our understanding of past greenhouse climates. Mid-latitude sea surface temperatures during the Campanian greenhouse period 78 million years ago were warmer and more seasonal than thought, suggests a reconstruction at monthly timescales based on individual oyster and rudist shells.

S easonal 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 2,3 , while being of considerable interest for projecting future climate and its impact on the ongoing biodiversity crisis [4][5][6] . 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 CO 2 levels 7,8 . Accurate reconstructions are imperative to evaluate climate model predictions under dissimilar climate states 9,10 . Particularly seasonal range is poorly constrained with little quantitative evidence. The warm, ice free Late Cretaceous period presents a valuable reference period to assess seasonal variability under greenhouse conditions 11,12 .
Reconstructions based on stable oxygen isotope ratios (δ 18 O c ) in marine carbonates and organic paleothermometry (e.g. TEX 86 ) indicate that Late Cretaceous global mean sea surface temperatures (SST) were~5-6°C warmer than today [11][12][13] with a reduced latitudinal temperature gradient (an "equable climate" 14 ), while exhibiting limited temperature seasonality 15,16 . However, the reliability of many past seasonal reconstructions is undermined by seasonal bias and poorly constrained assumptions of seawater composition (δ 18 O sw bias 17 ). This hampers our understanding of past warm climates and hinders accurate evaluation of climate models 16,18,19 .
Seasonal bias occurs if a proxy is interpreted as representing annual mean conditions but is in fact biased to a particular season 17 . Since fossil species producing the material that constitutes SST archives may not have a close modern relative for proxy calibration 20 , uncertainties about their growth seasons may unpredictably bias reconstructions 21 . Seawater oxygen isotope composition (δ 18 O sw ) is an important input parameter into the widely used carbonate δ 18 O c temperature proxy 22 , but it is highly variable across ocean basins (−3‰ to +2‰VSMOW [23][24][25] ) and remains poorly constrained across geological timescales 26,27 . Biases in assumed δ 18 O sw composition thus undermine SST reconstructions, especially those from highly variable epicontinental seas 28 .
The advent of carbonate clumped isotope (Δ 47 ) SST reconstructions on a seasonal scale promises to eliminate these two biases 17,29 . The clumped isotope thermometer yields accurate SST reconstructions independent of δ 18 O sw assumptions 30 . It also allows the reconstruction of δ 18 O sw , yielding information about the (local) hydrological cycle, an important aspect of climate rarely constrained in deep time, rectifying bias in the popular carbonate δ 18 O c 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 with negligible influence of isotope disequilibrium 31,32 (see Supplementary Discussion). The large sample sizes required for individual Δ 47 -based temperature estimates (>2 mg) have complicated paleoseasonality reconstructions using this accurate method 33 , but a recently developed statistical approach enables its use for seasonality reconstructions 17 .
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 O sw 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.3 Ma 34 ) coastal localities of the Kristianstad Basin in southern Sweden (46 ± 3°N paleolatitude 35 ; see Fig. 1 and "Methods"). We compare these reconstructions with fully coupled climate model simulations of the Campanian greenhouse (see "Methods") to explore their implications for Late Cretaceous greenhouse climate.

Results
All specimens showed clear seasonal δ 18 O c fluctuations of −2.0-0.0‰ in R. diluvianum, −2.0-0.0‰ in A. incurva and −2.7‰ to −1.0‰ in B. suecicus on which shell chronologies were based (see "Methods"). The assumption that periodic δ 18 O c fluctuations reflect seasonality is demonstrated to be a valid basis for constructing intra-shell chronologies in nearly all modern environments 17 . Seasonal δ 18 O c patterns show that the specimens record 3 (A. incurva and B. suecicus) to 6 (R. diluvianum) full years of growth. Clumped isotope analyses on small aliquots yielded Δ 47 ranges between 0.62-0.73‰ for R. diluvianum, 0.64-0.76‰ for A. incurva, and 0.63-0.75‰ for B. suecicus. Summaries of measurement results are displayed in Table 1.
Detailed step-by-step results of the data processing routine are shown in ref. 17 , in Supplementary Methods and Supplementary Figs. S2-3. Figure 2 and Table 1 show monthly Δ 47 , SST and δ 18 O sw 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 variations in growth rate ( Fig. 2A) caused differences in the sample size between monthly time bins, combining data from the same month in multiple growth years allowed reliable SST and δ 18 O sw reconstructions for all monthly time bins in each specimen. Calculations of mean annual temperature (MAT) and seasonality from these monthly averages eliminate seasonal bias due to growth rate variability. Statistically significant (p < 0.01) SST seasonality was observed in all specimens. Summer and winter temperatures, defined as mean temperatures of the warmest and coldest month, in A. incurva (13 ± 2-26 ± 4°C) and B. suecicus (14 ± 4-25 ± 3°C) are statistically indistinguishable (p > 0.2), while SST from R. diluvianum are significantly higher (20 ± 2-29 ± 2°C; p < 0.05). Significant δ 18 O sw seasonality was found in R. diluvianum (0.0 ± 0.3-1.1 ± 0.3‰VSMOW; p < 0.01) and B. suecicus (−1.8 ± 0.8-0.6 ± 0.5‰ VSMOW; p < 0.01), but not in A. incurva (−0.9 ± 0.2 to −0.4 ± 0.9‰VSMOW; p = 0.07; Fig. 2; Table 1). R. diluvianum records significantly higher δ 18 O sw values (p < 0.01) than the other specimens. In all specimens, monthly δ 18 O sw positively correlates with monthly SST (see Fig. 2).
We compare reconstructed SST from this and previous studies with local and global Campanian SSTs modeled using the HadCM3BL-M2.1aE model 36 . Our model has been improved from being highly utilized in IPCC intercomparison assessment reports and compares well with CMIP5-generation model for many variables, including surface temperature 36 . Importantly for this work, it is sufficiently computationally efficient to allow the long simulations required to reach close to equilibrium for paleoclimates 19 (see "Methods"). 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 pCO 2 simulations (see "Methods"). Model results are summarized in Supplementary Data 5. The modeled Campanian latitudinal SST gradient (difference between tropics and high-latitude MAT; 26°C in both simulations) resembles the modern (25°C gradient). Modeled 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 pCO 2 conditions and resembles the present (8.2°C) in the 4× preindustrial pCO 2 simulation, corroborating recent studies arguing against the hypothesis of reduced seasonality during greenhouse conditions 16,37 . Campanian modeled MAT is~18°C and~22°C under 2× and 4× preindustrial atmospheric pCO 2 , respectively, compared to~14°C in the modern ocean 38 , yielding an equilibrium climate sensitivity, or global warming per doubling of atmospheric CO 2 concentration, of~4°C 19 . Specifically, simulated seasonal SST ranges in the Campanian Kristianstad Basin of 7 ± 3-20 ± 2°C for 2× and 12 ± 2-26 ± 2°C for 4× preindustrial atmospheric pCO 2 forcing are significantly warmer than present day (3 ± 0.8-17 ± 0.4°C 38 ) and modeled pre-industrial local SST seasonality (−1.6 to +11.2°C 19 ).

Discussion
Comparison between specimens. Our novel Δ 47 -based monthly SST and δ 18 O sw reconstructions from A. incurva and B. suecicus are statistically indistinguishable from 4× preindustrial pCO 2 simulations (p > 0.05) and significantly warmer than the 2× preindustrial pCO 2 simulations (>4°C higher MAT, p < 0.05) of local SST seasonality (Fig. 3). Higher (p < 0.05) SST (+4-5°C) and δ 18 O sw (+1.0-1.5‰) in R. diluvianum are likely caused by local differences in its shallower, inter-tidal (<5 m) environment 35 . Temporary areal exposure during low tides could have elevated temperatures and δ 18 O sw recorded in R. diluvianum year-round by direct sunlight and evaporation, as in modern inter-tidal oyster species 39 . By comparison, the deeper (5-15 m) subtidal environments of A. incurva and B. suecicus were unaffected by these processes and may have received more water with an open marine δ 18 O sw signature (closer to the −1‰VSMOW assumed for ice-free oceans 40 ), especially in winter. Such local environmental differences are not resolved in the climate model simulations but show the unprecedented detail of local SST and δ 18 O sw reconstructions from clumped isotope analyses in bivalve shells (see Supplementary Discussion). The~1‰ δ 18 O sw seasonality shows that summers in the Campanian Kristianstad Basin either experienced excess evaporation, which increases δ 18 O sw by preferentially removing isotopically light seawater, or reduced precipitation, which supplies isotopically light meteoric water, reducing δ 18 O sw . Both processes lead to comparatively dry summers and wet winters.
Bias due to δ 18 O sw assumptions. Strong seasonal fluctuations in δ 18 O sw (up to 1.3‰ in B. suecicus) and regular deviations from the commonly assumed −1‰VSMOW δ 18 O sw value lead to large differences (up to 8.9°C in R. diluvianum) between SST estimates based on Δ 47 and δ 18 O c (Fig. 2). The risk of assuming constant δ 18 O sw is even more clearly illustrated by significantly (+3.5-6.0°C ) higher δ 18 O c -based seasonal temperature reconstructions for B. suecicus compared to A. incurva, while both specimens grew under similar SST seasonality conditions (Fig. 2B) 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 O c -based SST reconstructions due to inaccurate δ 18 O sw assumptions and seasonal bias 12,41 . It must be noted that non-carbonate temperature proxies (e.g. based on organic chemistry or palaeobotanical evidence) rely on assumptions other than δ 18 O sw which may also bias reconstructions.
Seasonal bias. 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, considerable differences in growth rate and δ 18 O sw seasonality between specimens would cause an unpredictable bias in MAT between −7.8°C and +1.4°C  Supplementary Fig. S4), but they highlight considerable variability between localities, even within the same study 25 (Supplementary Discussion). The proxies used in these studies are affected differently by either seasonal or δ 18 O sw bias, or both. An additional source of cold bias affecting carbonate microfossil (most notably chalk) records is mixing of biogenic carbonate with carbonate cements precipitated under cooler sea floor conditions during early diagenesis 41,43 .
Given the increase in frequency and duration of growth stops in modern mollusks with increasing latitude 21 , seasonal biases are likely more common in higher latitude environments. 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 present a tool for eliminating widespread biases related to seasonal variability and δ 18 O sw assumptions on SST reconstructions across time and space by combining long-term MAT reconstructions with snapshots of climate on the seasonal scale. The average seasonal range reconstructed from our three specimens (15-27°C range, MAT of 20°C) likely represents the most accurate SST seasonality reconstructions for the Campanian Boreal Chalk Sea to date. The reconstructions are supported by the remarkable agreement between Δ 47 -based SST ranges and climate model simulations. Late Cretaceous reconstructions from the same latitude yield similar terrestrial MATs with higher seasonality 37 , analogous to modern terrestrial-marine contrast (see Supplementary Discussion), and corroborate our findings of warmer, highly seasonal Late Cretaceous climate.
Data-model comparison. Robust agreement between our reconstructions and the 4× preindustrial pCO 2 model simulation down to the monthly scale provides strong evidence for considerably (up to~8°C) warmer higher latitudes during the Late Cretaceous greenhouse compared to the present day. Significant disagreement between summer, winter, and annual SST reconstructions from every specimen in this study and the 2× preindustrial pCO 2 simulation strongly favor higher (4× preindustrial pCO 2 ) radiative forcing (see Supplementary Data 4). Point-by-point data-model comparisons show that most previous Late Cretaceous SST reconstructions from the same latitudes yield lower temperatures with lesser data-model agreement    11,13 ), bulk analyses of fossil material with growth seasonality (e.g. mollusks and brachiopods 25,27 ) or a fixed growth season (e.g. planktic foraminifera 44 ) and organic proxies that may be seasonally biased (e.g. TEX 86 and U k' 37 8,45 ). In addition, our monthly δ 18 O sw 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 Δ 47based seasonality reconstructions enables the reconstruction of previously unknown high-resolution variability in salinity, local precipitation, and evaporation in past climates. Combined with longer-term, global-scale paleoclimate records and models, our new method for absolute monthly SST and δ 18 O sw 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.

Methods
Geological setting. The bivalve specimens used in this study were obtained from the Ivö Klack (R. diluvianum), Åsen (A. incurva) and Maltesholm (B. suecicus) localities in the Kristianstad Basin 35,46 in southern Sweden (56°2' N, 14°9' E; 46 ± 3°N paleolatitude based on paleorotation by; 47 see Fig. 1 and Supplementary  Data 7). The three distinct localities contain a rich (>200 species), well-preserved Campanian rocky shore fauna 35,46 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 35,46,48 . The tectonic quiescence of the region since the Late Cretaceous limited burial and promoted excellent shell preservation 34,35 . Burial of loosely compacted sediments of the studied localities was limited to a maximum of 40 m 35 . 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 49,50 (see Supplementary Methods). The Kristianstad Basin represents the highest latitude location for the occurrence of rudist bivalves known to date 46 .
Materials. Fossil R. diluvianum oysters were found in situ clinging to the sides of large boulders at a paleodepth of <5 m 46 , while the B. suecis rudist and A. incurva oysters were found in life position in a deeper setting (5-15 m) among skeletal fragments on the paleo-seafloor 51 (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 34 and 51 and summarized in Supplementary Methods and Supplementary  Fig. S1.
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 52 , 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 53 . A total of 145 samples was obtained from which 338 aliquots of~100 µg was 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 CO 2 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 CO 2 working gas (δ 13 C = −2.82‰; δ 18 O = −4.67‰), corrected for the pressure baseline [54][55][56] 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 31 . We measured more ETH-3 standards to better constrain uncertainties around expected Δ 47 values of our samples 57 . All isotope data were calculated using the new IUPAC parameters following 58 and Δ 47 values were projected to a 25°C acid reaction temperature with a correction factor of 0.071 ‰ 59 . Long-term Δ 47 reproducibility standard deviation was 0.04‰ (0.039‰ for MAT253 Plus and 0.045‰ for MAT253; combined uncertainty of 0.077‰ at the 95% confidence level) 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 O c compositions we applied an acid correction factor of 1.00871 22 and reported versus VPDB with a typical reproducibility below 0.13‰ (95% confidence level). Results were combined with δ 18 O c data previously measured in the same shells 34,51 (Supplementary Data 2) to improve the confidence of seasonal age models and the temporal resolution of SST and δ 18 O sw reconstructions.
Absolute paleoseasonality reconstructions. We reconstructed absolute SST seasonality by aligning Δ 47 data relative to the seasonal cycle observed in δ 18 O c using an age modeling routine 60 ( Supplementary Data 1 and 9). Note that while chronologies were based on seasonal oscillations in δ 18 O c records, the resulting age model is not compromised by unconstrained seasonal variability in δ 18 O sw (see discussion in ref. 17 , Supplementary Methods and Supplementary Figs. S2-3). Since only the shape of the seasonal oscillations in δ 18 O c is used for age modeling, age model results are independent from the absolute SST and δ 18 O sw seasonality and yield accurate results as long as the shape of the δ 18 O c curve exhibits annual cyclicity 17,60 (Supplementary Methods). We used a statistical approach to combine aliquots for Δ 47 -based seasonality reconstructions. A step-by-step explanation of our Δ 47 -δ 18 O c seasonality routine as well as a detailed evaluation of its precision and accuracy on a range of Δ 47 -δ 18 O c datasets is provided in 17 17 demonstrate that this method achieves the ideal compromise between eliminating bias and retaining high reproducibility while keeping SST and δ 18 O sw reconstructions independent of the δ 18 O c values on which the age model is based (see also Supplementary Methods). The clumped isotope temperature calibration by ref. 31 is statistically indistinguishable from the temperature relationship based on theoretical principles within the temperature range discussed in ref. 32 and is the culmination of recent convergence of measurement results between labs across the world and inter-lab standardization efforts 32,59,63 . 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 64 . 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 (HadCM3BL-M2.1aE) Campanian (78 Ma) simulation. Model boundary conditions (topography, bathymetry, solar luminosity) for the Campanian are the same as in 19 . We evaluate model simulations with radiative forcing (pCO 2 ) set to 560 ppmV (2× preindustrial concentration) and 1120 ppmV (4× preindustrial concentration), within the range of pCO 2 reconstructions for the Campanian as compiled by ref. 65 , and a modern astronomical configuration with dynamic vegetation. HadCM3 is, to our knowledge, the only model to run Campanianspecific boundary conditions with a range of pCO 2 forcing out to full equilibrium. This is critical as it was shown that the deep ocean and hence ocean circulation can take at least 5000 modeled years to fully equilibrate to the applied boundary conditions 19 , casting doubt on the validity of alternative model simulations that did not attain equilibrium. The model also compares well with CMIP5-generation model for many variables, including surface temperature 36 . Details on the HadCM3L model are provided in Supplementary Methods and in ref. 19 . Local seasonal SSTs are calculated for the paleorotated Kristianstad Basin 47 (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 m 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. diluvianum 66