Shell chemistry of the boreal Campanian bivalve Rastellum diluvianum (Linnaeus, 1767) reveals temperature seasonality, growth rates and life cycle of an extinct Cretaceous oyster

The Campanian age (Late Cretaceous) is characterized by a warm greenhouse climate with limited land-ice volume. This makes this period an ideal target for studying climate dynamics during greenhouse periods, which are essential for predictions of future climate change due to anthropogenic greenhouse gas emissions. Well-preserved fossil shells from the Campanian (±78 Ma) high mid-latitude (50^∘ N) coastal faunas of the Kristianstad Basin (southern Sweden) offer a unique snapshot of short-term climate and environmental variability, which complements existing long-term climate reconstructions. In this study, we apply a combination of high-resolution spatially resolved trace element analyses (micro-X-ray fluorescence – µXRF – and laser ablation inductively coupled plasma mass spectrometry – LA-ICP-MS), stable isotope analyses (IRMS) and growth modeling to study short-term (seasonal) variations recorded in the oyster species Rastellum diluvianum from the Ivö Klack locality. Geochemical records through 12 specimens shed light on the influence of specimen-specific and ontogenetic effects on the expression of seasonal variations in shell chemistry and allow disentangling vital effects from environmental influences in an effort to refine paleoseasonality reconstructions of Late Cretaceous greenhouse climates. Growth models based on stable oxygen isotope records yield information on the mode of life, circadian rhythm and reproductive cycle of these extinct oysters. This multi-proxy study reveals that mean annual temperatures in the Campanian higher mid-latitudes were 17 to 19 ^∘C, with winter minima of ∼13 ^∘C and summer maxima of 26 ^∘C, assuming a Late Cretaceous seawater oxygen isotope composition of −1 ‰ VSMOW (Vienna standard mean ocean water). These results yield smaller latitudinal differences in temperature seasonality in the Campanian compared to today. Latitudinal temperature gradients were similar to the present, contrasting with previous notions of “equable climate” during the Late Cretaceous. Our results also demonstrate that species-specific differences and uncertainties in the composition of Late Cretaceous seawater prevent trace element proxies (Mg∕Ca, Sr∕Ca, Mg∕Li and Sr∕Li) from being used as reliable temperature proxies for fossil oyster shells. However, trace element profiles can serve as a quick tool for diagenesis screening and investigating seasonal growth patterns in ancient shells.

1 Introduction

The Late Cretaceous was marked by a long cooling trend that brought global mean annual temperatures (MATs) down from the mid-Cretaceous climate maximum (±28 ^∘C surface ocean temperatures) in the Cenomanian and Turonian (±95 Ma) to slightly cooler temperatures (±22 ^∘C surface ocean temperatures) around the Campanian–Maastrichtian boundary (±72.1 Ma; Clarke and Jenkyns, 1999; Pearson et al., 2001; Huber et al., 2002; Friedrich et al., 2012; Scotese, 2016). This cooling trend was likely caused by a change in ocean circulation, initiated by the opening of the Equatorial Atlantic Gateway that separated the proto-North and proto-South Atlantic Ocean basins (Friedrich et al., 2009). The cooling trend is well documented in stable oxygen isotope records from the white-chalk successions of the Chalk Sea, which covered large portions of northwestern Europe during the Late Cretaceous period (Reid, 1973; Jenkyns et al., 1994; Jarvis et al., 2002; Voigt et al., 2010). The connection of the Chalk Sea to the (proto-)North Atlantic Ocean makes it an interesting area of study to constrain Late Cretaceous paleogeography and climate. These chalk successions are featured in various paleoclimate studies because they are accessible in good outcrops and consist predominantly of calcareous nannofossils which were thought to faithfully record sea surface conditions (e.g., Jenkyns et al., 1994). However, recent studies (e.g., O'Brien et al., 2017; Tagliavento et al., 2019) have shown that diagenetic overprinting likely biases these records towards cooler temperatures, resulting in the apparent cool tropics paradox (Pearson et al., 2001). Even though sea level changes seem to indicate possible small changes in land-ice volume during the Late Cretaceous, warm high-latitude paleotemperatures likely rule out the possibility of extensive polar ice sheets comparable in volume to modern ice caps (Barrera and Johnson, 1999; Huber et al., 2002; Jenkyns et al., 2004; Miller et al., 2005; Thibault et al., 2016).

Given its warm, essentially land-ice-free conditions, the Campanian serves as an interesting analogue for Earth's climate in the future, should anthropogenic and natural emissions continue to contribute to the rise in global temperatures and decrease in global ice volume on Earth (IPCC, 2013; Donnadieu et al., 2016). However, most Late Cretaceous climate reconstructions focus on reconstructing and modeling long-term evolutions of humid or arid conditions on land or past atmospheric and oceanic temperatures (DeConto et al., 1999; Thibault et al., 2016; Yang et al., 2018). Data on the extent of seasonal variability from this time period, especially outside the tropics, are scarce, although seasonality constitutes a fundamental component of the climate system (Steuber, 1999; Steuber et al., 2005; Burgener et al., 2018). Furthermore, many proxies used for paleoclimate reconstruction risk seasonal bias, and thus independent seasonality reconstructions serve as important tools to verify other climate reconstructions.

2 Background

2.1 Geological setting

The Kristianstad Basin, our study site, is located on the southeastern Baltic Sea coast of the southern Swedish province of Scania, which is presently located at 56^∘2^′ N, 14^∘9^′ E (see Fig. 1), somewhat higher than its Campanian paleolatitude, which is estimated at 50^∘ N (van Hinsbergen et al., 2015). The uppermost lower Campanian shallow-marine sediments deposited at Ivö Klack consist of sandy and silty nearshore deposits containing carbonate gravel (Christensen, 1975, 1984; Surlyk and Sørensen, 2010; Sørensen et al., 2015). The depositional setting is described as a rocky coastline that was inundated during the maximum extent of the Late Cretaceous transgression (Kominz et al., 2008; Csiki-Sava et al., 2015). Since the region has remained tectonically quiet since the Campanian, the deposits of Kristianstad Basin localities remain at roughly the same altitude as when they were deposited and have been subject to limited burial (Christensen, 1984; Surlyk and Sørensen, 2010).

The rocky-shore deposits of Ivö Klack are characterized by a diverse shelly fauna, consisting of well-preserved fossils and fragments of brachiopods, belemnites, echinoids and asteroids, polychaete worms, gastropods, corals, ammonites, and thick-shelled oysters, with a total of almost 200 different recognized species (Surlyk and Sørensen, 2010). In this diverse rocky-shore ecosystem, various habitat zones can be distinguished, each with their distinct suite of organisms adapted to local conditions of varying amounts of sunlight, sedimentation and turbulence (Surlyk and Christensen, 1974; Sørensen et al., 2012). This unique combination of marine biodiversity and preservation of original calcite shell material makes the localities in the Kristianstad Basin ideal for studying sub-annual variability in shell chemistry and reconstructing paleoseasonality and environmental change in the Campanian (Sørensen et al., 2015).

2.2 Bivalve shells

Fossil bivalve shells offer a valuable archive for studying past climates on a seasonal scale. Incremental measurements of shell chemistry in the growth direction (sclerochronology studies) potentially yield records of seasonal changes in the environment (Mook, 1971; Jones, 1983; Klein et al., 1996a; Schöne and Gillikin, 2013). Their distribution allows paleoseasonality reconstructions from bivalves across a wide range of latitudes (Roy et al., 2000; Jablonski et al., 2017), and the preservation potential of calcitic shell structures (especially in oyster shells) makes them ideal, if not unique, recorders of pre-Quaternary seasonality and sub-annual environmental change (Brand and Veizer, 1980; 1981; Al-Aasm and Veizer, 1986a, b; Immenhauser et al., 2005; Alberti et al., 2017). The incremental deposition of shell carbonate means, in theory, that the limits in terms of time resolution of the bivalve shell archive are governed by sampling resolution rather than the resolution of the record itself. While periods of growth cessation can occur (especially in high latitudes; Ullmann et al., 2010) and the true mechanisms of shell deposition at a very high (e.g., daily) temporal resolution are poorly constrained (see de Winter et al., 2020a, and references therein), in practice incremental shell deposition allows reconstructions of changes down to sub-daily timescales given the right sampling techniques (Schöne et al., 2005; Sano et al., 2012; Warter et al., 2018; de Winter et al., 2020a). Examples of chemical proxies used in such sclerochronology studies include stable carbon and oxygen isotope ratios and trace element ratios (e.g., Steuber et al., 2005; Gillikin et al., 2006; McConnaughey and Gillikin, 2008; Schöne et al., 2011; de Winter et al., 2017a, 2018).

2.3 Trace element proxies

The incorporation of trace elements and carbon isotopes into bivalve shells is influenced by so-called vital effects: biological controls on the incorporation of elements in the shell independent of the environment (Weiner and Dove, 2003; Gillikin et al., 2005). These vital effects have been shown to mask the characteristic relationships between shell trace element chemistry and the environment and appear to be distinct not only between different bivalve species but also between specimens of different ontogenetic age (Freitas et al., 2008). Differences between bivalve families mean that the trace element chemistry of some taxa (like scallops: family Pectinidae) are especially affected by vital effects (Lorrain et al., 2005; Freitas et al., 2008), while other families like oysters (family Ostreidae) seem to be more robust recorders of environmental conditions (Surge et al., 2001; Surge and Lohmann, 2008; Ullmann et al., 2010, 2013). The effect on shell composition and preservation of changes in microstructure and the amount of organic matrix present in different shell layers introduce uncertainty as to which parts of the shells are well-suited for reconstruction purposes (Carriker et al., 1991; Kawaguchi et al., 1993; Dalbeck et al., 2006; Schöne et al., 2010, 2013). The key to disentangling vital effects from recorded environmental changes lies in the application of multiple proxies and techniques on the same bivalve shells (the “multi-proxy approach”; e.g., Ullmann et al., 2013; de Winter et al., 2017a, 2018) and basing reconstructions on more than one shell (Ivany, 2012).

2.4 Stable isotope ratios

Because nearly all bivalves precipitate their shells at or near oxygen isotope equilibrium, the stable oxygen isotope ratio of bivalve shell carbonate is less susceptible to vital effects, such as growth kinetics (Uchikawa and Zeebe, 2012). Therefore, stable oxygen isotope ratios in bivalve shell carbonate (δ¹⁸O_c) are solely dependent on calcification temperature and seawater oxygen isotope composition (δ¹⁸O_sw), and this proxy is frequently used in sclerochronology studies as a paleothermometer (Kim and O'Neil, 1997; Schöne et al., 2005; Butler et al., 2013; Ullmann et al., 2013; Huyghe et al., 2015; de Winter et al., 2020b). Oxygen isotope records can function as a reference in the above-mentioned multi-proxy studies aimed at resolving vital effects and environmental and climatic changes. However, the weakness of this proxy lies in the fact that δ¹⁸O_sw is not always known, especially in deep-time settings (Veizer and Prokoph, 2015). As a result, a constant δ¹⁸O_sw of 0 ‰ VSMOW (Vienna standard mean ocean water) for modern icehouse climate conditions, or −1 ‰ VSMOW for an ice-free world (such as the Eocene or the Late Cretaceous; after Shackleton, 1986), is often assumed (e.g., Andreasson and Schmitz, 1996; Ivany and Runnegar, 2010; Huyghe et al., 2015). An offset of 1 ‰ between assumed δ¹⁸O_sw and actual δ¹⁸O_sw can result in a ∼4.6 ^∘C temperature offset in temperature reconstructions (Kim and O'Neil, 1997). This assumption may therefore introduce inaccuracies in absolute temperature reconstructions, but relative variations in δ¹⁸O_c can still yield important insights into high-resolution climate dynamics.

In marine mollusks, dissolved inorganic carbon (DIC) in the ambient seawater contributes the largest fraction of carbon (90  %) used for shell mineralization (McConnaughey, 2003; Gillikin et al., 2007) and therefore heavily influences δ¹³C values of shell carbonate. However, changes in respiration rates can alter the carbon budget of shell carbonate by adding or removing isotopically light respired carbon in the form of CO₂ (Lorrain et al., 2004). Environmental changes in DIC can also have a strong influence on this carbon budget, especially when bivalves grow in nearshore or estuarine conditions with large (seasonal) variations in environmental δ¹³C of DIC and organic carbon (Gillikin et al., 2006). Conceptual models exist that attempt to correlate shell δ¹³C in modern mollusks to environmental and physiological variations, but these require knowledge of ambient CO₂ pressures and δ¹³C values of DIC, gas ventilation rates in the animal, and CO₂ and O₂ permeabilities of calcifying membranes (McConnaughey et al., 1997), which are not available in fossil bivalve studies.

2.5 Aim

In this study, we present a detailed, multi-proxy comparison of the growth and chemistry of well-preserved fossil shells of the thick-shelled oyster Rastellum diluvianum (Linnaeus, 1767) recovered from the Ivö Klack locality on the northern edge of the Kristianstad Basin. We combine stable isotope proxies, conventional in sclerochronological studies (δ¹³C and δ¹⁸O; e.g., Goodwin et al., 2001; Steuber et al., 2005), with less well established trace element proxies (Mg∕Ca, Sr∕Ca, Mg∕Li and Sr∕Li; e.g., Bryan and Marchitto, 2008; Schöne et al., 2011; Füllenbach et al., 2015; Dellinger et al., 2018) and growth modeling based on δ¹⁸O seasonality (Judd et al., 2018) in an attempt to disentangle the effects of growth rate, reproductive cycle and environmental change on shell chemistry. The data gathered in this study allow a detailed discussion on seasonal changes in temperature and water chemistry in the coastal waters of the Kristianstad Basin in the late early Campanian as well as on the life cycle of R. diluvianum and its response to seasonal changes in its environment.

Figure 1 Paleogeographic map of the Boreal Chalk Sea (a) and the area of present-day southern Sweden (b) showing the location of Ivö Klack (modified after Sørensen et al., 2015).

3 Materials and methods

3.1 Sample acquisition and preparation

Complete valves of 12 individual R. diluvianum oysters were obtained from the Ivö Klack locality (see Fig. 2). Specimens of R. diluvianum were found in situ attached to the vertical sides of large boulders that characterized the rocky shore of Ivö Klack (Surlyk and Christensen, 1974) and are biostratigraphically assigned to the latest early Campanian Belemnellocamax mammilatus belemnite zone. The valves were cleaned and fully embedded in Araldite®2020 epoxy resin (Bodo Möller Chemie Benelux, Antwerp, Belgium). Dorsoventral slabs (±10 mm thick) were cut perpendicular to the hinge line using a water-cooled slow-rotating saw with a diamond-coated blade (thickness ±1 mm; Fig. 2). The surfaces cut on the central growth axis were progressively polished using silicon-carbide polishing disks (up to P2500, or 8.4 µm grain size). Polished surfaces were scanned at high (6400 dpi) resolution using an Epson Perfection 1650 flatbed color scanner (Seiko Epson Corp., Suwa, Japan). Resulting color scans of all polished R. diluvianum shell cross sections are provided in Fig. 2 and Fig. S1 from de Winter et al. (2020c, https://doi.org/10.5281/zenodo.2581304). Shell microstructures in R. diluvianum shells were studied in detail on high-resolution scans and by using reflected-light optical microscopy. Microstructural features were used to reconstruct the relative timing of shell growth (see Fig. 3). Fragments of visually well-preserved material from different microstructures in the shells were coated with gold and studied under a scanning electron microscope (SEM; Quanta 200 ESEM) and imaged at 1000×–2000× magnification (Fig. 3b–e). Chemical analyses were carried out sequentially on polished cross sections in order of destructive character of sampling (starting with the least-destructive measurements: micro-X-ray fluorescence – μXRF, laser ablation inductively coupled plasma mass spectrometry – LA-ICP-MS, microdrilling for IRMS and finally multi-cup inductively coupled plasma mass spectrometry – MC-ICP-MS – analysis on ∼26 mg samples).

Figure 2 Overview of the 12 Rastellum diluvianum shells used in this study. All shells are depicted at the same scale (see scale bar in center of image). Colors of the lines under sample names correspond to the colors of the lines in Figs. 4, 6, 8 and 10. Every shell is represented by an image of the inside of the valve analyzed as well as a color scan of the cross section through the shell on which high-resolution analyses were carried out. Dashed red lines show the location of these cross sections. The largest five shells (1–5, in top panel) were sampled for IRMS analyses (δ¹³C and δ¹⁸O). All shells were subjected to micro-X-ray fluorescence (µXRF), laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) and multi-cup inductively coupled plasma mass spectrometry (MC-ICP-MS) analyses. Full-size versions of the high-resolution color scans of shell cross sections are provided in de Winter et al. (2020c).

Figure 3 Overview image showing a high-resolution color scan of the cross section through R. diluvianum (a) on which the different shell textures as well as the directions of high-resolution analyses (in growth direction) are indicated. Thin blue lines denote parts of the shell that were deposited at the same time (growth increments). Areas marked in yellow represent sediment infilling below and within the valves. (b) and (c) show SEM images of the well-preserved foliated calcite in the shell. More porous structures in the shell (vesicular calcite) are depicted in SEM images shown in (d) and (e). Below are three XRF elemental maps of the same cross section: an RGB-colored map displaying the relative abundances of Fe, Si and Ca (f); a heatmap of Fe concentrations (g; see scale bar below map) and a heatmap of Sr concentrations (h; see scale bar below map). XRF mapping only yields relative (semi-quantitative) abundance of elements.

3.2 Micro-XRF mapping

Elemental abundance maps of all R. diluvianum shell cross sections were obtained using a Bruker Tornado M4 energy-dispersive micro-X-ray fluorescence scanner (µXRF; Bruker Nano GmbH, Berlin, Germany). All µXRF analyses carried out with the Bruker M4 Tornado are non-destructive. µXRF is equipped with a Rh filament metal–ceramic-tube X-ray source operated at 50 kV and 600 µA (30 W; maximum energy settings). The circular spot projected on the same surface is estimated to have a diameter of 25 µm (Mo Kα). A micrometer-precision XYZ translation stage allows for quick and precise sample movement such that a grid of 25 µm XRF spots can be measured on the sample surface by continuous scanning to construct elemental maps (3×10⁶–5×10⁶ pixels per map). Exposure times of the X-ray beam per sampling position in mapping mode (1 ms per pixel) are too short to gain an adequate signal-to-noise ratio for pixel-by-pixel quantification of elemental concentrations. Instead, processing of entire map surfaces using the Bruker ESPRIT™ software allows semi-quantitative elemental abundance maps to be created of the sample surface based on a mapping of the count rate in “Regions of Interest” of elements (see de Winter and Claeys, 2016; de Winter et al., 2017b; Fig. 3). XRF maps allow for a rapid assessment of the preservation state of original shell calcite based on variations in Si, Mn, Fe and Sr concentrations and guide the formulation of sampling protocols for further analyses (de Winter and Claeys, 2016; Fig. 3). Results of XRF mapping on all 12 R. diluvianum shell cross sections are provided in Fig. S2 from de Winter et al. (2020c).

3.3 Micro-XRF line scans

After XRF mapping, quantitative line scans were measured in the growth direction on shell cross sections. Dwell times of 60 s per measurement allowed sufficiently high signal-to-noise ratios for individual points in line scans to be quantified. This acquisition time was chosen so as to provide the optimal compromise between increasing run time (improved signal-to-noise ratio; enhanced reproducibility) and increasing the number of sampling positions (improving sampling density and allowing duplicate measurements) for the elements Mg, Al, Si, P, S, Ca, Ti, Mn, Fe, Cu, Zn and Sr (see discussion in de Winter et al., 2017b). The sampling interval of line scans was 50 µm, adding up to a total of 11 056 individual quantitative XRF spectra measured for this study. Spectra were quantified using the Bruker ESPRIT software calibrated using the matrix-matched BAS-CRM393 limestone standard (Bureau of Analysed Samples, Middlesbrough, UK), after which individual measurements were calibrated offline using seven matrix-matched certified reference materials (CCH1, COQ1, CRM393, CRM512, CRM513, ECRM782 and SRM1d), which were treated as samples (see Vansteenberge et al., 2020). R² values of calibration curves exceeded 0.99, and reproducibility standard deviations were better than 10  % relative to the mean. Even though line scans were positioned on well-preserved shell calcite based on the XRF map results, a second check was carried out in which individual points were rejected based on conservative thresholds for diagenetic recrystallization or detrital contamination ([Ca] <38 wt  %, [Si] >1 wt  %, [Mn] >200 µg g⁻¹ or [Fe] >250 µg g⁻¹; [Sr]∕[Mn] <100 mol mol⁻¹; see Al-Aasm and Veizer, 1986a; Sørensen et al., 2015). Concentrations of Ca, Mg and Sr in well-preserved shell sections were used to explore the potential of Mg∕Ca and Sr∕Ca molar ratios as paleoenvironmental proxies. Unprocessed results of XRF line scanning are provided in Fig. S3 from de Winter et al. (2020c).

3.4 LA-ICP-MS line scans

Spatially resolved elemental concentrations for Li, B, Mg, Si, P, Ca, Ti, V, Cr, Mn, Fe, Ni, Zn, Rb, Sr, Ba, Pb and U were calculated from a calibrated transient MS signal recorded during LA-ICP-MS line scanning in the growth direction (parallel to the XRF line scans) on shell cross sections. LA-ICP-MS measurements were carried out at the Atomic & Mass Spectrometry – A&MS – research unit of Ghent University (Ghent, Belgium) using a 193 nm ArF*excimer-based Analyte G2 laser ablation system (Teledyne Photon Machines, Bozeman, USA), equipped with a HelEx 2 double-volume ablation cell, coupled to an Agilent 7900 quadrupole-based ICP-MS unit (Agilent, Tokyo, Japan). Continuous scanning along shell transects using a laser spot with a diameter of 25 µm, scan speed of 50 µm s⁻¹ and detector mass sweep time of 0.5 s yielded profiles with a lateral sampling interval of 25 µm, amounting to a total of 9505 LA-ICP-MS data points gathered. The aerosol was transported using He carrier gas into the ICP-MS unit via the aerosol rapid introduction system (ARIS; Teledyne Photon Machines, Bozeman, USA). Elemental concentrations were calibrated using bracketed analysis runs on US Geological Survey (USGS) BCR-2G, BHVO-2G, BIR-1G, GSD-1G and GSE-1G and National Institute of Standards and Technology (NIST) SRM612 and SRM610 certified reference materials. Calcium concentrations (measured via ⁴³Ca) were used as an internal standard for data normalization and drift correction during the measurement campaign, and Ca concentrations of 38.5 wt  % were assumed for pristine shell carbonate. Coefficients of determination (R²) of a linear model fitted to the calibration curves were better than 0.99, and the standard deviation of reproducibility for elemental concentrations was better than 5 % relative to the mean value. Individual LA-ICP-MS measurements were inspected for diagenetic alteration or contamination by detrital material using the same thresholds as those used for XRF data (see above). LA-ICP-MS and µXRF measurements were combined to cover a wider range of elements, since some elements (e.g., S and Sr) were measured more reliably using µXRF, while others (e.g., Li or Ba) could only be quantified using LA-ICP-MS. Concentrations of Li, Mg and Sr were used to explore the potential of Mg∕Li and Sr∕Li molar ratios as proxies for paleoenvironmental change. Unprocessed results of LA-ICP-MS line scans are provided in Fig. S4 from de Winter et al. (2020c).

3.5 Isotope ratio mass spectrometry

A transect of powdered samples (±200 µg) was extracted for isotope ratio mass spectrometry (IRMS) analysis in the growth direction along well-preserved foliated calcite (Fig. 3) in the 5 largest of the 12 R. diluvianum shells (R. diluvianum 1–5; see Fig. 2) using a microdrill (Merchantek/Electro Scientific Industries Inc., Portland, OR, USA) equipped with a 300 µm diameter tungsten carbide drill bit, coupled to a microscope (Leica GZ6, Leica Microsystems GmbH, Wetzlar, Germany). A total of 531 IRMS samples were taken at an interspacing of 250 µm. Stable carbon and oxygen isotope ratios (δ¹³C and δ¹⁸O) were measured in a Nu Perspective IRMS equipped with a NuCarb carbonate preparation device (Nu Instruments, UK). The sample size (50–100 µg) allowed duplicate measurements to be carried out regularly (roughly once every 30 samples) to assess reproducibility. Samples were digested in 104 % phosphoric acid at a constant temperature of 70 ^∘C, and the resulting CO₂ gas was cryogenically purified before being led into the IRMS through a dual-inlet system. Isotope ratios were corrected for instrumental drift and fractionation due to variations in sample size, and the resulting values are reported relative to the Vienna Pee Dee Belemnite standard (‰ VPDB) using repeated measurements of the IA-603 stable isotope standard (International Atomic Energy Agency, Vienna, Austria). Reproducibility of δ¹⁸O and δ¹³C measurements on this standard were better than 0.1 ‰ and 0.05 ‰ (1σ; N=125), respectively. All stable isotope analysis results are provided in Fig. S5 from de Winter et al. (2020c), and plots of stable isotope and trace element records from all shells are given in Fig. S6 from de Winter et al. (2020c).

3.6 Growth and age modeling

Stable oxygen isotope curves measured in R. diluvianum were used to produce age models for the growth of the shell using a bivalve growth model written in MATLAB (MathWorks, Natick, MA, USA) which simulates δ¹⁸O curves using a combination of growth and temperature sinusoids (Judd et al., 2018). This simulation model was modified to calculate its temperatures based on calcite δ¹⁸O (following Kim and O'Neil, 1997) rather than from the aragonite δ¹⁸O–temperature relationship used in the original approach (after Grossman and Ku, 1986; see Judd et al., 2018). A value of −1.0 ‰ VSMOW was assumed for δ¹⁸O of Campanian ocean water (Shackleton, 1986; Thibault et al., 2016). Additional minor modifications in the source code allowed results of intermediate calculation steps in the model to be exported. The modified MATLAB source code is provided in Fig. S7 from de Winter et al. (2020c). Note that this model assumes that the shape and absolute value of δ¹⁸O curves depend solely on water temperature and growth rate (ignoring changes in δ¹⁸O_sw) and that a modeled year contains 365 d by construction (while this number should be slightly larger in the Late Cretaceous; de Winter et al., 2020a). Despite these caveats, shell chronologies reconstructed from seasonal patterns in δ¹⁸O should still be reliable, as they are only based on the shape of the δ¹⁸O curves. Uncertainties of modeled temperature curves were derived by propagating the measurement uncertainty of δ¹⁸O. Age models thus obtained for shells R. diluvianum 1–5 were used to align all proxy records on a common time axis. Age models for R.diluvianum 6–12 were constructed by extrapolating relationships between modeled seasonality and microstructures and trace element concentrations observed in R. diluvianum 1–5. Simultaneously deposited microstructural features in shell cross sections (see Fig. 3) were used to determine the actual dorsoventral height of the shells at different ages, linking shell height to the age and allowing the construction of growth curves for all 12 R. diluvianum shells. The total age and the season of spawning (or the start of shell growth) were determined by extrapolating the δ¹⁸O-based age models and by using the relationship between δ¹⁸O profiles and trace element records and growth increments observed in the shell.

3.7 Strontium isotope analysis

Samples (on average 26 mg) for strontium isotope analysis were obtained by drilling the well-preserved foliated calcite in all shells using a Dremel 3000 dental drill with a 0.5 mm tungsten carbide drill bit. Calcite samples were placed in Teflon beakers (Savillex LLC, Eden Prairie, MN, USA), dissolved in subboiled concentrated (14 M) nitric acid (HNO₃) at 120 ^∘C and left to dry at 90 ^∘C overnight, after which the residue was redissolved in 2 M HNO₃. Carbonate-bound strontium in the samples was purified following the ion-exchange resin chromatography method detailed in Snoeck et al. (2015). The ⁸⁷Sr∕⁸⁶Sr of purified Sr samples was determined using a Nu Plasma (Nu Instruments Ltd, Wrexham, UK) multi-collector (MC) ICP-MS unit in operation at the Université Libre de Bruxelles (ULB). During the measurement run, repeated analyses of NIST SRM987 standard solution yielded a ratio of 0.710250±0.000040 (2 SD; N=14), statistically consistent with the literature value of 0.710248±0.000058 (2 SE; McArthur et al., 2001; Weis et al., 2006). All results were corrected for instrumental mass discrimination by internal normalization and normalized to the literature value of NIST SRM987 (0.710248) through a standard-sample bracketing method. For each sample, ⁸⁷Sr∕⁸⁶Sr values are reported with 2σ uncertainty (Fig. S8 from de Winter et al., 2020c).

Figure 4 Plot showing the results of Sr-isotopic analyses with error bars (2σ) plotted on the Sr-isotope curve of McArthur et al. (2016; top of image). Numbers below the error bars indicate specimen number. Measurements from the 12 specimens of R. diluvianum are represented by parallelograms in different colors which match the probability distributions plotted below. The probability distribution curves in the lower panel show the distribution of uncertainty of each Sr-isotope measurement as well as the uncertainty of the Sr-isotope curve propagated to the age domain (colors of individual shells are the same as in Fig. 2). Inset shows schematically how uncertainties of the isotope measurements as well as the isotope curve are propagated into the age domain. The black curve shows the total uncertainty distribution function compiled from the 12 individual measurements following Barlow (2004), with the combined age estimate including uncertainty (2σ) shown above.

3.8 Strontium isotope dating

R. diluvianum specimens were independently dated by comparing ⁸⁷Sr∕⁸⁶Sr values measured in the samples with the Sr-isotope curve in the 2016 Geological Timescale (McArthur et al., 2016). Uncertainties in ⁸⁷Sr∕⁸⁶Sr measurements were propagated into dates by finding the closest date of the mean ⁸⁷Sr∕⁸⁶Sr value as well as the dates of the minimum (−2σ) and maximum (+2σ) ⁸⁷Sr∕⁸⁶Sr values by linearly interpolating ages in the ⁸⁷Sr∕⁸⁶Sr curve matching the measured ⁸⁷Sr∕⁸⁶Sr value, including the uncertainty estimated on the Sr-isotope curve itself. A composite age for the Ivö Klack deposits was obtained by combining the age uncertainty distributions of the individually dated ⁸⁷Sr∕⁸⁶Sr samples into a single age. Due to the non-linear shape of the ⁸⁷Sr∕⁸⁶Sr curve, uncertainties of the ⁸⁷Sr∕⁸⁶Sr ages were asymmetrical. Since no mathematical solution exists for the combination of asymmetric uncertainties, the asymmetric uncertainty in the total age had to be approximated through maximum likelihood estimation using the combined log likelihood function (Barlow, 2003). The approximation of the total uncertainty of combined ⁸⁷Sr∕⁸⁶Sr dating results in this study was carried out using the mathematical approach of Barlow (2004) in R (R Core Team, 2018; Roger Barlow, personal communication, 2019, code available at https://zenodo.org/record/1494909, last access: 27 May 2020). The uncertainty interval of the composite age is represented by 2 times the standard error (∼95.5 % confidence level). A plot of the uncertainty distributions of the individual specimens and the total uncertainty distribution is shown in Fig. 4. Raw ⁸⁷Sr∕⁸⁶Sr data are provided in Fig. S8 from de Winter et al. (2020c).

4 Results

4.1 Strontium isotope dating

Results of strontium isotope analyses are given in Fig. S8 from de Winter et al. (2020c). The mean strontium isotope ratio of all R. diluvianum specimens is 0.707552 (±0.000112; 95 % confidence level). The compilation of ⁸⁷Sr∕⁸⁶Sr results from 12 specimens of R. diluvianum (Fig. 4) shows how age estimates from individual specimens have considerable uncertainties (standard deviations around 1 Myr; see Fig. S8 from de Winter et al., 2020c), yet the uncertainty of the composite age is significantly smaller. The composite age for the Ivö Klack deposits is 78.14 Ma (±0.26; 95 % confidence level). This result places the age of the Ivö Klack deposits close to the early–late Campanian boundary when applying a 2-fold division of the Campanian and in the middle Campanian when applying a 3-fold division scheme (Ogg et al., 2016). This age estimate is similar to the age obtained when plotting the B. mammilatus zone on the recent integration schemes of the Campanian (Wendler, 2013). Earlier estimates (Christensen, 1997; Surlyk and Sørensen, 2010; Sørensen et al., 2015) yielded ages about 2–4 Myr older (80–82 Ma), but those relied on presently outdated and partly incorrect age models.

4.2 Shell structure and preservation

A combination of high-resolution color scans, SEM images and µXRF mapping of shell cross sections reveals that R. diluvianum shells consist of thin layers of dark, foliated calcite, interwoven with lighter, more porous carbonate layers (Fig. 3a–e). The latter are characterized by higher concentrations of Mn, Fe and Si and lower Sr concentrations (Fig. 3f–h). Foliated calcite layers are more densely packed on the inside of the shell, especially in the region of the adductor muscle scar, and at the shell hinge (Fig. 3a). They are characterized by high Sr concentrations and low concentrations of Mn, Fe and Si (Fig. 3f–h; Fig. S2 from de Winter et al., 2020c). Further away from the shell hinge and the inside of the valve, porous carbonate layers become more dominant. In these regions, µXRF mapping also clearly shows that detrital material (high in Si and Fe) is often found between the shell layers (Fig. 3f). SEM images show that the shell structure of R. diluvianum strongly resembles that of modern oyster species, as described in previous studies (Fig. 3b–e; Carriker et al., 1979; Surge et al., 2001; Ullmann et al., 2010, 2013; Zimmt et al., 2018). The major part of the shell consists of calcite, of which the foliated structures were sampled for chemical analyses in this study. As in modern oyster species, aragonite may originally have been deposited on the resilium of the shell, but this region is not considered for analyses (see outer tip of shell hinge in Fig. 3a; Stenzel, 1963; Carriker et al., 1979; Sørensen et al., 2012). Close similarities with modern oysters allow the inference that shell growth in R. diluvianum occurred in a similar way to modern oyster species like Ostrea edulis, Crassostrea virginica and Crassostrea gigas. On this basis we could estimate the total shell height at specific times in ontogeny from microstructural growth markers (purple lines in Fig. 3a; following Zimmt et al., 2018), linking growth to changes in shell chemistry.

4.3 Trace element results

The similarity in growth between R. diluvianum and modern oyster species is used to assess whether trace element variability in R. diluvianum can be interpreted in terms of environmental changes in a similar way to modern oyster shells (e.g., Surge and Lohmann, 2008; Ullmann et al., 2013; Mouchi et al., 2013; Dellinger et al., 2018). The combination of µXRF and LA-ICP-MS analyses on R. diluvianum shells resulted in records of changes in Mg∕Ca, Sr∕Ca (µXRF), Mg∕Li and Sr∕Li (LA-ICP-MS) as well as individual concentrations of trace elements such as Mg, Mn, Fe and Sr (Fig. 5). All chemical analyses were carried out on the dense foliated calcite exposed in cross sections close to the inner edge of the shell valve (Fig. 3a). High-resolution color scans and detailed recording of sampling positions allowed these records to be plotted on a common axis (see Figs. S6, S10 from de Winter et al., 2020c). In Fig. 5, results of chemical analyses of R. diluvianum specimens (including diagenetically altered parts) are compared with data from three other mollusk taxa (the belemnite B. mammilatus, the oyster Acutostrea incurva and the radiolithid rudist Biradiolites suecicus) from Ivö Klack (Sørensen et al., 2015) and data from extant oysters (Rucker and Valentine, 1961; Surge et al., 2001; Ullmann et al., 2013). Multi-proxy analysis revealed periodic variations in stable isotope and trace element ratios (see Fig. 6). The amplitudes of these variations plotted in Fig. 5 show that Mg and Sr concentrations measured in all three fossil bivalve taxa are similar, while concentrations in the belemnite rostra are much higher. Finally, plots in Fig. 5b and d demonstrate that diagenetic alteration (evident from elevated Mn concentrations) reduces the Sr concentration in carbonate of all four taxa. Stable oxygen isotope ratios of the shells are affected to a lesser degree (see below). The majority of measurements in all four taxa show very few signs of diagenetic alteration, with most measurements characterized by low (<100 µg g⁻¹) Mn concentrations (Fig. 5).

Figure 5 Cross plots summarizing the results of trace element and stable isotope analyses of the oysters R. diluvianum (blue), A. incurva (purple), the rudist Biradiolites suecicus (green) and the belemnite B. mammilatus (red; after Sørensen et al., 2015) from the Kristianstad Basin. Results from modern C. gigas (brown; Ullmann et al., 2013) and C. virginica (yellow; Rucker and Valentine, 1961; Surge et al., 2001) oysters are plotted for comparison. Points indicate individual data points; drawn polygons illustrate the range of the data, and crosses indicate the extent of seasonality (if present). (a) Strontium concentrations plotted against magnesium concentrations. (b) Strontium concentrations plotted against manganese concentrations. Arrows indicate the interpreted direction of diagenetic alteration, and the black dashed line shows the Sr∕Mn diagenesis threshold proposed for belemnite rostra by Sørensen et al. (2015; 100 mol mol⁻¹). (c) δ¹⁸O plotted against δ¹³C. Grey area indicates the range of stable isotope ratios measured in Campanian chalk deposits from the nearby Danish Basin by Thibault et al. (2016). (d) δ¹⁸O plotted against manganese concentrations, with arrows indicating proposed direction of diagenetic alteration.

Figure 6 Overview of stable oxygen (a) and carbon (b) profiles from all five shells in which stable isotope profiles were measured and plotted against ontogenetic ages. Black lines indicate five-point running averages through the time series. Vertical dashed lines separate years of growth. All vertical axes of δ¹⁸O and δ¹³C have the same scale. Color coding is the same as in Fig. 7.

4.4 Stable isotope records

An overview of stable isotope results of R. diluvianum (IRMS, only for shells 1–5; see Fig. 2) compared with the different taxa in Kristianstad Basin and modern oyster data is given in Fig. 5. Stable isotope ratios of the rudist and oyster shells overlap, while belemnite rostra of the species B. mammilatus are characterized by much lower δ¹³C and higher δ¹⁸O values. This suggests that δ¹³C values in belemnite rostra from this species are impacted by vital effects, while heavier δ¹⁸O values of the belemnites suggest that belemnites lived most of their life away from the coastal environment (in deeper waters), as previously suggested by Sørensen et al. (2015). In contrast, stable isotope ratios recorded in the bivalve shells overlap and match the isotope ratios measured in Campanian chalk deposited in the neighboring Danish Basin (Thibault et al., 2016). Records of δ¹³C and δ¹⁸O in the growth direction through R. diluvianum shells exhibit periodic variations (Fig. 6). These variations are much more regular in δ¹⁸O records, which show extreme values below −3 ‰ and up to 0 ‰ VPDB (Fig. 6a). Some shells, such as R. diluvianum 3 (Fig. 7), exhibit longer-term trends on which these periodic variations are superimposed. These trends suggest the presence of multi-annual cyclicity with a period on the order of 10–20 years, but the length of R. diluvianum records (max. 10 years) is smaller than the estimated period of these changes and is therefore insufficient to statistically validate the presence of this cyclicity.

The extreme values in δ¹⁸O records translate into tentative temperatures in the range of extremes of 12 to 26 ^∘C when assuming a constant δ¹⁸O_sw value of −1.0 ‰ (e.g., Thibault et al., 2016) and applying the temperature relationship of Kim and O'Neil (1997). However, the assumption of constant δ¹⁸O_sw may add bias to the temperature reconstructions, as seawater composition may not have been constant or reflect the marine value year-round in the nearshore Ivö Klack setting. Carbon isotope ratios (δ¹³C) do not always follow the same trends as δ¹⁸O records (Fig. 6b). In many parts of R. diluvianum shells, there is a clear covariation between the two isotope ratios, suggesting that δ¹³C is affected by seasonal changes. However, in other parts this correlation is less clear, suggesting that other (non-seasonal) factors play a role in determining the δ¹³C of shell material.

4.5 Age models

Modeling the growth of R. diluvianum bivalves from seasonal variations in δ¹⁸O profiles yielded age models, growth rate estimates and reconstructions of water temperature variations during the lifetime of the bivalves. Due to the clear seasonal patterns in δ¹⁸O records (Figs. 6a, 7), modeled δ¹⁸O profiles closely approximated the measured δ¹⁸O profiles (total R²=0.86, N=412; see Figs. S5 and S9 from de Winter et al., 2020c), lending high confidence to shell age models (see example in Fig. 7). Modeling allowed all proxies measured in the shells of R. diluvianum to be plotted against shell age, clearly revealing the influence of seasonal variations in environmental parameters on shell chemistry (Fig. S10 from de Winter et al., 2020c). The age models reveal clear, statistically significant (p<0.05) ontogenetic trends in Mg∕Li, Sr∕Li and δ¹³C in nearly all specimens of R. diluvianum (see Table 1). In three out of five shells, δ¹³C increases with age (see Fig. 6b and Table 1), with ontogenetically older specimens (e.g., R. diluvianum no. 2) yielding overall higher δ¹³C values (Fig. 6b). The distribution of slopes of ontogenetic trends in Li∕Ca cannot be distinguished from random variation. Therefore, no widely shared ontogenetic trends were found for Li proxies in R. diluvianum shells.

Table 1 Overview of the slopes of ontogenetic trends in Li∕Ca and δ¹³C records. p values on the bottom of the table show that the distribution of Li∕Ca slopes is statistically indistinguishable from random.

Figure 7 Example of multi-proxy records measured in R. diluvianum specimen 3 plotted against distance in growth direction (see image on top and Fig. 3 for reference). From top to bottom, records of Sr∕Li (light blue), Mg∕Li (light green), Sr∕Ca (dark blue), Mg∕Ca (dark green), δ¹³C (red), δ¹⁸O (blue dots with error bars) and modeled growth rate (light green shading) are plotted. The shaded blue curve plotted underneath the δ¹⁸O record illustrates the result of growth and δ¹⁸O modeling and its propagated error (vertical thickness of curve, 2 SD). The dashed blue curve plotted on top of the δ¹⁸O record shows the observed multi-annual trend in the data. Vertical dashed lines separate growth years. Multi-proxy plots for all specimens are given in Fig. S10 from de Winter et al. (2020c).

5 Discussion

5.1 Preservation

The relative lack of burial and tectonic activity in the Kristianstad Basin has provided ideal circumstances for the nearly immaculate preservation of R. diluvianum shells in the Ivö Klack locality (Kominz et al., 2008; Surlyk and Sørensen, 2010). The excellent state of these shells is evident in the preservation of original (porous and foliated) microstructures that closely resemble those reported for several species of modern ostreid shells (Carriker et al., 1979; Surge et al., 2001; Ullmann et al., 2010, 2013; Zimmt et al., 2018; Figs. 2–3). High-magnification SEM images demonstrate the excellent preservation of foliated and vesicular calcite structures in R. diluvianum shells (Fig. 3b–d). The preservation state of R. diluvianum shells meets the SEM-based preservation criteria for robust stable isotope analysis set by Cochran et al. (2010). Micro-XRF mapping reveals that the foliated calcite in the shells is characterized by high Sr concentrations and low concentrations of Mn, Fe and Si, elements which are generally associated with diagenetic alteration (e.g., Brand and Veizer, 1980; Al-Aasm and Veizer, 1986a; Immenhauser et al., 2005; Fig. 3b–h). Trends in bulk Mn and Sr concentrations observed in all fossil species from Ivö Klack (Fig. 5b; including less well preserved parts) likely point towards a diagenetic process affecting a less well preserved subset of the data. The lack of covariation between Mn concentration and δ¹⁸O shows that there is little evidence for meteoric diagenesis in these shells (Fig. 5d; Ullmann and Korte, 2015). Instead, these patterns are best explained by early marine cementation of porous carbonate structures from seawater with similar temperature and δ¹⁸O to the living environment (see also Sørensen et al., 2015).

Typically, a Mn concentration threshold of 100 µg g⁻¹ is applied, below which Cretaceous low-magnesium carbonates are assumed to be suitable for chemical analysis (Steuber et al., 2005; Huck et al., 2011). Strontium concentrations above 1000 µg g⁻¹ have also been used as markers for good preservation, since diagenetic processes can cause strontium to leach out of carbonates (e.g., Brand and Veizer, 1980; Huck et al., 2011; Ullmann and Korte, 2015). Previous studies of belemnites in Kristianstad Basin proposed a molar Sr∕Mn threshold of 100 for belemnite rostra (Sørensen et al., 2015). However, the height of thresholds used for diagenetic screening differs widely in the literature due to variation between species, fossil matrices and burial histories (e.g., Veizer, 1983; Steuber et al., 1999; Ullmann and Korte, 2015; de Winter and Claeys, 2016). Applying these thresholds risks introducing biases to chemical datasets from fossil shells and may not be an ideal method for diagenetic screening. Furthermore, large variation in the in vivo incorporation of Mn and Sr in mollusk shell carbonate and a strong dependence on the diagenetic setting can make the interpretation of shell preservation from trace element ratios alone highly ambiguous (Ullmann and Korte, 2015). The complex patterns in multi-proxy datasets in this study (Fig. 5) merit great care in applying simple, general thresholds for grouping different processes of carbonate diagenesis. Therefore, in this study, a multi-proxy approach is applied for diagenetic screening in which data are excluded based on a combination of Si, Ca, Mn, Fe and Sr concentrations; δ¹⁸O values; and SEM and visual observations of the shell structure at the location of measurement.

5.2 Dating the Ivö Klack locality

While previous attempts at dating Campanian strata mainly focused on relative dating using magneto- and biostratigraphy (Montgomery et al., 1998; Jarvis et al., 2002; Voigt et al., 2010), integration of cyclostratigraphic approaches has recently allowed constraining the age of the Campanian deposits more precisely (Voigt and Schönfield, 2010; Thibault et al., 2012; Wendler, 2013; Thibault et al. 2016). Unfortunately, these rarely cover the time interval in which the Ivö Klack sediments were deposited (e.g., Wendler, 2013; Perdiou et al., 2016). Strontium isotope dating places the Ivö Klack deposits at 78.14±0.26 Ma (Fig. 4), slightly above the early–late Campanian boundary, while the B. mammilatus biozone is defined as late early Campanian (Wendler, 2013). An influx of radiogenic strontium-rich weathering products from the nearby Transscandinavian Igneous Belt may bias age estimates from strontium isotope ratios (Högdahl et al., 2004). However, studies of modern strontium isotope ratio variability (Palmer and Edmond, 1989) and the potential bias of strontium isotope ratios in shallow-water carbonates (Kuznetsov et al., 2012; El Meknassi et al., 2018) show that the effect of such inputs on strontium isotope dating results is generally negligible, except in semi-confined shallow-marine basins characterized by considerable freshwater input and low salinities (<7 psu; psu – practical salinity units). No evidence exists for such exceptional conditions at Ivö Klack (see Sect. 2.1). We therefore conclude that our strontium isotope age estimate, together with biostratigraphic constraints, places the Ivö Klack locality in the earliest late Campanian. The refined dating of the Ivö Klack deposits and fossils allows the results of sclerochronological investigations presented in this work to be placed in the context of longer-term climate reconstructions with improved precision.

5.3 Trace element variability

Extracted ranges in seasonal-scale periodic variability in Mg∕Ca, Sr∕Ca, Mg∕Li and Sr∕Li in all 12 R. diluvianum shells (Fig. 8) show that interpreting these records in terms of temperature changes is not straightforward. Some of this difficulty arises from the large inter-shell variability in trace element ranges, mostly expressed in strong positive correlations between Sr∕Li and Mg∕Ca (R²=0.76) and between Sr∕Li and Mg∕Li (R²=0.93). The benthic-foraminifera proxy transfer function for Mg∕Li (Bryan and Marchitto, 2008) does not work for R. diluvianum (temperatures >50 ^∘C), presumably due to typically lower Mg concentrations in foraminifera compared to bivalves (Yoshimura et al., 2011). Bivalve-specific temperature relationships of Mg∕Ca (Surge and Lohmann, 2008; based on Crassostrea virginica), Sr∕Li (Füllenbach et al., 2015; based on Cerastoderma edule) and Mg∕Li (Dellinger et al., 2018; based on Mytilus edulis) yield temperatures in the same range as those reconstructed from stable oxygen isotope measurements (10–20 ^∘C). However, Sr∕Li-based temperature trends are opposite to those based on Mg proxies, suggesting that they cannot both be applicable to R. diluvianum. Poorly constrained changes in seawater chemistry (Mg∕Ca and Sr∕Ca ratios of ocean water) also hinder these trace element-based reconstructions (Lear et al., 2003; Coggon et al., 2010; Rausch et al., 2013). The strong Mg∕Li–Sr∕Li correlation indicates that both proxies are likely strongly affected by the specimen-specific ontogenetic trends in Li∕Ca described in Table 1. This, together with the large inter-specimen variability, shows that both Li proxies cannot be used as temperature proxies in R. diluvianum. An annual stack of all proxies shows that the positive correlation between Mg∕Ca and δ¹⁸O (Fig. 9) is opposite to the temperature relationships found in modern oyster species (Surge and Lohmann, 2008; Mouchi et al., 2013; Ullmann et al., 2013). This, together with the reduced seasonal variability (1.2 mmol mol⁻¹ versus 4–10 mmol mol⁻¹ in modern oysters; Surge and Lohmann, 2008; Mouchi et al., 2013) and the large (>3 mmol mol⁻¹; Fig. 8) inter-specimen variability both in mean value and seasonal Mg∕Ca range, rules out Mg∕Ca as a reliable temperature proxy in R. diluvianum. This result demonstrates that earlier successful attempts to establish calibration curves for Li- and Mg-based temperature proxies (e.g., Füllenbach et al., 2015; Dellinger et al., 2018) are probably strictly limited to these bivalve species or close relatives. The same conclusion was also drawn by Dellinger et al. (2018) based on Li∕Mg and Li isotope ratio measurements in biogenic carbonates. The lack of Mg∕Li or Sr∕Li calibrations in modern oyster shells limits the interpretation of our results for these element ratios. Establishing such calibrations using modern oysters in cultured experiments may allow these proxies to be used for reconstructions from fossil oyster shells in the future.

While not a promising candidate for reconstructing temperature (Gillikin et al., 2005; Schöne et al., 2013; Ullmann et al., 2013), seasonal Sr∕Ca fluctuations and relationships between Sr∕Ca and δ¹⁸O are consistent between individuals (Figs. 8–9; see also Fig. S6 from de Winter et al., 2020c). This allows Sr∕Ca ratios to be used together with microstructural observations of growth increments as a basis for seasonal-scale age models in shells for which no δ¹⁸O measurements were performed. Both mean Sr∕Ca values and seasonal variability in R. diluvianum are consistent with those observed in the same microstructure in modern Crassostrea gigas growing in a similar environment (0.8–1.0 mmol mol⁻¹; Ullmann et al., 2013), suggesting a consistent incorporation of Sr by different oyster taxa over time. It must be noted that one should be cautious when directly comparing trace element concentrations in biogenic calcite between different time periods, as the seawater composition of Late Cretaceous oceans (e.g., concentrations of Mg, Ca, Sr and especially Li) may have been different from that of the present-day ocean (Stanley and Hardie, 1998; Coggon et al., 2010; Rausch et al., 2013). Local enrichments in seawater Sr concentrations at Ivö Klack driving increased Sr content in R. diluvianum are unlikely, since Sr∕Ca ratios exhibit only small (2 %–3 %) lateral variability in the world's oceans (De Villiers, 1999). Because Sr∕Ca ratios in Late Cretaceous oceans were half those of the modern ocean, one would expect that Sr concentrations in Late Cretaceous biogenic carbonate would be half those in carbonates formed in the modern ocean if the partition coefficient between seawater concentrations and shell concentrations remained constant (Stanley and Hardie, 1998; de Winter et al., 2017a). The fact that this reduction in Sr concentrations relative to the modern ocean is not observed in R. diluvianum may mean that there is a physiologically fixed concentration of Sr in an oyster's shell that is independent of ambient Sr concentrations.

Figure 8 Cross plots showing the extent of interpreted seasonality observed in records of four trace element proxies in all 12 R. diluvianum specimens. Colors of lines of individual shells correspond to colors indicated in Fig. 2. Temperature conversions from previously published regressions of the proxies with temperature are shown on opposite axes, with grey dashed lines corresponding to major tick marks on the temperature scale. (a) Sr∕Li plotted against Mg∕Ca, showing a strong significant intra-shell correlation. (b) Sr∕Li plotted against Mg∕Li, showing a strong significant intra-shell correlation due to dominant variations in Li concentration. Note that two different Mg∕Li temperature calibrations were explored. (c) Sr∕Ca plotted against Mg∕Ca, showing weak but significant intra-shell correlation. (d) Sr∕Ca plotted against Mg∕Li, showing a weakly significant intra-shell correlation. Data for this plot are found in Fig. S13 from de Winter et al. (2020c).

Figure 9 Composite of multi-proxy records from all R. diluvianum shells stacked and plotted on a common time axis of 1 year to illustrate the general phase relationships between various proxies in the shells. Records were colored as in Fig. 7. Annual stacks plotted in this figure were produced by applying age models on all multi-proxy records, plotting all results against their position relative to the annual cycle and applying 20-point moving averages.

5.4 Temperature seasonality

An annual stack of all R. diluvianum proxy records shows a δ¹⁸O_c-based temperature variability in Ivö Klack of 16–21 ^∘C when assuming a constant δ¹⁸O_sw of −1 ‰ VSMOW (Fig. 6). However, comparison with δ¹⁸O seasonality in individual specimens shows that the annual stack severely dampens seasonality due to small phase shifts in maximum and minimum temperatures, small uncertainties in the age models between years and specimens, and inter-annual differences and longer-term trends in temperature (see Fig. 6). A more accurate estimate of the seasonal extent is obtained by calculating the seasonal range from the coolest winter temperatures (12.6 ^∘C in R. diluvianum 4; Fig. 6 and Fig. S10 from de Winter et al., 2020c) with the warmest recorded summer temperature (26 ^∘C in R. diluvianum 1; Fig. 6 and Fig. S10 from de Winter et al., 2020c), which yields an extreme maximum seasonal sea surface temperature (SST) range of ±13.4 ^∘C.

A complication of these reconstructions is the assumption of constant δ¹⁸O_sw of −1 ‰ VSMOW, which is unlikely to be completely accurate in the nearshore Ivö Klack locality. Comparison with data from Crassostrea gigas growing in a similar nearshore environment (Ullmann et al., 2010; German North Sea coast, 54^∘ N) shows that such an environment away from large river mouths can typically experience seasonal salinity fluctuations of ∼4 psu, resulting in a dampening of the seasonal δ¹⁸O_c cycle by ∼0.5 ‰ VPDB. Such a salinity effect would reduce our reconstructed 13–26 ^∘C seasonal temperature range by ∼2 ^∘C to 14–25 ^∘C.

In addition, mean annual δ¹⁸O_sw values can be considerably lower than the global average seawater composition (e.g., −1 ‰ to −1.5 ‰ VSMOW compared to global ocean mean of 0 ‰ VSMOW in Ullmann et al., 2010). Considering such a deviation would reduce reconstructed temperatures by 4–6 to 10–21 ^∘C, much colder than open marine reconstructions of the Boreal Chalk Sea by Thibault et al. (2016). This result would be in strong disagreement with a recent study by Tagliavento et al. (2019) in which clumped isotope analyses (which do not rely on the assumption of constant δ¹⁸O_sw) were used to correct the δ¹⁸O_c-based reconstructions of the boreal chalk and yielded higher temperatures (∼26 ^∘C MAT for open marine SST) and a correction of δ¹⁸O_sw towards 1 ‰–1.5 ‰ heavier values (resulting in a Campanian δ¹⁸O_sw of −0.5 ‰–0 ‰ VSMOW). Another caveat is that salinity effects on local δ¹⁸O_sw strongly depend on the local δ¹⁸O_sw of the local freshwater source (riverine or precipitation), which in the present-day higher mid-latitudes is around −7 ‰ VSMOW to −8 ‰ VSMOW (e.g., Ullmann et al., 2010), but this is impossible to constrain at Ivö Klack during the Campanian within the scope of this study.

If local δ¹⁸O_sw values at Ivö Klack were indeed 1 ‰–1.5 ‰ reduced with respect to those in the fully marine Boreal Chalk Sea, and marine δ¹⁸O_sw was around 0-0.5 ‰ VSMOW rather than the assumed −1 ‰ VSMOW, the effects of these two biases cancel each other out, and the best estimation of the extreme seasonal SST range at Ivö Klack based on this study's data would be 14–25 ^∘C with a MAT of 19 ^∘C. This MAT is comparable to the MAT of the late early Campanian Boreal Chalk Sea waters of 17–19 ^∘C calculated from coccolith δ¹⁸O_c (Lowenstam and Epstein, 1954; Jenkyns et al., 2004; Friedrich et al., 2005; Thibault et al., 2016) and slightly warmer than mean annual air temperatures from this paleolatitude based on phosphate-δ¹⁸O reconstructions (±15 ^∘C; Amiot et al., 2004). However, Ivö Klack SSTs are ∼6 ^∘C colder than the clumped isotope-based reconstructions from marine chalk samples (Tagliavento et al., 2019). The latter could indicate that coastal SSTs and air temperatures were much colder than marine temperatures in the Campanian higher latitudes, but such temperature differences are highly unusual compared to modern climates. Alternatively, this difference could highlight a severe temperature bias in (both phosphate and carbonate) δ¹⁸O-based reconstructions, which should be further investigated using independent marine temperature proxies such as the clumped isotope paleothermometer (e.g., de Winter et al., 2018; Tagliavento et al., 2019).

Modeled growth rates in R. diluvianum peak near the end of the low-temperature season, and average growth rates are lowest shortly after the temperature maximum (Fig. 9). This phase shift between temperature and growth rate could indicate that growth in R. diluvianum in this setting was not limited by low temperatures, as observed in modern mid- to high-latitude oysters (Lartaud et al., 2010). High temperature extremes (>25 ^∘C) may have slowed or stopped growth, as recorded in modern low-latitude settings (Surge et al., 2001). Heat shock has been shown to limit the growth of modern oysters (Crassostrea gigas; Li et al., 2007), although the relatively moderate SST seasonality suggests that very high (>25 ^∘C) temperatures were not common at the Ivö Klack locality (Fig. 6). In this respect it is important to recognize that the use of δ¹⁸O records from multiple specimens reduces the effects of growth cessations of individuals on seasonal SST reconstructions and allows the full seasonal range in SST to be resolved.

The reconstructed MAT for Ivö Klack is 7–8^∘ warmer than the present-day seasonal SST range in the North and Baltic seas at similar latitude (e.g., 2–18 ^∘C monthly seasonality range in Baltic Sea, Karlskrona, Sweden, 56^∘ N; and 4–18 ^∘C monthly seasonality range in North Sea, Esbjerg, Denmark, 55^∘ N; International Research Institute for Climate and Society, 2020). Seasonal SST range at Ivö Klack (11 ^∘C) is significantly lower than the 14–16 ^∘C maximum seasonal range in temperature that occurs in the present-day Baltic and North seas (International Research Institute for Climate and Society, 2020). Data on temperature seasonality in the Late Cretaceous are scarce, especially in high-latitude settings. However, comparison with data on Cretaceous temperature seasonality between 15 and 35^∘ N paleolatitude (Steuber et al., 2005) shows that while MAT at 50^∘ N was significantly lower than at lower latitudes (19 and 25–30 ^∘C, respectively), the seasonal temperature range during cooler periods in the Late Cretaceous was remarkably similar between latitudes (10–15 ^∘C in subtropical latitudes vs. ±14 ^∘C in this study). This observation contrasts with the present-day situation in northern Africa and Europe, in which seasonal temperature ranges are generally much higher at mid-latitudes to high latitudes (30–50^∘ N) than at lower latitudes (10–30^∘ N; Prandle and Lane, 1995; Rayner, 2003; Locarnini et al., 2013; NOAA, 2020). Our SST reconstructions also show that Late Cretaceous latitudinal temperature gradients and mid- to high-latitude seasonality were larger than previously assumed based on climate model results (Barrera and Johnson, 1999; Hay and Floegel, 2012; Upchurch et al., 2015).

5.5 Shell growth and ontogeny

5.5.1 Growth curves

Growth curves of individual R. diluvianum specimens clearly converge to a general growth development curve for the species (Fig. 10). Considering that these growth curves are based on δ¹⁸O and Sr∕Ca transects in different stages of life in different specimens (large age variation; see Fig. 8), individual growth curves are remarkably similar. The growth of R. diluvianum takes the typical shape of the asymptotic von Bertalanffy curve, in which shell height (H_t) development with age (t) is related to a theoretical adult size H_max and a constant k in the equation $H_t\left[\mathrm{mm}\right]=H_{\max }\cdot \left(\mathrm{1}-e^{k\cdot (t\left[\mathrm{yr}\right]-t_{\mathrm{0}})}\right)$, with t₀ representing the time at which the growth period started (always zero in this case; von Bertalanffy, 1957). When this formula is regressed over all modeled growth data of all shells (1 data point per day, 15 146 points in total), the fit with an H_max of ±120.3 mm and a k value of ±0.32 is very good (R²=0.79; see Fig. 10).

The consistency in growth curves between individuals of R. diluvianum is surprising, as modern oyster species are known to exhibit large variations in growth rates and shell shapes as a function of their colonial lifestyle, which often limits the growth of their shells in several directions (Galtsoff, 1964; Palmer and Carriker, 1979). The strong resemblance of growth between individuals and the close fit of the idealized von Bertalanffy growth model suggests that growth of R. diluvianum at Ivö Klack was relatively unrestricted in space. This hypothesis is consistent with the apparent mode of life of R. diluvianum in Ivö Klack, cemented next to each other in loose groups, subject to strong wave action and turbulence, such that individuals received the same (high) supply of food (Surlyk and Christensen, 1974; Sørensen et al., 2012). The shape of the growth curve of R. diluvianum resembles that of modern Chesapeake Bay oysters (Crassostrea virginica), which exhibit a slightly larger modeled maximum height (150 mm) and a slightly smaller K value (0.28). A larger subset of R. diluvianum specimens studied by Sørensen et al. (2012) demonstrates that these bivalves could grow up to 160 mm in height. The curvature of the growth curve of R. diluvianum (K value) is also similar to that found for other modern shallow-marine bivalve species (e.g., Macoma balthica, K=0.2–0.4; Bachelet, 1980; Pinna nobilis, K=0.33–0.37; Richardson et al., 2004) and significantly higher than in growth curves of deep shelf-dwelling bivalves (e.g., Placopecten magellanicus, K=0.16–0.24; MacDonald and Thompson, 1985; Hart and Chute, 2009) or bivalves from cold habitats (e.g., North Atlantic Arctica islandica, K=0.06; Strahl et al., 2007). This reflects the high growth rates (steeper growth curves, higher K values) of shallow-marine bivalves compared to species living in less favorable (colder or deeper) habitats, with R. diluvianum clearly being part of the former group.

Figure 10 Shell height plotted against age for all R. diluvianum records (see Fig. 8 for color legend of lines representing individuals). The similarity between growth curves of different specimens allows a von Bertalanffy curve to be fitted to the data with high confidence. Sinusoidal patterns superimposed on all growth curves are caused by seasonal variability in growth rate (see Figs. 6–7). Data found in Fig. S9 from de Winter et al. (2020c).

5.5.2 Seasonal growth

To study variability in minimum growth temperature (Fig. 11a), length of the growth season and time of year at which maximum growth occurs (Fig. 11b), we isolated individual growth years from all age models of the five shells in which δ¹⁸O curves were measured (Fig. 11). The onset and end of each growth year correspond to maxima in δ¹⁸O values (minima in temperatures). The onset of the first growth year in each shell at its precise position relative to the seasonal temperature cycle showed the season in which the individual settled and started growing its shell (Fig. 11b). All data used to create plots in Fig. 11 are provided in Fig. S14 from de Winter et al. (2020c). Relationships between these growth parameters are summarized in Table 2.

Figure 11 Overview of statistical evaluation of growth parameters of R. diluvianum derived from age modeling in shells 1–5. (a) Histogram of minimum temperatures of growth (equivalent to the time at which growth stops or the minimum yearly temperature) in R. diluvianum, showing that the temperature at which growth slows coincides with the time of first growth at post-larval settlement (p=0.717). (b) Histogram of the season of maximum growth relative to the δ¹⁸O seasonality cycle shows no significant concentration towards a favorable growing season, while moments of first growth (following settlement after the larval stage) are significantly concentrated towards the high-δ¹⁸O season.

Table 2 Overview of statistical evaluation of growth parameters of R. diluvianum derived from age modeling in shells 1–5. Coefficients of determination (R²) and p values were determined for relationships between temperature seasonality, average temperature, age of the bivalve, length of the season, minimum growth temperatures, and annual average and maximum growth rates. Values in italics indicate strong correlations, while values in bold indicate the absence of a correlation. Data reported in Fig. S14 from de Winter et al. (2020c).

The growing season is shorter than 365 d in all but 5 modeled years, demonstrating that growth stops or slowdowns did occur in R. diluvianum. Minimum growth temperatures (temperatures at which growth stopped) are concentrated around 17 ^∘C (χ²=0.0088; Fig. 11a) and correlate strongly to MAT (R²=0.57; Table 2), suggesting that growth halts in R. diluvianum were not forced by an absolute temperature threshold but rather by timing relative to the seasonality (circadian rhythm). If there had been a fixed temperature threshold (e.g., 6 ^∘C or 10 ^∘C for Crassostrea gigas; Lartaud et al., 2010; Ullmann et al., 2010) the length of the growth season would have positively correlated with annual mean temperature, but this is not the case (R²<0.1). Other authors have suggested that growth in modern C. gigas does not actually cease completely but rather slows down significantly, which may also have been the case in R. diluvianum. On average, the moment of minimum growth in R. diluvianum occurs right after the highest temperatures of the year are reached (early autumn; Fig. 9).

The season of first growth (after larval stage) is concentrated in the 2 last months before the δ¹⁸O maximum (first half of winter), when modeled water temperatures are ±17 ^∘C (Fig. 11b). Note that only three of the five shells allowed sampling of the first month of growth, and extrapolated records for the other two shells are likely unreliable. Comparing Fig. 11a and Fig. 11b shows that growth halts and settling occur at similar temperatures (16.9±0.7 and 17.0±0.3 ^∘C, respectively; p=0.72), suggesting that these events occur simultaneously or on either side of a seasonal growth halt (if it occurs). The timing of spawning and settling in R. diluvianum differs from that in modern oysters, which generally spawn during the spring season, with the young oyster spat settling in the summer (e.g., for Crassostrea gigas: Fan et al., 2011). If the larval stage of R. diluvianum lasted as long, this would put the spawning season of R. diluvianum in summer. In the case of modern oysters, it is known that reproduction requires a relatively warm minimum temperature (ideally around 22 ^∘C for C. gigas; Cognie et al., 2006) and that a combination of salinity and temperature conditions is important (Fan et al., 2011), while extreme temperatures (>28 ^∘C; Surge et al., 2001) can induce shock. Similar temperature conditions occurred at Ivö Klack in summer, possibly explaining the delay of the spawning season of R. diluvianum.

Figure 11b shows that the distribution of months with fastest growth rate is random (p(χ²)=0.055, <95 % confidence). However, in 27 of the 58 years, the growth peak occurs in the season with decreasing δ¹⁸O values (“spring season”). Table 2 shows that the extent of temperature seasonality (difference between minimum and maximum δ¹⁸O converted to temperature) significantly influences the length of the growing season (R²=0.40). However, total annual growth and maximum growth rates are independent of SST (both seasonal range and MAT), and ontogenetic age of the organism does not predict a significant part of any of the above-mentioned growth and seasonality parameters (Table 2). If temperature controlled the growth of R. diluvianum shells, larger temperature seasonality would increase the chance of crossing thresholds of temperature tolerance, which would shorten the length of the growing season. All this suggests that temperature seasonality may not have been the dominant factor causing growth cessations in R. diluvianum. This hypothesis is supported by the observation that temperatures at which growth cessations occur (16.9±0.7 ^∘C; Fig. 11a) show large variability and do not correspond significantly to the lowest temperatures of the year. A possible explanation might be that these growth cessations occurred in response to attempted predation or extreme weather events, which are not necessarily paced to the seasonal cycle.

These observations do not necessarily show that R. diluvianum tolerated larger temperature differences than modern oyster taxa because the extent of seasonality (14–25 ^∘C) causes neither the lower (∼10 ^∘C) nor the upper (∼28 ^∘C) temperature limit between which shell growth occurs in modern oysters to be reached. If temperature tolerance of R. diluvianum did resemble that of its closest modern relatives, then the mild seasonal temperatures at Ivö Klack might have provided the ideal temperature range for its growth. Perhaps these favorable conditions partly explain why biodiversity and abundance of invertebrates at Ivö Klack were so high (Surlyk and Sørensen, 2010).

5.5.3 Productivity

Shell growth in R. diluvianum may not have been governed by temperature but rather by changes in productivity. The observation that peak growth rates and settling both occur during the early spring or late autumn season (before or after the growth cessation; Fig. 11b) supports this hypothesis. Spring or autumn productivity blooms caused by increases in nutrient-rich fresh water from onshore (Arthur et al., 1983; Krantz et al., 1987) or storm-induced mixing of more nutrient-rich deeper waters are common in present-day mid- and high-latitude marine ecosystems (e.g., Waniek, 2003; Danielsson et al., 2008). An increase in seasonal freshwater influx would cause longer growth cessations to occur in the spring season, reducing the length of the growing season while also dampening the temperature seasonality reconstructed from δ¹⁸O due to the influx of isotopically light fresh water which dampens the seasonal cycle (see Sect. 5.5.2). This explains the correlation found between these two parameters (Table 2). At the same time, this freshwater input would increase reconstructed MAT by increasing δ¹⁸O values in R. diluvianum shells, explaining the weak positive correlation between MAT and length of the growing season (Table 2).

The occurrence of spring blooms is supported by weak 0.5 ‰–1.0 ‰ seasonal variability in δ¹³C (Fig. 6). Seasonal changes in productivity and/or salinity will cause changes in DIC in the environment, which are apparent in the δ¹³C of the shell above the ontogenetic trends and the effect of respiration on δ¹³C (see Sect. 2.4; Table 1). The fact that a clear seasonality in δ¹³C is absent from the stack in Fig. 9 shows that these productivity peaks do not occur at regular times of the year and that their effect on δ¹³C is obscured by ontogenetic trends. The 0.5 ‰–1.0 ‰ shifts in δ¹³C that appear to be seasonal are much smaller than those in modern oyster records (2 ‰–3 ‰ in low-latitude estuarine Crassostrea virginica; Surge et al., 2001, 2003; Surge and Lohmann, 2008). Instead, the determined shifts more closely resemble the 0.5 ‰ variability in δ¹³C observed in modern Crassostrea gigas from the same approximate latitude as Ivö Klack in the North Sea (Ullmann et al., 2013). The extreme isotopic shifts in the estuarine C. virginica specimens have been shown to be caused by large shifts in freshwater input due to large seasonal variations in rainfall over southern North America (Surge et al., 2003), while smaller variations in C. gigas from the North Sea are produced by DIC changes due to seasonal changes in productivity (e.g., spring blooms; Ullmann et al., 2013), indirect effects of temperature variability (Chauvaud et al., 2011) and changes in food source (Marchais et al., 2015). The closer resemblance of R. diluvianum to the North Sea condition shows that the Ivö Klack paleoenvironment did not experience large seasonal shifts in freshwater input and may have seen productivity peaks in spring season. The latter interpretation is in agreement with the coincidence of negative δ¹³C excursions (in parts of the records not affected by ontogenetic trends and respiration), with the high-δ¹⁸O season (winter or spring; Fig. 6 and Fig. S6 from de Winter et al., 2020c), the settling of larvae and the onset of shell growth, and a peak in growth rates in the spring season (much like in wild modern oysters; Berthelin et al., 2000; Figs. 9, 11a).

5.5.4 Ontogeny

A part of the variation in δ¹³C is explained by the presence of ontogenetic trends. These trends are known to occur in marine and freshwater bivalves, including bivalves with symbionts (Klein et al., 1996b; Watanabe et al., 2004; Gillikin et al., 2007; McConnaughey and Gillikin, 2008). The scale and direction of the trends in δ¹³C are not consistent between individual R. diluvianum shells, which is also the case in other bivalve species (see Sect. 4.5; Table 1; McConnaughey and Gillikin, 2008, and references therein). Studies of modern bivalves show that in larger (older) bivalves, the contribution of respired CO₂ to carbon in the shell is larger (up to 40 %; Gillikin et al., 2007). This finding explains common trends of reducing δ¹³C with age in bivalve shells, since respired carbon is isotopically lighter than environmental DIC. Since ontogenetic trends are likely caused by changes in the amount of respired carbon entering the shell, but the direction of these trends is inconsistent in R. diluvianum, the contribution of respired CO₂ to R. diluvianum shells likely did not strictly increase with age. While this complicates the interpretation of δ¹³C records, the relative contribution of environmental changes to δ¹³C variability in R. diluvianum shells does appear to be highest at the positive end of the ontogenetic trend.

6 Conclusions

The highly biodiverse marine invertebrate community at Ivö Klack in the Kristianstad Basin in southern Sweden offers a unique opportunity to recover a wealth of information about Campanian climate and environment in high latitudes and the ecology and life of extinct invertebrate species that lived under these conditions. The lack of burial and tectonic activity in the region favored Rastellum diluvianum fossil shells from Ivö Klack that were well preserved, as is evident from the excellent preservation of growth structures typical for ostreid shells as well as from the limited evidence for geochemical changes associated with diagenetic alteration. This excellent preservation allows the shells of R. diluvianum to be used to accurately and precisely constrain the age of the Ivö Klack locality using strontium isotope stratigraphy (78.14±0.26 Ma; ⁸⁷Sr∕⁸⁶Sr $=\mathrm{0.707552}\pm \mathrm{0.000112}$). Furthermore, R. diluvianum shells reveal sub-annual-scale variability in temperature, local environment and growth rates through a multi-proxy geochemical approach. The combination of trace element and stable isotope measurements with growth modeling based on δ¹⁸O records in the shells allows all measured proxies to be aligned on the same time axis. Application of transfer functions for potential Mg∕Ca, Mg∕Li and Sr∕Li temperature proxies established in modern invertebrates yields temperatures consistent with those calculated from δ¹⁸O records. However, close examination of the seasonal-phase relationships between these proxies reveals that the sub-annual variability in these trace element ratios is not controlled by temperature changes alone. This observation supports previous studies that found the expression of trace element proxies to be highly variable among species and even among different specimens of the same species. If trace element proxies are to be used for seasonality reconstructions in pre-Quaternary times, a more robust, non-species-specific model for the incorporation of trace elements by bivalves is required. Establishing such a model requires culture experiments with different bivalve species in which multiple parameters influencing trace element composition can be controlled (e.g., temperature, salinity, food intake and microstructure).

Stable isotope records in R. diluvianum shells reveal a MAT of 19 ^∘C with a seasonal water temperature range of ±11 ^∘C (14–25 ^∘C) at Ivö Klack. This value for MAT is lower than full marine SST of Boreal Chalk, recently re-evaluated with the clumped isotope thermometer. The difference highlights potential biases in temperature reconstructions based on δ¹⁸O values and argues for re-evaluations of these proxy records with more accurate techniques such as clumped isotope analysis. Comparing the seasonal temperature range reconstructed from R. diluvianum shells with other Late Cretaceous seasonality records from lower latitudes reveals that both latitudinal gradients and SST seasonality outside the tropics were much higher than predicted by climate models. This disagreement between data and models clearly illustrates the need for more data on Late Cretaceous seasonality outside the (sub-)tropical latitudes and highlights how important such proxy-based reconstructions are for improving our understanding of the dynamics in temperature variability in both space and time during greenhouse climates.

Finally, the coupled modeling and multi-proxy approach applied in this study sheds light on the effects of environmental changes on the life cycle and sub-annual growth of R. diluvianum shells. This study reveals that growth curves of R. diluvianum strongly resemble those in modern shallow-marine bivalves that grow in coastal high-latitude environments. However, ontogenetic changes in growth rate of our boreal oysters seem unrelated to temperature, in contrast to modern, high-latitude oysters that tend to lower their growth rate and cease mineralization below a certain cold threshold. We conclude that growth cessations and sub-annual changes in growth rate in R. diluvianum were most likely not caused by intolerable temperatures but rather by the circadian rhythm tied to the seasonal cycle and seasonal changes in sea surface productivity, driven by nutrient-rich freshwater inputs.