Introduction

Seasonal extremes were of vital importance for the evolution and distribution of life over geological history1. The effects of greenhouse warming on seasonal variability in temperature and the hydrological cycle are still poorly constrained2,3, while being of considerable interest for projecting future climate and its impact on the ongoing biodiversity crisis4,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 CO2 levels7,8. Accurate reconstructions are imperative to evaluate climate model predictions under dissimilar climate states9,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 conditions11,12.

Reconstructions based on stable oxygen isotope ratios (δ18Oc) in marine carbonates and organic paleothermometry (e.g. TEX86) indicate that Late Cretaceous global mean sea surface temperatures (SST) were ~5–6 °C warmer than today11,12,13 with a reduced latitudinal temperature gradient (an “equable climate”14), while exhibiting limited temperature seasonality15,16. However, the reliability of many past seasonal reconstructions is undermined by seasonal bias and poorly constrained assumptions of seawater composition (δ18Osw bias17). This hampers our understanding of past warm climates and hinders accurate evaluation of climate models16,18,19.

Seasonal bias occurs if a proxy is interpreted as representing annual mean conditions but is in fact biased to a particular season17. Since fossil species producing the material that constitutes SST archives may not have a close modern relative for proxy calibration20, uncertainties about their growth seasons may unpredictably bias reconstructions21. Seawater oxygen isotope composition (δ18Osw) is an important input parameter into the widely used carbonate δ18Oc temperature proxy22, but it is highly variable across ocean basins (−3‰ to +2‰VSMOW23,24,25) and remains poorly constrained across geological timescales26,27. Biases in assumed δ18Osw composition thus undermine SST reconstructions, especially those from highly variable epicontinental seas28.

The advent of carbonate clumped isotope (Δ47) SST reconstructions on a seasonal scale promises to eliminate these two biases17,29. The clumped isotope thermometer yields accurate SST reconstructions independent of δ18Osw assumptions30. It also allows the reconstruction of δ18Osw, yielding information about the (local) hydrological cycle, an important aspect of climate rarely constrained in deep time, rectifying bias in the popular carbonate δ18Oc 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 disequilibrium31,32 (see Supplementary Discussion). The large sample sizes required for individual Δ47-based temperature estimates (>2 mg) have complicated paleoseasonality reconstructions using this accurate method33, but a recently developed statistical approach enables its use for seasonality reconstructions17.

Here we use clumped isotope analyses on microsampled (~100 µg) profiles through fossil bivalve shells to obtain, for the first time, absolute SST and δ18Osw 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 Ma34) coastal localities of the Kristianstad Basin in southern Sweden (46 ± 3°N paleolatitude35; 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.

Fig. 1: Campanian (78 million years ago) paleogeography and geological setting.
figure 1

A Global paleogeography used in climate model19 B Northern Europe, black star indicates the Kristianstad Basin C Southern Sweden with Kristianstad Basin (in dark gray, 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 gray indicates emergent land. Maps (B) and (C) are adapted from ref. 46.

Full size image

Results

All specimens showed clear seasonal δ18Oc 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 δ18Oc fluctuations reflect seasonality is demonstrated to be a valid basis for constructing intra-shell chronologies in nearly all modern environments17. Seasonal δ18Oc 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.

Table 1 Overview of analytical results (δ18Oc and Δ47) and reconstructions.
Full size table

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 δ18Osw 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 δ18Osw 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 δ18Osw 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 δ18Osw values (p < 0.01) than the other specimens. In all specimens, monthly δ18Osw positively correlates with monthly SST (see Fig. 2).

Fig. 2: Paleoseasonality reconstructions.
figure 2

A From bottom to top: relative monthly sampling frequencies reflecting growth rate variability (bar chart), monthly average clumped isotope value (Δ47), sea surface temperature (SST) and seawater oxygen isotope value (δ18Osw) reconstructions from R. diluvianum (orange), A. incurva (purple) and B. suecicus (green). Shaded envelopes indicate 95% confidence levels. Red and blue dots respectively indicate warmest and coldest months. B SST and δ18Osw reconstructions of warmest month (red symbols), coldest month (blue symbols) and annual average (symbols in color of specimen). Thin crosses indicate 95% confidence level uncertainties. Vertical bars on the right indicate summer, winter, and mean annual temperature (MAT) estimates from carbonate oxygen isotope values (δ18Oc; assuming constant δ18Osw of −1‰VSMOW40). Cross-sections through specimens are drawn with horizontal 10 mm scale bars.

Full size image

We compare reconstructed SST from this and previous studies with local and global Campanian SSTs modeled using the HadCM3BL-M2.1aE model36. 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 temperature36. Importantly for this work, it is sufficiently computationally efficient to allow the long simulations required to reach close to equilibrium for paleoclimates19 (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 pCO2 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 pCO2 conditions and resembles the present (8.2 °C) in the 4× preindustrial pCO2 simulation, corroborating recent studies arguing against the hypothesis of reduced seasonality during greenhouse conditions16,37. Campanian modeled MAT is ~18 °C and ~22 °C under 2× and 4× preindustrial atmospheric pCO2, respectively, compared to ~14 °C in the modern ocean38, yielding an equilibrium climate sensitivity, or global warming per doubling of atmospheric CO2 concentration, of ~4 °C19. 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 pCO2 forcing are significantly warmer than present day (3 ± 0.8–17 ± 0.4 °C38) and modeled pre-industrial local SST seasonality (−1.6 to +11.2 °C19).

Fig. 3: Comparison between model and reconstructions.
figure 3

A Campanian latitudinal sea surface temperature (SST) gradient with vertical orange, purple and green bars showing seasonality reconstructions from this study and dashed black lines indicating modeled mean annual temperatures (560 ppmV  = 2× preindustrial and 1120 ppmV = 4× preindstrial CO2 pressure) with gray envelopes representing seasonality. Black symbols and bars show previous SST reconstructions11,12,13,15,25,42,68,69. The large range in bulk mollusk clumped isotope (Δ47)-based SST estimates is discussed in Supplementary Discussion. The shaded yellow envelope indicates modern seasonal SST range (NOAA, 2021). Horizontal dashed lines mark modern and Campanian global mean annual temperature (MAT). B Monthly SST reconstructions (in orange, purple and green) and local model outputs (in gray) 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.

Full size image

Discussion

Comparison between specimens

Our novel Δ47-based monthly SST and δ18Osw 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 (Fig. 3). Higher (p < 0.05) SST (+4–5 °C) and δ18Osw (+1.0–1.5‰) in R. diluvianum are likely caused by local differences in its shallower, inter-tidal (<5 m) environment35. Temporary areal exposure during low tides could have elevated temperatures and δ18Osw recorded in R. diluvianum year-round by direct sunlight and evaporation, as in modern inter-tidal oyster species39. 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 δ18Osw signature (closer to the −1‰VSMOW assumed for ice-free oceans40), especially in winter. Such local environmental differences are not resolved in the climate model simulations but show the unprecedented detail of local SST and δ18Osw reconstructions from clumped isotope analyses in bivalve shells (see Supplementary Discussion). The ~1‰ δ18Osw seasonality shows that summers in the Campanian Kristianstad Basin either experienced excess evaporation, which increases δ18Osw by preferentially removing isotopically light seawater, or reduced precipitation, which supplies isotopically light meteoric water, reducing δ18Osw. Both processes lead to comparatively dry summers and wet winters.

Bias due to δ18Osw assumptions

Strong seasonal fluctuations in δ18Osw (up to 1.3‰ in B. suecicus) and regular deviations from the commonly assumed −1‰VSMOW δ18Osw value lead to large differences (up to 8.9 °C in R. diluvianum) between SST estimates based on Δ47 and δ18Oc (Fig. 2). The risk of assuming constant δ18Osw is even more clearly illustrated by significantly (+3.5–6.0 °C) higher δ18Oc-based seasonal temperature reconstructions for B. suecicus compared to A. incurva, while both specimens grew under similar SST seasonality conditions (Fig. 2B). Similarly, δ18Oc-based 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). This illustrates that the constant δ18Osw assumption is only valid in settings with negligible δ18Osw seasonality and where δ18Osw is known. Low-latitude Tethyan SST seasonality reconstructions based on rudist δ18Oc15 agree with model simulations, which may indicate that δ18Osw seasonality is less important in open marine settings. However, data–model agreement is by no means solid evidence for correct δ18Osw 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 δ18Oc-based SST reconstructions due to inaccurate δ18Osw assumptions and seasonal bias12,41. It must be noted that non-carbonate temperature proxies (e.g. based on organic chemistry or palaeobotanical evidence) rely on assumptions other than δ18Osw 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 δ18Osw seasonality between specimens would cause an unpredictable bias in MAT between −7.8 °C and +1.4 °C (Fig. 2A; Supplementary Data 3). Indeed, our Campanian mid-latitude SST ranges (~15–27 °C, MAT of 20 °C) are significantly higher than previous SST reconstructions of the same paleolatitude based on fish tooth δ18Oc (15–20 °C42), chalk δ18Oc (12–15 °C13), bulk mollusk Δ47 (5–12 °C25), TEX86 (15–20 °C12) and sub-annual mollusk δ18Oc (15–22 °C;34 Fig. 3A). On average, these previous studies yield lower MAT (~15 °C) than our seasonally controlled reconstructions and model (see Supplementary Discussion and Supplementary Fig. S4), but they highlight considerable variability between localities, even within the same study25 (Supplementary Discussion). The proxies used in these studies are affected differently by either seasonal or δ18Osw 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 diagenesis41,43.

Given the increase in frequency and duration of growth stops in modern mollusks with increasing latitude21, 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 δ18Osw 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 seasonality37, 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 pCO2 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 pCO2 simulation strongly favor higher (4× preindustrial pCO2) 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 (Supplementary Discussion). Bio-archives from mid to high latitudes are likely much more sensitive to δ18Osw and seasonality bias than low-latitude records21,25, contributing to the flawed paradigm of shallow latitudinal temperature gradients during greenhouse climates. Our results concur with the recent trend of converging data and model reconstructions yielding modern-scale Late Cretaceous latitudinal temperature gradients16,42, thereby challenging the hypothesis of “equable climate” during greenhouse periods14 (see Supplementary Discussion). Moreover, our unique absolute monthly SST reconstructions and model simulations corroborate growing evidence against the hypothesis of reduced temperature seasonality in greenhouse climates37. 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 and method with which to explore latitudinal seasonality gradients throughout the Mesozoic and across different ocean basins15.

Conclusions

Our new absolute temperature seasonality reconstructions merit critical evaluation of classical paleoclimate records that risk bias, such as those based on δ18Oc (assuming constant δ18Osw11,13), bulk analyses of fossil material with growth seasonality (e.g. mollusks and brachiopods25,27) or a fixed growth season (e.g. planktic foraminifera44) and organic proxies that may be seasonally biased (e.g. TEX86 and Uk’378,45). In addition, our monthly δ18Osw 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 precipitation, and evaporation in past climates. Combined with longer-term, global-scale paleoclimate records and models, our new method for absolute monthly SST and δ18Osw 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 Basin35,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 fauna35,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 chemostratigraphy35,46,48. The tectonic quiescence of the region since the Late Cretaceous limited burial and promoted excellent shell preservation34,35. Burial of loosely compacted sediments of the studied localities was limited to a maximum of 40 m35. We can therefore conclude that burial temperatures never exceeded 80 °C and that solid-state reordering did not affect clumped isotope results from these specimens49,50 (see Supplementary Methods). The Kristianstad Basin represents the highest latitude location for the occurrence of rudist bivalves known to date46.

Materials

Fossil R. diluvianum oysters were found in situ clinging to the sides of large boulders at a paleodepth of <5 m46, 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-seafloor51 (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 in34 and51 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 targeted52, 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 alteration53. 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 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 (δ13C = −2.82‰; δ18O = −4.67‰), corrected for the pressure baseline54,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 values31. We measured more ETH-3 standards to better constrain uncertainties around expected ∆47 values of our samples57. All isotope data were calculated using the new IUPAC parameters following58 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 δ18Oc compositions we applied an acid correction factor of 1.0087122 and reported versus VPDB with a typical reproducibility below 0.13‰ (95% confidence level). Results were combined with δ18Oc data previously measured in the same shells34,51 (Supplementary Data 2) to improve the confidence of seasonal age models and the temporal resolution of SST and δ18Osw reconstructions.

Absolute paleoseasonality reconstructions

We reconstructed absolute SST seasonality by aligning Δ47 data relative to the seasonal cycle observed in δ18Oc using an age modeling routine60 (Supplementary Data 1 and 9). Note that while chronologies were based on seasonal oscillations in δ18Oc records, the resulting age model is not compromised by unconstrained seasonal variability in δ18Osw (see discussion in ref. 17, Supplementary Methods and Supplementary Figs. S2–3). Since only the shape of the seasonal oscillations in δ18Oc is used for age modeling, age model results are independent from the absolute SST and δ18Osw seasonality and yield accurate results as long as the shape of the δ18Oc curve exhibits annual cyclicity17,60 (Supplementary Methods). We used a statistical approach to combine aliquots for Δ47-based seasonality reconstructions. A step-by-step explanation of our Δ47-δ18Oc seasonality routine as well as a detailed evaluation of its precision and accuracy on a range of Δ47-δ18Oc datasets is provided in17 and in Supplementary Methods. R functions used to calculate seasonality from the sample size optimization routine are compiled into the documented R package seasonalclumped and deposited in the open-source R code repository CRAN61. 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). SSTs are calculated from Δ47 values in monthly time bins (1/12th of the seasonal cycle) using the temperature calibration by62 recalculated in ref. 31, and δ18Osw is reconstructed from Δ47-SST and δ18Oc following22 (Supplementary Methods and Supplementary Data 3). Tests on a diverse group of modern datasets in ref. 17 demonstrate that this method achieves the ideal compromise between eliminating bias and retaining high reproducibility while keeping SST and δ18Osw reconstructions independent of the δ18Oc 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 efforts32,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 definitions64. 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 in19. 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 ref. 65, and a modern astronomical configuration with dynamic vegetation. HadCM3 is, to our knowledge, the only model to run Campanian-specific boundary conditions with a range of pCO2 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 conditions19, 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 temperature36. Details on the HadCM3L model are provided in Supplementary Methods and in ref. 19. Local seasonal SSTs are calculated for the paleorotated Kristianstad Basin47 (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. diluvianum66,67. For comparison, modern SST data come from the National Oceanic and Atmospheric Administration38 (Supplementary Data 6 and Supplementary Methods).