The fossil bivalve Angulus benedeni benedeni: a potential seasonally resolved stable-isotope-based climate archive to investigate Pliocene temperatures in the southern North Sea basin

Bivalves record seasonal environmental changes in their shells, making them excellent climate archives. However, not every bivalve can be used for this end. The shells have to grow fast enough so that micrometre- to millimetre-sampling can resolve sub-annual changes. Here, we investigate whether the bivalve Angulus benedeni benedeni is suitable as a climate archive. For this, we use ca. 3-million-year-old specimens from the Piacenzian collected from a temporary outcrop in the Port of Antwerp area (Belgium). The subspecies is common in Pliocene North Sea basin deposits, but its lineage dates back to the late Oligocene and has therefore great potential as a high-resolution archive. A detailed assessment of the preservation of the shell material by micro-X-ray fluorescence, X-ray diffraction, and electron backscatter diffraction reveals that it is pristine and not affected by diagenetic processes. Oxygen isotope analysis and microscopy indicate that the species had a longevity of up to a decade or more and, importantly, that it grew fast and large enough so that seasonally resolved records across multiple years were obtainable from it. Clumped isotope analysis revealed a mean annual temperature of 13.5 ± 3.8 ^∘C. The subspecies likely experienced slower growth during winter and thus may not have recorded temperatures year-round. This reconstructed mean annual temperature is 3.5 ^∘C warmer than the pre-industrial North Sea and in line with proxy and modelling data for this stratigraphic interval, further solidifying A. benedeni benedeni's use as a climate recorder. Our exploratory study thus reveals that Angulus benedeni benedeni fossils are indeed excellent climate archives, holding the potential to provide insight into the seasonality of several major climate events of the past ∼ 25 million years in northwestern Europe.

1 Introduction

Bivalves can be powerful climate archives. They add growth increments to their shells on timescales as short as sub-daily timescales (Schöne et al., 2002) and often live for many years – up to 500 years in certain species (Butler et al., 2013). Furthermore, most bivalve species precipitate their shells isotopically in near equilibrium with water, and their isotopic composition is little affected by so-called “vital effects” that skew temperature calibrations (Weiner and Dove, 2003; Huyghe et al., 2022). They therefore have the potential to record environmental changes at a seasonal resolution. This is a higher temporal resolution than can be obtained from, for example, deep-sea benthic oxygen isotope curves (Westerhold et al., 2020), and resolving seasonality changes of the past is vital to our understanding of modern climate change: climate models predict an enhanced seasonal contrast in both temperature and precipitation (Seneviratne et al., 2021). Short-term events such as heatwaves and extreme precipitation have already been observed to have increased in frequency and are predicted to become even more common in the future (Seneviratne et al., 2021).

The tellinid bivalve genus Angulus has not previously been studied for the purpose of climate reconstruction, but it holds the potential to be a valuable climate archive. Seasonality studies have previously been carried out successfully on other tellinids, such as Donax variabilis (Jones et al., 2005) and Donax obesulus (Warner et al., 2022). Here, we study a representative subspecies of this genus, Angulus benedeni benedeni (Nyst and Westendorp, 1839), in order to assess its suitability as a climate archive. This subspecies is the final taxon of an endemic North Sea basin lineage that occurred from the late Oligocene to the Late Pliocene (Marquet, 2004, 2005). The lineage therefore spans several major climatic events such as the Oligocene–Miocene transition, the Middle Miocene Climatic Optimum, and the mid-Piacenzian Warm Period (mPWP: Westerhold et al., 2020). The generic attribution of the lineage of A. benedeni benedeni is still the subject of research; it may be more appropriately referred to as Tellina. The investigated specimens either originate from the mPWP (3.264–3.025 Ma; Dowsett et al., 2010) of the southern North Sea basin, Belgium, or slightly post-date this period. During the mPWP, global average temperatures were 2–4 ^∘C higher than today (Haywood and Valdes, 2004; Dowsett et al., 2012; Haywood et al., 2020), and atmospheric CO₂ concentrations were 340–380 ppm (Pagani et al., 2010; de la Vega et al., 2020). While current atmospheric CO₂ concentrations have already exceeded 400 ppm (Monastersky, 2013), mPWP temperatures are similar to what has been predicted for moderate IPCC pathways in 2100 (IPCC, 2021).

In order to extract seasonal temperature records from a bivalve shell, a high growth rate and a lifespan of several years are required. This facilitates sampling multiple data points per year over multiple seasons. Meeting this requirement in A. benedeni benedeni was achieved through oxygen isotope analysis of shells sampled at millimetre resolution. Subsequent temperature reconstruction was carried out through carbonate clumped isotope analysis in conjunction with oxygen isotope analysis. Clumped isotope analysis allows for temperature reconstructions without knowledge of the oxygen composition of the seawater (δ¹⁸O_sw), which is required for the more traditional palaeothermometer based only on δ¹⁸O of carbonates (δ¹⁸O_c) (Ghosh et al., 2006; Schauble et al., 2006). This independence from δ¹⁸O_sw is important in estimating temperature for the mPWP, as the size of the ice sheets – which store light ¹⁶O and thus drive up the δ¹⁸O of the oceans – is poorly constrained for this period (Dowsett et al., 2016). It is also relevant for the shallow–marine waters in which bivalves often live, as the δ¹⁸O_sw can be influenced by changing influxes of isotopically light river water (Schöne et al., 2004; Chauvaud et al., 2005; Johnson et al., 2009). For both high-resolution sampling and stable isotope analysis, preservation is key. Stable-isotope-based temperature reconstructions necessitate the preservation of the original isotopic signal. Therefore, we applied a range of methods for analysing shell mineralogy, microstructure, and chemical composition to assess how best to screen for shell preservation in Angulus. In addition to screening for widespread diagenetic alteration using micro-X-ray fluorescence (XRF) and X-ray diffraction (XRD), a detailed study on potential diagenetic alteration of A. benedeni benedeni's shell was carried out using electron backscatter diffraction (EBSD). This method has as an additional advantage that it provides an overview of the shell's microstructures. Our results show that A. benedeni benedeni's shell records well-preserved multiannual records that can be sampled at a seasonal scale, making this subspecies a promising candidate for further palaeoclimatic study.

2 Geological setting

Three specimens of Angulus benedeni benedeni (SG-125, SG-126, and SG-127) were obtained from the Pliocene sequence of Belgium in the southern North Sea basin. During this period, the North Sea basin was only open to the Atlantic Ocean on the north side. Modern-day France and England were connected by the Weald–Artois land bridge to the south (Gibbard and Lewin, 2016). Proxy data indicate strong warming in the North Atlantic region during the mid-Piacenzian, with global mean sea-surface temperature (SST) anomalies of up to +4–7 ^∘C instead of the global +2–4 ^∘C (Dowsett et al., 2012, 2013). The specimens were collected in 2013 from a shell bed near the top of the Oorderen Member of the Pliocene Lillo Formation. The exact location was the temporary outcrop for the construction of the Deurganck Dock Lock (Deurganckdoksluis), now Kieldrecht Lock (Kieldrechtsluis; 51^∘16^′44^′′ N, 4^∘14^′52^′′ E), in the port of the Antwerp area (Fig. 1a, b). Photos of the collecting locality with an in situ specimen of A. benedeni benedeni can be found in Fig. 6 of Deckers et al. (2020). The marine Lillo Formation (De Meuter and Laga, 1976) was deposited during the late Zanclean to the Piacenzian (ca. 3.7–2.6 Ma) in the southern North Sea basin (Fig. 1c). The Oorderen Member is an unconsolidated, greyish fine-grained shell-rich unit containing numerous mollusc shells. Its burial depth likely did not exceed 100 m (Johnson et al., 2022). Angulus benedeni benedeni is common in upper bioturbated clayey interval. The Oorderen Member was deposited in a neritic environment (Louwye et al., 2004; De Schepper et al., 2009). Marquet (2004) estimated a water depth of 35–45 m based on bivalve depth ranges but did not differentiate between the lower sandy (Atrina and Cultellus) and the upper clayey (benedeni) levels, which all differ in their lithology and sedimentary structures. Gaemers (1975) estimated a minimum depth of 10–20 m for the Oorderen Member, based on molluscs, otoliths, and sedimentary structures. The upper clayey level must have had the shallowest palaeodepth at the time of deposition. The Oorderen Member has previously been interpreted as having been deposited during the mPWP (Louwye and De Schepper, 2010; Valentine et al., 2011), defined as the interval from 3.264 to 3.025 Ma by Dowsett et al. (2010). However, despite ongoing efforts it remains complex to date these proximal sandy deposits. There are few reliable biostratigraphic dates available. Based on these, the Oorderen is older than 2.72 Myr, but younger than 3.21 Myr (De Schepper et al., 2009; Vervoenen et al., 2014; Wesselingh et al., 2020). These sparse age constraints do not exclude a slightly younger age for the top of the Oorderen Member from which the shells were collected (Fig. 1c).

Figure 1 Location of the modern North Sea basin and the stratigraphic level from which the A. benedeni benedeni specimens originated. (a) Overview map of the location of Belgium and the present-day extent of the North Sea within northwestern Europe. The dark rectangle indicates the location of map (b). Adapted after Louwye et al. (2020). (b) The geographical extent of the Lillo Formation in Belgium, after Louwye et al. (2020), and the location of the city of Antwerp. (c) Lithological column of the Lillo Formation and its members with their chronostratigraphic position, after Wesselingh et al. (2020). The Oorderen is drawn here in its longest possible range of deposition, bound by the age of the underlying Luchtbal Member (< 3.21 Ma) and biostratigraphic constraints within the Oorderen itself (> 2.72 Ma). Age boundaries for the mPWP are taken from Dowsett et al. (2010).

3 Methods

3.1 Light microscopy

Light microscopy was employed on all specimens to analyse the macrostructures of A. benedeni benedeni's shell and to count growth increments. Two specimens (SG-126, SG-127) were partially or fully embedded in Araldite 2020 epoxy, and high-polish thick (approximately 5 mm) sections were made, with a final polish using 0.3 µm aluminium oxide. These thick sections were viewed under a Keyence VHX-5000 digital microscope at ×250 and ×1000 magnification. Light microscopy composite images were made with the Fiji stitching ImageJ plugin (Preibisch et al., 2009). These composites were used to count the growth increments. The macrostructural layers were identified and named following Bieler et al. (2014), where M+ and M− denote layers external and internal relative to the pallial myostracum (M) (Fig. 2). The equivalent of these labels to another common naming scheme (Popov, 1986; Milano et al., 2017) is shown in Fig. 2d.

Figure 2 Digital microscopy images of shell layers in the A. benedeni benedeni shells. Layers are marked M+2 to M−2; growth increments are indicated with dark blue dotted lines; dog – direction of growth. (a) Section of SG-127 indicating four layers, magnification of ×250. From outer to inner shell, these are layer M+2 with crescent-shaped growth increments; layer M+1 with low-angle growth increments and prism-like crystals parallel to the layer thickness; layer M−1 with horizontal, parallel growth increments; and layer M−2, dark brown with horizontal, parallel growth increments (difficult to see in image a due to a similar tone to the background). (b) Section of SG-126 showing a part of layer M+2 in more detail, magnification of ×500. (c) Section of SG-126 showing a part of layer M−2 at a better contrast, magnification of ×250. (d) Schematic overview of the four shell layers and their positions relative to each other and the pallial myostracum (M). In brackets is the equivalent of the layer name in the terminology of Popov (1986, 2014).

3.2 X-ray diffraction

To check for diagenetic remineralisation of original aragonite into secondary calcite, powder X-ray diffraction (XRD) analysis of a representative section of one A. benedeni benedeni specimen (SG-125) was carried out on a Bruker D8 Advance diffractometer at Utrecht University. The A. benedeni benedeni pattern was compared to XRD spectra of pure aragonite and calcite samples. The calcite spectrum was recovered from the RRUFF database (Lafuente et al., 2015; RRUFF ID: R040070). For the pure aragonite spectrum, 1.12 g of shell from a modern common cockle (Cerastoderma edule) collected at the North Sea coast near Hoek van Holland (51^∘59^′14^′′ N, 4^∘05^′56^′′ E) was prepared and measured under the same conditions as the fossil samples. Further measurement details are found in Appendix A1 in Table A1.

3.3 Micro-X-ray fluorescence

For additional screening of diagenetic alteration, concentrations of major and trace elements (including Fe, Mn and Sr) in the three specimens of A. benedeni benedeni (SG-125, SG-126, SG-127) were determined using non-destructive, energy-dispersive micro-X-ray fluorescence (μXRF) analysis. A benchtop Bruker M4 Tornado μXRF scanner was used for this, available at the Analytical, Environmental, and Geochemistry Research Group (AMGC) of the Vrije Universiteit Brussel (Brussels, Belgium). Online calibration (initial quantification) of the spectra was done with the matrix-matched standard BAS-CRM393 (Bureau of Analysed Samples, Middlesbrough, UK) in the M4 software (de Winter et al., 2017a). Following this, an offline calibration was carried out using a set of matrix-matched carbonate standards (Vellekoop et al., 2022). Further details are found in appendix A2. Several line scans were made per shell, after which elemental concentrations were calculated for the entire shell.

3.4 Electron backscatter diffraction

To map potential early diagenesis that was not detected in bulk XRD or μXRF analysis, as well as to determine which microstructures are present in A. benedeni benedeni, electron backscatter diffraction (EBSD) was carried out on thin sections (ca. 40 µm) of all three specimens. EBSD involves the diffraction of electrons, generated by a SEM-based electron beam, off the crystal planes of a sample, forming a diffraction pattern related to the crystal lattice parameters. This allows for the precise determination of crystal orientation, shape, and size (Cusack, 2016). EBSD analysis was carried out using an Oxford Instruments Symmetry EBSD detector attached to a Zeiss Gemini 450 SEM at Utrecht University. Thin sections of the samples were mechanically polished with 0.3 µm aluminium oxide suspension and finished with chemical Syton^® polish. The thin sections were covered with a thin (several nm) carbon coating to keep charge build-up during the measurements at a minimum. Beam and map acquisition settings (voltage, beam current, dwell time, step size, high/low vacuum) were varied to obtain the optimal results for each map; the acquisition settings for the map shown here are indicated in Fig. 5. Data processing was done in the Oxford Instruments AZtecCrystal software. Data clean-up consisted of wild spike removal followed by filling in unclassified areas using an iterative nearest neighbour procedure from eight down to six nearest neighbours. The grain boundary threshold was set to 10^∘. Contoured pole figures were plotted to visualise the main crystallite orientations (Fig. 5f). For these pole figures, multiple-of-uniform-density (MUD) values were calculated. The MUD value is an indication of the degree of clustering of the poles relative to a random distribution: the higher the MUD, the stronger the co-orientation (or preferred orientation) of the crystals.

3.5 Stable isotope analysis

3.5.1 Sample preparation

Samples of approximately 100–300 µg each were taken from the outer surface of two specimens (SG-126 and SG-127) by drilling into it from the exterior along commarginal paths. Specimen SG-125 was not sampled for isotope analysis as there was too little material left following XRD analysis. A Dremel hand-held drill was used, which was equipped with a tungsten carbide drill bit of 0.8 mm diameter (Fig. B1). During sampling, care was taken not to drill too deep, as this would lead to mixing of material from different shell layers that was precipitated at different points in time (Fig. B2, see also Sect. 4.1.1). The distance from the umbo to the sample track at the axis of maximum growth was measured in planar view with digital callipers (Fig. B1, Table B1). The entire shell height of specimen SG-126 was sampled (Fig. B1a). A total of 55 samples were taken from specimen SG-126, and a total of 30 samples were taken from specimen SG-127 (Fig. B1b). The sampling resolution was 0.10 to 1.07 mm (mean: 0.48 mm, 1σ: 0.24 mm; no significant difference between the two specimens, p>0.05). Samples were not homogenised prior to analysing.

For stable isotope analysis, multiple measurements were carried out on material from each drilled path. For clarity, the following terminology is used: one “sample” refers to the total powder amount that was drilled in one sample track, which runs parallel to the growth increments. One sample therefore represents a certain amount of time. One “aliquot” represents a small amount of powder taken from a “sample” which was then analysed for stable isotopes. One “measurement” refers to the measurement of one such aliquot. Multiple measurements were taken from a single sample. “Replicates” and “duplicates” refer to different aliquots taken from the same sample.

3.5.2 Oxygen isotope mass spectrometry

Drilled samples from specimen SG-126 were first analysed for carbonate oxygen isotopes (δ¹⁸O_c) on a Thermo Scientific GasBench II gas preparation system coupled to a Thermo Scientific Delta V isotope ratio mass spectrometer; 50–100 µg of powdered material was weighed for each sample and standard. Two in-house standards, Naxos marble and Kiel carbonate, were used. Their accepted values are given in Table B2. Each sample was flushed with helium for 5 min to remove atmospheric air. The samples were reacted with 103 % phosphoric acid at 70 ^∘C. A Naxos standard was measured approximately every 10 samples during the run. From the bivalve samples, every 10th sample was weighed in duplicate. These duplicates were measured at the end of the run. The carbon isotope data (δ¹³C_c) are reported here only when it is relevant for the reproducibility of the analyses. Sample δ¹⁸O_c data were corrected for intensity-dependent fractionation using a linear regression between δ¹⁸O_NAXOS and mass 44 intensity (Fig. B3). Standard deviations and uncertainties were calculated after this correction was carried out. No significant change in the Naxos standard values over time (also called “drift”) was observed during the measurement. The δ¹⁸O_c and δ¹³C_c data are reported relative to Vienna Pee Dee Belemnite (VPDB). The internal standard deviation of conventional stable isotope analyses was 0.04 ‰ (1σ) for δ¹³C and 0.06 ‰ (1σ) for δ¹⁸O_c (averaged over all measurements, n=64). External reproducibility was 0.05 ‰ for δ¹³C and 0.07 ‰ for δ¹⁸O_c based on repeated measurements of the Naxos standard (N=11). The Kiel carbonate standard showed external reproducibility of 0.09 ‰ for δ¹³C and 0.06 ‰ for δ¹⁸O_c (N=4).

3.5.3 Clumped isotope mass spectrometry

Clumped isotope analysis on carbonate involves analysing the abundance of ¹⁸O–¹³C bonds. The deviation of the occurrence of these heavy ¹⁸O–¹³C bonds from the expected stochastic distribution is thermodynamically determined and thus independent of the isotopic composition of oxygen of the seawater (δ¹⁸O_sw) and carbon of the dissolved inorganic carbon (δ¹³C) (Ghosh et al., 2006; Schauble et al., 2006). This ¹⁸O–¹³C deviation is termed Δ₄₇. Clumped isotope values (Δ₄₇) were analysed using a Kiel IV carbonate device coupled to a Thermo Scientific MAT 253 Plus isotope ratio mass spectrometer (available at Utrecht University) with a long-integration dual-inlet (LIDI) workflow (Hu et al., 2014; Müller et al., 2017) and following the protocol described by Meckler et al. (2014). For specimen SG-126, analysis was focussed on samples with δ¹⁸O_c values in the ranges 2.0 ‰–3.3 ‰ (close to the maximum) and 0.4 ‰–1.2 ‰ (close to the minimum), based on the previously measured aliquots (Sect. 3.5.2). From specimen SG-127, all samples were analysed; 75–95 µg of powdered material was weighed out. Carbonate standards of similar weight were measured in an approximately 1:1 ratio to the standards in each run (22 samples, 24 standards; Kocken et al., 2019). ETH-3 was measured more often as it is closest in composition to the samples (Kocken et al., 2019). These standards were two control standards – Merck CaCO₃ (synthetic; product code 1.02059.0050) and IAEA-C2 (Bavarian travertine) – as well as ETH-1, ETH-2, and ETH-3 (Bernasconi et al., 2018, 2021). Their composition is given in Table B2. The ETH standards, which have varying Δ₄₇, δ¹⁸O_c, and δ¹³C compositions, were used to transfer the sample Δ₄₇ values to the ICDES-90 reference frame (Bernasconi et al., 2021). Merck and IAEA-C2 were not involved in data processing (treated as samples) but used to assess long-term measuring uncertainty.

Carbonates were reacted with 105 % phosphoric acid at 70 ^∘C. The produced CO₂ gas was purified through two cryogenic traps (−170 ^∘C) and a Porapak trap (−50 ^∘C). It then entered the mass spectrometer, and masses 44–49 were measured against a reference gas (δ¹⁸O = −4.67 ‰, δ¹³C = −2.82 ‰). To correct for negative intensity baseline in the Faraday cups of the mass spectrometer, we performed scans of the shape of the mass spectrometry peaks using reference gas measurements before each run. The baseline on masses 47, 48, and 49 was identified and corrected by regressing the size of the negative intensity outside the peak against intensity of mass 44 (He et al., 2012). The calculated δ values were then corrected for ¹⁷O using the values presented in Brand et al. (2010). The Δ₄₇ values were calculated and averaged over the 40 pulses as per the LIDI system. Correction of raw Δ₄₇ values to I-CDES values (Bernasconi et al., 2021) was done in two steps. Initially all measurements are corrected using bracketing ETH-3 measurements in order to correct for variability during single runs (i.e. over the span of several hours). The offset of these bracketing ETH-3 raw Δ₄₇ values from their I-CDES values is added to the samples. In a second step an empirical transfer function (ETF) was constructed by regressing the raw ETH-1, ETH-2, and ETH-3 values over their accepted Δ₄₇ values and using the resulting linear regression line to transfer the sample Δ₄₇ values to the I-CDES90 ^∘C reference frame (Fig. B4). For the ETF, 403 ETH standards from multiple runs measured over the span of two months were used for linear regression (53 ETH-1, 64 ETH-2, 286 ETH-3; Fig. B5).

The δ¹⁸O_c and δ¹³C data of the same measurements were corrected through a 15-point running average based on the ETH standards, in order to eliminate long-term trends present within the δ¹⁸O_c and δ¹³C raw data. The δ¹⁸O_c and δ¹³C data from three runs were not used for the running average correction nor were they used in the δ¹⁸O_c records of A. benedeni benedeni. These runs show large deviations from expected values and exhibit a strong a correlation between δ¹⁸O_c and δ¹³C (runs 466, 468, 489; these runs are marked in Fig. B5). It seems that mass-dependent fractionation occurred during the gas preparation. The Δ₄₇ results appear not to have been influenced by this fractionation as is evidenced by the measured standards, which show no systematic Δ₄₇ offset (Fig. B5). Inclusion of Δ₄₇ values from these runs did not significantly alter the temperature reconstructions beyond decreasing the confidence level on the results. Final δ¹⁸O_c and δ¹³C values were reported relative to VPDB. Samples that showed a strong drift during the 40 LIDI measurement pulses, had a low intensity, a high standard deviation, or showed signs of contamination, were deemed unreliable and were removed from the dataset. The intensity cut-off was < 9.0 V. The cut-off for standard deviations in Δ₄₇ data within the 40 measurement pulses was > 0.10 ‰ for both standards and samples. Contamination evident from increased intensity on the mass 49 cup was quantified using the 49 parameter – the ratio between mass 49 and mass 44 intensities – and a cut-off of > 0.1 ‰ was used for standards and > 0.2 ‰ for samples. After removing erroneous samples, 103 sample measurements remained. The external reproducibility for the clumped isotope analyses was 38 ppm for IAEA-C2 (n=23) and 60 ppm for Merck (n=15). The external reproducibility of IAEA-C2 was used as the standard deviation on Δ₄₇ for the samples, as it is closer in composition to the samples than Merck (Fig. B5). After corrections, most of the Δ₄₇ data ranged between 0.55 ‰ and 0.70 ‰ for both specimens.

For oxygen isotopes, the external reproducibility of the standards does not reflect the reproducibility of the samples based on duplicate measurements. The samples from each sampling track were not homogeneous, as they likely contain multiple thin growth increments and were not homogenised after drilling. Therefore, different aliquots from the same sample sometimes resulted in considerable differences when measured, in both the GasBench and the MAT 253 PLUS analyses. Instead of using the external reproducibility on the Naxos standard, the uncertainty on the shell data was determined as follows: for each sample for which replicates were analysed, the standard deviation of these replicates results was taken, and all these standard deviations were averaged. To this goal, δ¹⁸O_c data from GasBench and MAT 253 PLUS analyses were combined. This resulted in an uncertainty of 0.15 ‰ for SG-126 (n=28 replicates). There were not enough duplicate δ¹⁸O_c measurements for SG-127, so the value of 0.15 ‰ was also applied to this specimen.

3.6 Palaeotemperature reconstruction

All Δ₄₇ values were compiled and converted to a mean annual temperature (MAT) using the temperature transfer function of Meinicke et al. (2020, 2021):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\mathrm{\Delta }}_{\mathrm{47}}\left[\mathrm{‰}{\mathrm{ICDES}}_{\mathrm{90}{}^{\circ }\mathrm{C}}\right]& =\mathrm{0.0397}\pm \mathrm{0.0011}\times {\mathrm{10}}^{\mathrm{6}}\\ & \times T^{-\mathrm{2}}\left[K\right]+\mathrm{0.1518}\pm \mathrm{0.0128}.\end{array}\end{array}$$

This temperature transfer function is based on foraminifera and therefore more suited to the low temperature range in which bivalves live. Recently, it was shown that this clumped isotope thermometer yields the most accurate temperature reconstructions in aragonitic bivalves (de Winter et al., 2022). Since the temperature function for Δ₄₇ is not linear, errors associated with the Δ₄₇ data were propagated using a Monte Carlo simulation (n = 10⁵) of uncertainty in slope and intercept of the temperature transfer function as well as measurement uncertainty on Δ₄₇ values assuming a normal uncertainty distribution. The standard deviation, standard error, and 95 % confidence level were then calculated from this normally distributed simulated temperature dataset. Uncertainties on temperatures were calculated after temperature conversion. Since the temperature transfer function for clumped isotopes is not linear, this introduces a bias, specifically towards warmer temperatures. The temperature means were therefore calculated from averaging the Δ₄₇ values instead, and the Monte Carlo-generated errors were transferred to these mean temperatures.

The δ¹⁸O of the seawater (δ¹⁸O_sw) was calculated from Δ₄₇-based MAT and δ¹⁸O_c data using the temperature transfer function for aragonite of Grossman and Ku (1986), modified by Dettman et al. (1999):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{sw}}\left[\mathrm{‰}\mathrm{VSMOW}\right]& ={{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}_{\mathrm{c}}\left[\mathrm{‰}\mathrm{VPDB}\right]\\ & -{\displaystyle \frac{\mathrm{20.6}-T\left[{}^{\circ }\right]}{\mathrm{4.34}}}+\mathrm{0.2}.\end{array}\end{array}$$

VSMOW stands for Vienna Standard Mean Ocean Water. This average δ¹⁸O_sw, the Δ₄₇-based MAT, and minimum and maximum δ¹⁸O_c values were then re-inserted into this equation to obtain δ¹⁸O_c-based summer and winter temperatures.

3.7 Growth models

To gain insights into the nature of the growth of A. benedeni benedeni, two types of growth models were fitted to the growth data of specimen SG-126: the Von Bertalanffy growth function (VBGF), a logarithmic model, and the Gompertz equation, a logistic model. Both are widely used for estimating and analysing growth in modern species (Lee et al., 2020). The built-in nls function in R (R Core Team, 2021) was used to fit the data to both growth functions (see code availability). The VBGF and the Gompertz equations are restrictive in terms of how much their shape can change, and more sophisticated approaches are available today (Lee et al., 2020). However, as this study has only one specimen available for growth model fitting, and since it concerns a fossil species, we believe that this simple approach is appropriate. The VBGF has the form of Eq. (3), where H is the body size of the specimen at time t (here: shell height measured as distance from umbo), H_asymp is the asymptotic height (i.e. the theoretical maximum shell height), k is the growth coefficient, t is the time (here: in years), and t0 is the theoretical time where H equals 0.

$$\begin{array}{}\text{(3)}& H\left(t\right)=H_{\mathrm{asymp}}\left(\mathrm{1}-e^{-k\left(t-t_{\mathrm{0}}\right)}\right)\end{array}$$

The Gompertz equation has the form of Eq. (4):

$$\begin{array}{}\text{(4)}& H\left(t\right)=H_{\mathrm{asymp}}{e}^{-{\mathrm{be}}^{-\mathrm{ct}}},\end{array}$$

where H is again the shell height at time t, H_asymp is the theoretical maximum shell height, and b and c are coefficients. Each parameter, its associated standard error (SE), and p value were calculated.

4 Results

4.1 Light microscopy

4.1.1 Structural analysis

Angulus benedeni benedeni shells consist of four layers, M+2 to M−2 (Fig. 2), labelled from the outer to the inner surface. M+2 and M+1 are located outside of the pallial myostracum and form the outer shell layers. M+2 is a thin (ca. 100 µm) white-to-grey coloured layer with crescent-shaped growth increments. M+1 is a ca. 200 µm thick layer (measured at the thickest point, halfway along the shell height) with low-angle growth increments and white–beige in colour. M−1 and M−2 are located inside of the pallial myostracum and form the inner shell layers. M−1 is a thick (ca. 1000 µm) layer with clear growth increments parallel to the layer boundaries and a white–beige colour. M−2 is a thin (ca. 50 µm) brown-coloured layer with horizontal, parallel growth increments similar to M−1.

4.1.2 Growth increments

In specimen SG-126, a maximum of 115 growth lines were counted in layer M−1 at the hinge (Fig. 3a, c). In specimen SG-127, a maximum of 130 growth lines were counted just ventral of the hinge (Fig. 3b, c), also in layer M−1. The growth lines were not all clearly visible, and in some areas in the shell, fewer growth lines could be discerned. As it was deemed easier to miss some growth lines due to the limited resolution than to overestimate the number, the microscope images with the maximum counted number of growth lines are shown here. Although still an underestimate, these are thought to be the best estimate of the actual number of growth lines.

Figure 3 Growth lines in the A. benedeni benedeni specimens. (a) Counted growth lines in specimen SG-126, 115 in total. Counted in shell layer M−1. (b) Counted growth lines in specimen SG-127, 130 in total. Counted in layer M−1. (c) Schematic figure of the shell indicating where the growth lines were counted.

4.2 X-ray based analyses

X-ray diffraction analysis (XRD) and micro-X-ray fluorescence (μXRF) point analysis indicate that the aragonite in the shells of A. benedeni benedeni was not diagenetically altered. XRD analysis revealed that the shell of A. benedeni benedeni consists of 100 % original aragonite, as its spectrum is identical to that of pure aragonite (Fig. 4a). Micro-X-ray fluorescence (μXRF) point analysis supports this interpretation, as the specimen is high in Sr while low in Fe and Mn (Fig. 4b). Higher concentrations of Fe and Mn and lower concentrations of Sr in fossil carbonates are generally associated with pore water interaction during diagenetic alteration (e.g. partial recrystallisation) of the carbonate (Brand and Veizer, 1980; Al-Aasm and Veizer, 1986; Hendry et al., 1995). The absence of high Fe and Mn concentrations and presence of high Sr suggest that this process did not take place in these specimens.

Figure 4 X-ray-based analysis results to assess diagenetic alteration in the specimens. (a) X-ray diffraction analysis of A. benedeni benedeni compared with patterns of pure aragonite and calcite, indicating an aragonitic composition of the shell. (b) Violin plot showing the results of the X-ray fluorescence analysis of A. benedeni benedeni for the elements Fe, Mn, and Sr. The “violin” shape depicts the kernel density function for the elemental concentration, plotted vertically and mirrored. The specimen is high in Sr and shows a broad distribution within these higher values. While there are some data points with high Fe and Mn values, the distribution is strongly skewed towards low values, suggesting overall low Fe and Mn values.

4.3 Electron backscatter diffraction

4.3.1 Diagenesis assessment

Electron backscatter diffraction (EBSD) analysis shows no evidence of diagenetic alteration in the shells. In the EBSD maps, 0 % to 0.2 % of the area was initially classified as calcite. Upon inspection, these “calcite” grains are wild spikes that were subsequently removed in the data clean-up (Fig. C1). There was no local secondary mineral growth visible in any map (e.g. Fig. 5b–e), which would have presented itself as larger crystals (“blocky calcite”) that do not follow the surrounding structures (Cusack, 2016; Casella et al., 2017). The pole figures show a clear preferred orientation (Fig. 5f), which is indicative of original growth structures (Casella et al., 2017).

Figure 5 EBSD map of specimen SG-127. Map acquisition specifics: 20 kV, high vacuum, 0.2 µm step size. (a) Band contrast map with scale, orientation, growth lines, direction of growth (dog), and shell layers M+2, M+1, and M−1, separated by dashed lines. (b) EBSD inverse pole figure (IPF) X map, showing grain orientations parallel to the X0 axis. Unit cells of aragonite with the most common orientations are shown, the colour of which corresponds to the colour coding of the map. (c) EBSD IPF Y map showing grain orientations parallel to the Y0 axis, with unit cells. (d) EBSD IPF Z map showing grain orientations parallel to the Z0 axis, with unit cells. (e) Inverse pole figure (IPF) colour key for the EBSD orientation maps. (f) Contoured pole figures for the three main crystallographic directions and the MUD, for the different shell layers. (g) Grain size distribution chart for the entire map and the different shell layers. Note the logarithmic y axis. (h) Schematic figure indicating the approximate location of the map on the shell. (i) Shell layers from map IPF Z zoomed in and with added grain (black) and twinning (yellow) boundaries.

4.3.2 Grain size and orientation

The aragonitic shell material is fine grained, and its crystals show a strong preferred orientation as well as twinning. Most grains have an area of a few µm² or smaller (Fig. 5g). Layer M−1 has a narrower grain size distribution than M+1 and M+2, with no grains larger than 25 µm² present in the analysed section. The multiples-of-uniform-density (MUD) values indicate a strong preferred orientation (Fig. 5f). Layer M−1 has an increased MUD of 104.07 compared to ca. 40 for the other two layers, indicating a stronger preferred orientation (Fig. 5f). This is due to the strong preferred orientation of the 100 axis, while the 001 and 010 axes show a weaker preferred orientation compared to the other layers. The pole plots show double maxima of the 001 and 010 axes, which rotate around the 100 axes (M+2 and M+1, Fig. 5f). This is observed in layers M+2 and M+1. In layer M−1, the 001 and 010 axes form a girdle around the 100 axis rather than two distinct maxima (M−1, Fig. 5f). The rotational 100 axis is generally oriented parallel to the growth direction. This is especially evident from the lamellae in layer M−1 of Fig. 5b–e.

The angle of ca. 64^∘ between the 001 and 010 axes is characteristic of aragonite polysynthetic (110) twinning (Griesshaber et al., 2013). Calculated twinning boundaries indicate that the two dominant orientations, as indicated by the inverse pole figure (IPF) colour coding, in layers M+2 and M−1 represent twins (Fig. 5i). In layer M+1, the two dominant orientations that make up the lamellae are not a result of twinning. Instead, each of the two orientations (orange and purple) consists of their own set of twins (light and dark orange and purple, respectively; Fig. 5i).

4.3.3 Microstructures

Within the EBSD maps, several different structures are observed (Fig. 6). Layer M+2 consists of a compound composite prismatic structure, formed by bundles of crystals that diverge from the centre of the layer (Fig. 6a) (Popov, 2014). Layer M+1 consists of a crossed lamellar structure. This is indicated by two sets of crystal orientations that are oriented perpendicular to the growth lines and that alternate and interfinger, similar to a zebra pattern (Fig. 6b) (Popov, 2014; Crippa et al., 2020). Layer M−1 shows multiple structures. In the hinge region, a hatched pattern is visible (Fig. 6c), corresponding to a complex crossed lamellar structure (Popov, 2014). In other parts of M−1, elongated, columnar crystals with their long side oriented parallel to the growth direction are observed (Fig. 6d), indicative of a prismatic structure (Popov, 2014). Around this prismatic structure, there are small grains with no clear preferred shape orientation (Fig. 6d). This has been interpreted as a homogenous structure (Popov, 2014). However, it might also be a different structure that cannot be analysed at this resolution due to the fine-grained nature of the aragonite. Layer M−2 was not captured on any EBSD map.

Figure 6 Different shell structures identified in the EBSD maps of specimen SG-127. (a) Compound composite prismatic structure in layer M+2. (b) Crossed lamellar structure in layer M+1. (c) Complex crossed lamellar structure in layer M−1 at the hinge of the shell. (d) Prismatic structure in layer M−1. There are also areas in this layer (e.g. top right in panel d) without a clear grain shape structure, which are interpreted as a homogeneous structure. (e) Schematic figure of the shell indicating the approximate location of the maps.

4.4 Stable isotope analysis

4.4.1 Oxygen isotopes

The δ¹⁸O_c data of specimens reveal sinusoidal records that range between +0.3 ‰ and +3.3 ‰ (Fig. 7a, b) and suggest multi-annual growth, with δ¹⁸O_c peaks corresponding to winters and valleys corresponding to summers. Specimen SG-126 shows at least 8–9 years of growth (Fig. 7a). The first and last years of this specimen may not have been recovered during sampling, as the first and last few millimetres were not sampled. In addition, the sampling resolution was likely too coarse to capture all years in the ontogenetically oldest section of the shell, as growth increments generally become thinner with age. In SG-126, high δ¹⁸O_c values generally correspond to darker bands and vice versa, which suggests that the dark–light couplets represent the winter and summer seasons. The δ¹⁸O_c variability decreases near the ventral margin of the shell. The sampled section of SG-127 shows six cycles in δ¹⁸O_c values (Fig. 7b). In contrast to SG-126, the relationship between δ¹⁸O_c value and shell colour is not consistent. There are many thin light and dark bands on SG-127, and samples likely contain material from both light and dark bands (Fig. 7b; not all indicated on Fig. 7b as the exact location relative to the samples was not determined).

Figure 7 δ¹⁸O_c records of the A. benedeni benedeni shells, plotted relative to the distance from the umbo. (a) Record of specimen SG-126, consisting of 55 samples each measured 1 to 3 times. (b) Record of specimen SG-127, consisting of 30 samples. The average was calculated from the different measurements of a single sample. If there was just one measurement, that value was taken; 1 and 2 standard deviations are plotted around the average. The standard deviation was calculated as described in Sect. 4.4.1 (0.15 ‰). In the background, the approximate colour of the growth band on the outer surface of the shell is indicated (see also Fig. B1). The warm and cold δ¹⁸O_c “bins” are marked in red and blue bands. These bins encompass all data points below 0.9 ‰ and above 2.4 ‰. See also Sect. 3.4.3 and Fig. 9. Raw data are plotted in Fig. B6.

4.4.2 Growth curves

A growth curve was reconstructed for the specimen SG-126 using the δ¹⁸O_c record (Fig. 8a–b). This was not done for specimen SG-127 as only a small section of the entire shell was sampled. As there appear to be several smaller fluctuations in the δ¹⁸O_c record, a “year” was defined as a fluctuation that consists of 2 or more data points. This leaves out a small peak at 9.97 mm from the umbo and a larger peak at 21.02 mm (Fig. 8). An alternative growth model that counts this larger peak as a year is shown in Fig. B7. Figure 8 shows the fitted Von Bertalanffy and Gompertz growth models. For the Von Bertalanffy model, the parameters H_asymp, k, and t0 were estimated to be 52.35 ± 7.64 SE (p<0.05), 0.09 ± 0.02 SE (p<0.05), and 0.03 ± 0.16 SE (p≫0.05) respectively. H_asymp is the estimated maximum shell height in millimetres that A. benedeni benedeni would have reached, ca. 52 mm. The interpretation of growth coefficient k is problematic, and besides having a large standard error and low confidence, t0 is merely a mathematical artefact that does not have any biological meaning (Lee et al., 2020). For the Gompertz model, parameters H_asymp, b, and c were estimated to be 34.11 ± 1.98 SE (p<0.05), 2.58 ± 0.14 SE (p<0.05), and 0.29 ± 0.03 SE (p<0.05), respectively. The Gompertz model thus suggests a shorter maximum shell height of 34.11 mm. The standard error on the residuals is 0.51 for the VBGF and 0.63 for the Gompertz equation. The R² value for both fits is 0.99, but as both functions are non-linear, this is not a reliable indicator of goodness of fit (Spiess and Neumeyer, 2010). This rather preliminary growth model fitting suggests a maximum size of A. benedeni benedeni of 34–52 mm.

Figure 8 Growth curves for specimen SG-126. (a) Growth curves constructed by inserting the counted years versus distance-from-umbo data points into two different fitting algorithms: the asymptotic Von Bertalanffy growth function and the logistic Gompertz function. (b) δ¹⁸O_c record with the inferred growth years marked by lines connecting to plot (a).

4.4.3 Δ₄₇ and palaeotemperature

A mean annual temperature (MAT) of 13.5 ± 3.8 ^∘C has been calculated through averaging all Δ₄₇ values and converting this average to a single temperature, with the 95 % confidence level (CL) determined through Monte Carlo error propagation. Reliable Δ₄₇-based summer and winter temperatures could not be determined due to a limited amount of data. This issue is illustrated in Fig. 9. To obtain summer and winter datasets, we selected samples with respectively low and high δ¹⁸O_c values from both specimens and subsequently averaged their associated Δ₄₇ values (similar to the “optimisation” approach as discussed in de Winter et al., 2021a). This grouping by δ¹⁸O_c value is necessitated by the large uncertainty on individual Δ₄₇ measurements, which means that Δ₄₇ data points cannot be split into meaningful “warm” and “cold” groups directly. The choice of a cut-off value for warm and cold δ¹⁸O_c values is a trade-off between signal and noise: a cut-off close to the average δ¹⁸O_c value results in a higher sample size and a narrower confidence interval for the resulting temperature reconstruction. However, this also results in averaging samples from a large part of the year and thus a dampening of the seasonality signal. For a cut-off that only includes the lowest and highest δ¹⁸O_c values, the opposite is true. To determine whether there was an acceptable compromise between signal and noise, the average warm and cold Δ₄₇ values and their associated 95 % confidence levels were plotted for a range of δ¹⁸O_c-based cut-offs (Fig. 9). Oxygen isotope data from both specimens were combined to improve the statistics. Combining δ¹⁸O_c and Δ₄₇ data to reconstruct seasonality is the recommended strategy when a sampling resolution of 1 data point per month cannot be obtained, as was the case on the small shells of A. benedeni benedeni (Zhang and Petersen, 2023).

Figure 9 Overview of averages (round and diamond symbols) and 95 % confidence intervals (dashed lines) for different-sized groups of summer and winter Δ₄₇ aliquots. Sizes of groups increase towards the middle of the plot. The larger symbols with vertical error bars highlight the issue of having insufficient aliquot measurements. The bin sizes are small, and so they should come close to the true (non-averaged) seasonality. However, the number of aliquots is reduced to ca. 20 at this point. This translates to large errors, making these cold and warm averages statistically indistinguishable from each other.

Figure 9 shows that an acceptable compromise between signal and noise to resolve seasonality is not possible with the limited amount of data obtained here. The δ¹⁸O_c cut-offs for < 0.9 ‰ and > 2.4 ‰ are highlighted to illustrate this. At these cut-off points, which are arbitrarily chosen, the number of data points (n=20 and n=21, respectively) is large enough to reduce the 95 % CL range. Further improvement within this dataset is hardly possible as the error bars plateau shortly after, corresponding to ca. n=30. However, the two datasets are not statistically different (p>0.05, Student's t test) already in sample sizes starting from 20 aliquots, due to both the wide δ¹⁸O_c range included in these cut-offs and the large 95 % CL range. Increasing the number of clumped isotope measurements that correspond to the more extreme δ¹⁸O_c values would decrease the large errors and would allow for making the cut-offs narrower. This way, they better represent the seasonal extremes. That this approach works is exemplified by the MAT of 13.5 ± 3.8 ^∘C. Due to the large number of measurements (N=103), the error on this temperature has been reduced to a few degrees Celsius. The same is possible for seasonal temperatures if the dataset is large enough. Besides the large analytical error associated with clumped isotopes, seasonal variability will also add to the standard deviation on the MAT. Even with more data, the degree to which this error can be reduced is therefore limited. Clumped-isotope-based summer and winter temperatures are not calculated here and will not be discussed further, as the data generated in this study are not sufficient to do so in a meaningful way. Instead, we discuss an alternative approach to obtain seasonal temperatures below.

4.4.4 δ¹⁸O_sw, δ¹⁸O_c, and temperature

The average δ¹⁸O_sw is 0.10 ± 0.88 ‰ VSMOW (95 % CL), based on the average Δ₄₇-based temperature of 13.5 ± 3.8 ^∘C and the average δ¹⁸O_c of both shells combined (1.53 ± 0.015 VSMOW, 95 % CL). The Grossman and Ku (1986) equation was applied to the highest and lowest δ¹⁸O_c values and the annual mean δ¹⁸O_sw to obtain winter and summer temperatures (Fig. 10a). To obtain these highest and lowest δ¹⁸O_c values, the three lowest (0.30 ‰, 0.42 ‰, 0.43 ‰) and highest (3.07 ‰, 3.16 ‰, 3.29 ‰) data points were averaged. These three points were chosen to capture the largest seasonal range that was recorded while not relying on a single measurement. This resulted in summer and winter temperatures of 18.5 ± 3.9 and 6.4 ± 3.9 ^∘C (95 % CL, Fig. 10b). Note that this approach does not take into account seasonal changes in δ¹⁸O_sw.

Figure 10 δ¹⁸O_c-based summer and winter temperatures. (a) Schematic illustrating how the δ¹⁸O_c-based temperatures have been calculated and what values were used. All values after ± are 95 % CLs. (b) Mean δ¹⁸O_c-based summer (red circle) and winter (blue diamond) temperatures and their 95 % CL. In the background, the Δ₄₇-based mean temperature and confidence interval is shown (grey band). The δ¹⁸O_c values used are the average of the lowest and highest three data points.

5 Discussion

5.1 Diagenesis

No signs of diagenesis of the shells were found through micro-X-ray fluorescence (μXRF), X-ray diffraction (XRD), and electron backscatter diffraction (EBSD) analyses. The μXRF analyses indicated that the shell material was high in Sr and low in Fe and Mn; these latter two elements are regularly used to pinpoint diagenesis (Fig. 2b; de Winter and Claeys, 2017). The XRD pattern of specimen SG-125 consisted of 100 % aragonite. No blocky calcite was observed in the EBSD images of any of the three specimens (Cusack, 2016; Casella et al., 2017). This alone does not preclude diagenetic alteration, as calcite is not formed in the earliest stages of diagenesis (Cochran et al., 2010; Marcano et al., 2015; Ritter et al., 2017). Further analysis of the EBSD maps did not, however, produce any evidence for diagenetic alteration. Microstructures that are similar to those found in modern Tellinidae were observed in the EBSD maps (Popov, 2014). In altered aragonitic bivalves, the structures as observed from EBSD maps are more homogenous and have a lower MUD compared to pristine material (Casella et al., 2017; study of Arctica islandica). Even the homogeneous structure in A. benedeni benedeni (Fig. 6d) appears less chaotic than altered aragonite (Casella et al., 2017). Combined with the shallow burial depth (≤ 100 m; Johnson et al., 2022) and the unconsolidated and fine-grained nature of the surrounding sediment, this evidence suggests that no significant diagenetic alteration took place. The isotope values presented here therefore record the formation temperature of the biogenic aragonite, assuming equilibrium fractionation, and support the interpretation of the stable isotope data in this study.

Previous studies have indicated good preservation of bivalves from the Oorderen Member (Valentine et al., 2011; Johnson et al., 2022), from other members of the Lillo Formation (Johnson et al., 2022), and from the Early Pliocene Ramsholt Member of the Coralline Crag Formation (Johnson et al., 2009; Vignols et al., 2019). Pliocene shells from the North Sea basin are expected to yield reliable temperature results via geochemical analysis, as long as preservation is not obviously suspect (e.g. in the Sudbourne Member of the Coralline Crag Formation in southeastern England, whose aragonitic shells have nearly all been dissolved; Balson et al., 1993). However, it remains imperative to screen for early diagenesis when analysing shells from new localities and older stratigraphic intervals. This is especially relevant for seasonality studies. During early diagenesis, the intra-shell δ¹⁸O_c variations that reflect seasonal variability can be homogenised, dampening the original signal. Evidence that this has taken place is limited to subtle changes in scanning electron microscope images (Moon et al., 2021). To the authors' knowledge, whether these changes are detectable with EBSD has not been studied yet.

Based on the results in this study as well as practical considerations, we recommend a combination of light microscopy and trace element analysis (XRF) for diagenetic screening in general, while the addition of XRD analysis is advised only for aragonitic shells (such as A. benedeni benedeni). While EBSD is a powerful method for observing and classifying shell microstructures and assessing their preservation, its limited accessibility and comparatively time-consuming data processing renders this method somewhat unwieldy for routine preservation screening.

5.2 Angulus benedeni benedeni shell structure

The shell of A. benedeni benedeni consists of four macrolayers, of which the outer three show microstructures that were identified using EBSD. From outer (M+2) to inner (M−1), these are compound composite prismatic, crossed lamellar, complex crossed lamellar–prismatic–homogeneous (Popov, 2014). These structures have all been described previously in other tellinid genera such as Macoma, Peronidia, and the lucinid Megaxinus (Popov, 2014). Tellinidae generally have three layers, sometimes with sub-layers present. The M−1 layer in Tellinidae is usually homogeneous, as it appears to be in our specimens (Popov, 2014). Beyond this, there exists a wide range of variation in tellinid bivalve shell microstructures (Popov, 2014). Angulus benedeni benedeni is similar in structure to the closely related A. nysti (Popov, 2014). Layer M−2 was not studied with EBSD here, but light microscopy suggests that it might be structurally similar to the neighbouring M−1 layer (Fig. 2).

Aragonite polysynthetic (110) twinning – indicated by a 64^∘ misorientation angle between the 001 and 010 axes – is frequently observed in bivalves (Kobayashi and Akai, 1994; Griesshaber et al., 2013; Crippa et al., 2020) and other molluscs (Schoeppler et al., 2019). It contributes to rapid shell growth as it is more efficient at filling up space than non-twinned growth (Schoeppler et al., 2019). This twinning is observed in all analysed layers of A. benedeni benedeni. As the double 001 and 010 maxima are stronger in the external layers M+2 and M+1 than in the internal layer M−1, twinning may be more significant in these two layers. This may, in turn, be linked to more rapid precipitation. This study supports previous investigations in highlighting the utility of EBSD for analysing shell structures in bivalves (Cusack, 2016; Checa et al., 2019; Crippa et al., 2020), which remains a rather new application of this technique.

5.3 Angulus benedeni benedeni growth history

Angulus benedeni benedeni experienced slower growth in winter, could live for up to a decade or longer (Figs. 7, 8), and likely reached lengths of 34–52 mm. The maximum (winter) δ¹⁸O_c values decrease, while the minimum (summer) values stay similar throughout the oxygen isotope record of specimen SG-126. This amplitude reduction may reflect growth breaks during the colder months. This is supported by the correlation of high δ¹⁸O_c values with dark growth bands on the shell's exterior (Fig. 7). Darker bands in bivalves often have high organic matter content related to decreased mineralisation rates (Lutz and Rhoads, 1980, cited in Carré et al., 2005). The pronounced ontogenetic growth rate decline that is often seen in bivalves (McConnaughey and Gillikin, 2008) is not apparent from the growth curve of specimen SG-126. The growth curve only becomes slightly less steep in the last few measured years, and there is no marked decrease in the wavelength of the δ¹⁸O_c record (Fig. 8).

Simple growth modelling has provided a preliminary range of maximum shell height for A. benedeni benedeni. The VBGF yields a higher maximum shell height estimate than the Gompertz equation (52 vs. 34 mm, respectively). Due to the nature of the respective models, the VBGF tends to overestimate the maximum shell height, while the Gompertz equation tends to underestimate it (Urban, 2002). These two estimates therefore yield a plausible range for the maximum shell height of A. benedeni benedeni. The heights of specimens SG-126 and SG-127 (30–36 mm) are close to the H_asymp of the Gompertz model, and the other specimens in the collection of the Royal Belgian Institute of Natural Sciences generally do not exceed 40 mm in height. The shape of the VBGF and Gompertz curves suggest that specimen SG-126 had not yet reached its maximum shell height. A. benedeni benedeni may have lived longer than 8–9 years.

Bivalves can form growth increments on a variety of timescales, and their growth rhythm can be influenced by solar and lunar cycles (Tran et al., 2011). The number of growth lines counted in A. benedeni benedeni was divided by the inferred age of the specimen to determine its growth rhythm. The sampled intervals span several years: around 9 years for SG-126, which represents most of the shell height, and around 6 years for SG-127, which represents only around one-third of the shell height. In both shells, 100+ growth lines were counted. In SG-126, the maximum number of years counted (9) and the maximum number of growth lines counted (115) correspond to ca. 13 growth lines per year. This is an interval of approximately 29 days, which matches the periodicity of the monthly synodic lunar–tidal cycle. However, both growth lines and years are likely to be underestimated. Their recordings are both limited by either the microscopy or the sampling resolution. Additionally, the shells likely experienced winter growth stops, so the growth increments do not record a full year. It is difficult to put constraints on these uncertainties due to limited data availability. These uncertainties work in different directions (underestimating the number of years results in more time represented within a growth increment; underestimating the number of growth increments and the amount of time missing due to winter stops results in less time within a growth increment). As two out of three factors push the duration of a growth increment towards lower values, each increment may represent less than 29 days. This fits with the nature of lunar–tidal forcing. There is little difference between the two spring–neap cycles within a monthly cycle (Kvale, 2006). If tides exerted influence, the fortnightly spring–neap cycle is expected to be present as well due to its larger amplitude. Monthly lunar cycles have been observed in bivalves (Pannella and MacClintock, 1968), but these are present as bundled spring–neap tide pairs with one being more dominant than the other, not as individual growth lines as observed here. Actual monthly forcing of environmental conditions linked to growth increments is therefore unlikely.

It is also possible that the formation of growth lines is driven by aperiodic or quasi-periodic processes such as storms – which would skew the correlation between growth years and growth lines – or is a result of an internally controlled rhythm rather than being externally forced. Internal forcing is a strong possibility. Finally, many bivalves form growth increments on smaller timescales as well (e.g. (semi-)diurnal or circatidal; Judd et al., 2018, and references therein). We cannot rule out the presence of such ultradian cyclicity in A. benedeni benedeni. It must be noted that the growth increment analysis here is very simple. The primary aim of analysing its growth history was inferring whether this subspecies grew at a sufficient rate to allow for sub-annual resolution sampling. Oxygen isotope analysis already confirmed this was the case. Therefore, more sophisticated techniques, such as using cathodoluminescence to identify increments or applying spectral analysis to determine their periodicity, were beyond the scope of this study. However, a more detailed sclerochronology study with more specimens would be a good avenue for further work.

5.4 Angulus benedeni benedeni as a climate archive

Angulus benedeni benedeni shows promise as an archive for climate reconstructions. It can live for up to a decade or longer and enables the reconstruction of multiannual climate records at a seasonal resolution. Both δ¹⁸O_c and Δ₄₇ analyses can be successfully applied to this subspecies, and with a larger dataset it would be possible to obtain clumped-isotope-based summer and winter temperatures. We estimate that 3–5 specimens would be required to obtain confidence intervals on Δ₄₇-based seasonal temperatures that are ≤ 3.0 ^∘C. As confidence intervals scale with 1 $/$ N², progressively more aliquots would be needed to bring down the error, and so reducing errors beyond this threshold would be challenging. The relatively thin shell in combination with the long lifespan makes it somewhat challenging to sample at the required high temporal resolution.

The subspecies is especially suited for reconstructions of climate in the Pliocene of the North Sea basin, where it is common (De Meuter and Laga, 1976). However, the genus Angulus dates back to the late Eocene (Marquet, 2004, 2005) and may be used to reconstruct older climates as well, given that these older specimens are well-preserved. Potential target intervals include the Eocene–Oligocene transition, the Oligocene–Miocene transition, the Middle Miocene Climatic Optimum, and the mid-Piacenzian Warm Period (Westerhold et al., 2020). This first palaeoclimatological study on Angulus highlights the opportunities for using this genus as a climate archive and confirms earlier findings that tellinid bivalves are good subjects for palaeoclimatic research (Jones et al., 2005; Warner et al., 2022). Avenues of research are not limited to stable isotope analysis but may extend to minor and trace elemental analyses by, for example, extended μXRF analyses and laser-ablation ICP-MS profiles (de Winter et al., 2017b, 2020).

5.5 Clumped- and oxygen-isotope-based palaeotemperatures

The reconstructed mean annual temperature (MAT) of 13.5 ± 3.8 ^∘C based on clumped isotope thermometry is 3.5 ^∘C higher than that of the pre-industrial North Sea (Mackenzie and Schiedek, 2007; Emeis et al., 2015). Our temperature reconstructions for the Piacenzian in the North Sea are similar to a previous estimate of ca. 13 ^∘C for the North Sea from the somewhat older shallow marine Coralline Crag Formation in southeastern England (Dowsett et al., 2012). Although the complex stratigraphic context cannot exclude an age younger than the mid-Piacenzian Warm Period (mPWP), modelled anomalies for the mPWP North Sea area of around +3.5 ^∘C are in close agreement with our finding (Haywood et al., 2020). Compared to global estimates for the mPWP, a 3.5 ^∘C warming is on the higher end of the range given by proxy data (+2–4 ^∘C, sea surface temperature, Dowsett et al., 2012) and model data (+1.7–5.2 ^∘C, surface air temperature, Haywood et al., 2020). Overall, the reconstructed MAT agrees with both proxy and model data for this stratigraphic interval, supporting the validity of the clumped isotope data obtained from the shell of A. benedeni benedeni.

Seasonal temperature ranges based on δ¹⁸O_c corrected with mean annual δ¹⁸O_sw from clumped isotope analysis show a range of 12.1 ^∘C, from 6.4 ± 3.9 ^∘C in winter to 18.5 ± 3.9 ^∘C in summer (Fig. 10b). This is higher than the modern seasonal range of ca. 7–9 ^∘C (7–8 ^∘C in winter to 15–16 ^∘C in summer; Kooij et al., 2016). The average δ¹⁸O_sw of 0.10 ± 0.88 ‰ VSMOW is within the range of ca. −0.4 ‰ to +0.4 ‰ for the modern North Sea, although it is on the heavier side for coastal areas which are closer to −0.4 ‰ to −0.2 ‰ (Harwood et al., 2008). Regardless of the value of temperature reconstructions from this single locality for representing the mean North Sea condition, the calculated range is likely an underestimation. This is the result of inherent limitations of both the methodology and the bivalve-archive itself. Limits to the sampling resolution and growth breaks or stops during winter can result in some of the highest and lowest temperatures not being analysed. Furthermore, summer stratification of the water column would have dampened summer temperatures at the seafloor where A. benedeni benedeni lived. The upper Oorderen stratigraphic interval from which the specimens were collected is interpreted to have been shallow and strongly influenced by tidal currents, which is inconsistent with summer stratification. Still, as the water depth is not well constrained and the thermocline depth is unknown, summer stratification cannot be ruled out. Finally, the calculated seasonal temperatures assume a constant δ¹⁸O_sw. A δ¹⁸O_sw that varies throughout the year can either dampen or amplify the δ¹⁸O_c signal. A dampening corresponds to a positive correlation between δ¹⁸O_sw and temperature, e.g. due to enhanced evaporation during summer and enhanced precipitation during winter–spring. During summer, high temperatures will drive δ¹⁸O_c to lower values. This is dampened by a high ambient δ¹⁸O_sw due to evaporation, resulting in overall higher δ¹⁸O_c values in summer. The reverse is true in the winter. High summer evaporation and high winter–spring precipitation is expected in mid-latitude settings such as the North Sea, and a strong positive correlation between δ¹⁸O_sw and temperature is observed in the eastern North Sea today (Ullmann et al., 2010). These four factors – limited sampling resolution, possible growth slow-down in winter, possible stratification in summer, and possible dampening due to δ¹⁸O_sw fluctuations – all contribute to a potential underestimation of the actual seasonality. With the exception of seasonal δ¹⁸O_sw variability, these limitations also apply to clumped isotope seasonality reconstructions. Clumped isotope palaeothermometry has the additional caveat that it will almost always be necessary to combine many summer and winter data points that may not all represent “peak” summer and winter. This too can lead to some degree of seasonal dampening. Any further study on bivalve-derived seasonality in the region should address these concerns. In this context, the seasonal range recorded by these A. benedeni benedeni specimens may have been larger than that of the North Sea today. This enhanced seasonality is in line with several bivalve-based studies of the Pliocene southern North Sea (Johnson et al., 2009; Valentine et al., 2011; Johnson et al., 2022) and supports a more extensive study of Pliocene seasonality on Angulus bivalves.

6 Conclusions

Bivalves are valuable climate archives as they have the potential to record environmental changes on a seasonal scale. Here, we show that the bivalve Angulus benedeni benedeni shows promise as a climate archive. Based on microscopy and oxygen isotope measurements, this subspecies lived for multiple years and grew fast enough to allow for sampling its shell at a seasonal resolution. It likely experienced slower growth in winter and may not have recorded the coldest months. The reliability of the stable isotope results is confirmed by a diagenetic assessment of the shells using XRD, XRF, and EBSD, all of which attest to the pristine nature of the shell material. Furthermore, it reveals that the shell of A. benedeni benedeni consists of the following microstructures: compound composite prismatic, crossed lamellar, complex crossed lamellar, prismatic, and homogeneous. This is comparable to what has been observed in other tellinid bivalves. The clumped isotope analysis revealed that the mean annual temperature in the mid-Piacenzian southern North Sea basin was 13.5 ± 3.8 ^∘C, 3.5 ^∘C warmer than the pre-industrial average for this region. This is in line with model and proxy estimates for the stratigraphic interval from which the specimens originated. Overall, the findings of this study open up new avenues of research into the seasonality of past climates using A. benedeni benedeni, as well as other Angulus bivalves.

Appendix A: X-ray analyses

A1 X-ray diffraction

A part of specimen SG-125 was analysed on a Bruker D8 Advance diffractometer at Utrecht University and was ground to < 10 µm in a Retsch McCrone mill with zirconium oxide grinding elements prior to front-loading it into a PMMA sample holder (diameter: 25 mm; depth: 1 mm). The instrument settings are given in Table A1. The instrument was calibrated with a corundum crystal.

A2 Micro-X-ray fluorescence

Major and trace elements were measured on specimens SG-125, SG-126, and SG-127 using μXRF line scanning. Analysis was carried out on a benchtop Bruker M4 Tornado μXRF scanner at the Analytical, Environmental, and Geochemistry Research Group (AMGC) of the Vrije Universiteit Brussel (Brussels, Belgium). The Bruker M4 Tornado is equipped with a 30 W rhodium anode metal-ceramic X-ray tube operated at 50 kV and 600 µA. A polycapillary lens focuses an X-ray beam on the polished sample surface on a spot with a diameter of 25 µm (calibrated for Mo kα radiation). Fluorescing X-ray radiation is detected using two 30 mm² silicon drift detectors with a spectral resolution of 135 eV (calibrated for Mn-kα radiation). The system is operated under near-vacuum conditions (20 mbar) to allow lower-energy radiation of lighter elements (e.g. Mg) to be measured. More details on the μXRF measurement set-up and acquisition are provided in de Winter and Claeys (2017) and Kaskes et al. (2021).

X-ray spectra were produced using a spatial resolution of 25 µm and a dwell time of 60 s. This setting allows for enough time to reach the time of stable reproducibility and provides the optimal balance between longer measurement time (resulting in better defined XRF spectra) and higher sample sizes (see discussion in de Winter et al., 2017a). Spectra were quantified using online calibration with the matrix-matched standard BAS-CRM393 (Bureau of Analysed Samples, Middlesbrough, UK) in the M4 software (following de Winter et al., 2017a). After quantification, the trace element results were calibrated using an offline calibration constructed using a set of matrix-matched carbonate standards (CCH-1; Université de Liège, Belgium, COQ-1; US Geological Survey, Denver, CO, USA, CRM393, CRM512, CRM513, ECRM782; Bureau of Analysed Samples, Middlesbrough, UK; and SRM-1d; National Institute of Standards and Technology, Gaithersburg, MD, USA; see e.g. de Winter et al., 2021b).

Table A1 Instrument settings for X-ray diffraction analyses.

Appendix B: Stable isotopes

Table B1 Position of each sample taken from specimens SG-126 and SG-127, measured from the umbo to the centre line of the shell (along the axis of maximum growth).

Table B2 Accepted values for stable isotope standards, both oxygen and clumped.

Figure B1 Sampling procedure of the shells for stable isotope analysis. (a) Sample tracks on specimen SG-126, indicated by dotted lines; 55 in total, each fifth sample track has been labelled. The figure shows the sampled half (top) and that same half prior to sampling (bottom). (b) Sample tracks on specimen SG-127, indicated by dotted lines, 30 in total. This specimen was for the large part covered by epoxy, so only a small section at the upper surface could be sampled.

Figure B2 Illustration of the importance of not drilling too deep. In scenario 1, the two outermost layers (M+2, M+1) are drilled, whose material was precipitated at the same time. In scenario 2, slightly younger material from M+1 is drilled as well. In scenario 3, material from below the pallial myostracum which separates M+1 and M−1 is drilled. As a result of the way bivalves grow, this results in mixing materials of very different precipitation ages. In addition, in this scenario materials are precipitated from different internal fluids (dorsal and ventral of the pallial line). In this study, scenario 1 was the goal. Some sampling according to scenario 2 may have been present, but scenario 3 was avoided.

Figure B3 Linear correlation between the Naxos δ¹⁸O data (11 standards measured) and the mass 44 intensity that was used to correct the δ¹⁸O data. The correction was required due to the recommended weight range for each aliquot as per the laboratory’s protocol it was quite large (50–100 µg). This results in a broad range of mass 44 values as well.

Figure B4 Empirical transfer function (ETF) constructed from a linear regression between accepted and measured Δ₄₇ values from standards ETH-1, ETH-2, and ETH-3. This function was used to correct all samples and standards and transfer them to the I-CDES90C reference frame. All standards plotted here have already been corrected for potential drift during the measurement by the ETH-3 bracketing correction. The standards were measured during 28 runs over a period of 2 months.

Figure B5 Corrected Δ₄₇ data for external (ETH-3, ETH-2, ETH-1) and internal (IAEA-C2, Merck) standards. (a)–(c): ETH-3 Δ₄₇, 720 δ¹⁸O, and δ¹³O data. Each point represents the measurement of a single aliquot. (d)–(f): ETH-2 Δ₄₇, δ¹⁸O, and δ¹³O data. (g)–(i): ETH-1 Δ₄₇, δ¹⁸O, and δ¹³O data. (j)–(l): IAEA-C2 Δ₄₇, δ¹⁸O, and δ¹³O data. (m)–(o): Merck Δ₄₇, δ¹⁸O, and δ¹³O data. δ¹⁸O and δ¹³O measurements affected by mass-dependent fractionation are marked as red open circles (see Sect. 3.5.3).

Figure B6 Raw δ¹⁸O_c data for specimens SG-126 (top) and SG-127 (bottom), measured on the two different instruments as detailed in Sect. 3.5.2 and 3.5.3. The raw GasBench data shown here are not corrected for mass 44 variability yet (see Fig. B3), which explains the offset between the two records. All analyses described in the main text and subsequent plots concern the corrected GasBench data.

Figure B7 Alternative growth curves for specimen SG-126, counting the peak at 21 mm as a year (marked with red arrow in panel b). (a) Growth curves constructed by inserting the counted years versus distance-from-umbo data points into two different fitting algorithms: the asymptotic Von Bertalanffy growth function and the logistic Gompertz function. Compared to the growth curves that do not include this peak at 21 mm, the residuals are slightly worse for both models (0.51 vs. 0.58 for the Von Bertalanffy model, 0.63 vs. 0.84 for the Gompertz model). The estimated maximum shell height is similar for the Gompertz model (34.11 vs. 31.34 mm), but significantly shorter for the Von Bertalanffy model (52.35 vs. 40.34 mm). (b) δ¹⁸O_c record with the inferred growth years marked by lines connecting to plot (a).

Appendix C: Electron backscatter diffraction

Figure C1 EBSD map of specimen SG-125 indicating the absence secondary calcite crystals. Map acquisition specifics: 10 kV voltage, 1 µm step size. Measured at low vacuum due to excessive charging. Classification: 54.6 % aragonite, 45.2 % unclassified, 0.2 % calcite. The map consists of the band contrast overlain by the classified pixels, which are either aragonite (blue) or calcite (red). M+2, M+1, and M−1 correspond to the shell layers as described in Sect. 4.1.1. No clean-up was done on this map. The pixels identified as calcite do not show any of the characteristics of actual secondary calcite (large crystals compared to the fine-grained aragonite, breaking up the aragonitic structure), are always single pixels, and are removed through AZtecCrystal's wild spike removal algorithm. Several dust shadows are also seen in this image; these areas are unclassified. There are also several scratches present on the surface.