Multi-proxy assessment of brachiopod shell calcite as a potential archive of seawater temperature and oxygen isotope composition

Most of our knowledge of past seawater temperature history is based on δ¹⁸O values of calcium carbonate fossil shells. However, the determination of past temperatures using this proxy requires the knowledge of past seawater δ¹⁸O values, which is generally poorly constrained. Other paleothermometers using carbonate archives, such as $\mathrm{Mg}/\mathrm{Ca}$ ratios and clumped isotopes (Δ₄₇), have been developed to allow for paleotemperatures to be estimated independently and to allow past ocean δ¹⁸O values to be calculated using various groups of calcifying organisms. Articulated brachiopod shells are some of the most commonly used archives in studies of past oceanic geochemistry and temperature. They are abundant in the fossil record since the Cambrian, and for decades, their low Mg–calcite mineralogy has been considered relatively resistant to diagenetic alteration. Here, we investigate the potential of brachiopod shells as recorders of seawater temperatures and seawater δ¹⁸O values using new brachiopod shell geochemical data by testing multiple well-established or suggested paleothermometers applied to carbonate archives.

Modern articulated brachiopod shells covering a wide range of temperatures (−1.9 to 25.5 ^∘C), depths (5 to 3431 m) and salinities (33.4 to 37.0 PSU) were analysed for their stable isotope compositions (δ¹³C, δ¹⁸O and Δ₄₇) and their elemental ratios ($\mathrm{Mg}/\mathrm{Ca}$, $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Na}/\mathrm{Ca}$ and $\mathrm{Li}/\mathrm{Ca}$). Our data allowed us to propose a revised oxygen isotope fractionation equation between modern-brachiopod shell calcite and seawater:

$$\begin{array}{}\text{(1)}& T=-\mathrm{5.0}(\pm \mathrm{0.2})({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}-{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sw}})+\mathrm{19.4}(\pm \mathrm{0.4}),\end{array}$$

where δ¹⁸O_c is in ‰ VPDB, δ¹⁸O_sw is in ‰ VSMOW, and T is in ^∘C. Our results strongly support the use of clumped isotopes as an alternative temperature proxy but confirm significant offsets relative to the canonical relationship established for other biogenic and abiogenic calcium carbonate minerals. Brachiopod shell $\mathrm{Mg}/\mathrm{Ca}$ ratios show no relationship with seawater temperatures, indicating that this ratio is a poor recorder of past changes in temperatures, an observation at variance with several previous studies. Despite significant correlations with brachiopod living temperature, brachiopod shell $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Na}/\mathrm{Ca}$ and $\mathrm{Li}/\mathrm{Ca}$ values indicate the influence of environmental and biological factors unrelated to temperature, which undermines their potential as alternative temperature proxies. Kinetic effects (growth rates) could explain most of the deviation of brachiopod shell calcite from expected isotopic equilibrium with seawater and part of the distribution of $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Na}/\mathrm{Ca}$ and $\mathrm{Li}/\mathrm{Ca}$ ratios.

1 Introduction

Oxygen isotope ratios (δ¹⁸O) of brachiopod shell calcite constitute the most extensive record of marine temperatures over the Phanerozoic eon (Prokoph et al., 2008; Veizer and Prokoph, 2015). Nevertheless, our limited knowledge of δ¹⁸O values of past seas and oceans prevents a confident interpretation of this record, as carbonate δ¹⁸O values (δ¹⁸O_c) are dependent on both shell growth temperature and living-water δ¹⁸O (δ¹⁸O_sw; Epstein et al., 1953; Kim and O'Neil, 1997; Kim et al., 2007; Brand et al., 2013, 2019). Other paleotemperature proxies, such as trace element ratios and clumped isotopes (Δ₄₇), may be the key to building a more confident record of Phanerozoic ocean temperatures. Owing to the determination over the last decades of reliable oxygen isotope fractionation equations between calcite and seawater, such independent temperature estimates can be combined with δ¹⁸O_c measurements to reconstruct δ¹⁸O_sw values and hence to characterize oceanic or local hydrography.

Alternative carbonate archive thermometers have been developed for different calcifying organisms and have been proposed in the literature as potentially applicable to brachiopod shells. Element $/$ Ca ratios such as $\mathrm{Mg}/\mathrm{Ca}$, $\mathrm{Sr}/\mathrm{Ca}$ and $\mathrm{Li}/\mathrm{Ca}$ are the most common in the literature. $\mathrm{Mg}/\mathrm{Ca}$ thermometry of foraminiferal calcite (Nürnberg, 1995; Anand et al., 2003) is a widely used proxy that is the basis for the reconstruction of changes in δ¹⁸O_sw and ice volumes over the Cenozoic era (Lear, 2000; Billups and Schrag, 2003; Miller et al., 2020). Building on the pioneering work of Lowenstam (1961), Brand et al. (2013, 2019) documented a positive correlation between seawater temperature and the MgCO₃ content of brachiopod shells sampled at a worldwide scale. Plus, following the conclusions of Jiménez-López et al. (2004), Brand et al. (2013) developed an isotopic fractionation equation that considers the amount of Mg²⁺ substituted for Ca²⁺ in the crystal lattice of the low-Mg calcitic shell which was abandoned in their more recent works (Brand et al., 2019). The MgCO₃–temperature relationship is described by a curve named the Global Brachiopod Mg Line (GBMgL), yet Brand et al. (2013, 2019) acknowledge species-specific deviations from this concept that fall well above the said line. Plus, brachiopod shell Mg content has been found to have a strong taxonomic trend, with shells of thecideid and craniid brachiopods being made of high-Mg calcite in contrast to the low-Mg calcite of rhynchonellide and terebratulide brachiopod shells (Brand et al., 2003; Ullmann et al., 2017). This undermines the brachiopod shell MgCO₃–temperature relationships which uses a dataset that consists of thecideid, rhynchonellide and terebratulide (Brand et al., 2013, 2019).

While $\mathrm{Sr}/\mathrm{Ca}$ ratios are mostly used to assess the shell preservation of fossil specimens (Brand and Veizer, 1980), this ratio also constitutes a potential marine paleothermometer applicable to brachiopod shells, as is the case for corals (McCulloch et al., 1994; Shen et al., 1996; Marshall and McCulloch, 2002; Swart et al., 2002; Ayling et al., 2006; DeLong et al., 2011). Previous studies already suggested $\mathrm{Sr}/\mathrm{Ca}$ ratios as marine paleothermometers for brachiopod shells owing to their correlations with oxygen isotope ratios and seawater temperature (Lowenstam, 1961; Mii and Grossman, 1994; Pérez-Huerta et al., 2008; Butler et al., 2015; Ullmann et al., 2017). Similarly, the brachiopod $\mathrm{Li}/\mathrm{Ca}$ ratio was also suggested as a potential temperature proxy in the seminal study of Delaney et al. (1989). Subsequent studies have also confirmed that $\mathrm{Li}/\mathrm{Ca}$ may constitute a potentially accurate thermometer in brachiopods (Dellinger et al., 2018; Rollion-Bard et al., 2019; Washington et al., 2020) and in corals (Montagna et al., 2014; Marchitto et al., 2018). For this latter group, however, the $\mathrm{Li}/\mathrm{Mg}$ ratio has been considered to be a much more precise thermometer than $\mathrm{Li}/\mathrm{Ca}$ (Montagna et al., 2014; Marchitto et al., 2018) despite some identified caveats such as the scrambling role of organic matter (Cuny-Guirriec et al., 2019).

On top of element $/$ Ca proxies, clumped isotope (Δ₄₇) thermometry has been developed more recently in relation to different carbonate materials (Ghosh et al., 2006; Zaarur et al., 2013). This parameter measures the anomaly of ¹³C–¹⁸O bonds within the carbonate lattice relative to their stochastic abundance. Δ₄₇ values are strongly correlated to crystallization temperature and are independent from δ¹⁸O_sw. This proxy has been applied to different marine calcifying organisms to estimate δ¹⁸O_sw in the deep past (Petersen et al., 2016; Bergmann et al., 2018; Henkes et al., 2018; Wierzbowski et al., 2018; Price et al., 2020; Vickers et al., 2020, 2021; de Winter et al., 2021; Meckler et al., 2022). Although promising, brachiopod shell Δ₄₇ composition was only used in a couple of studies to estimate past temperatures (Henkes et al., 2013; Came et al., 2014). Indeed, the few existing modern-brachiopod Δ₄₇–temperature data are in good agreement with the canonical Δ₄₇–temperature relationship established for other biogenic and abiogenic calcium carbonate minerals (Henkes et al., 2013; Came et al., 2014). Recent compilation of Δ₄₇–T calibration data from different laboratories using the new carbonate-based standardization (Bernasconi et al., 2018, 2021) concluded with a unique Δ₄₇–T relationship for synthetic and natural carbonates (Anderson et al., 2021). The establishment of such a unique calibration equation strengthens the good agreement observed between natural and synthetic Δ₄₇–T relationships and would suggest that this calibration is valid when applied to brachiopod shells, especially as empirical calibrations are in very good agreement with theoretical predictions (Jautzy et al., 2020). However, Δ₄₇ measurements from brachiopod shells are absent from this compilation, and other studies suggest that they may be subjected to significant deviations from Δ₄₇–T relationships that are likely related to differences in growth rates (Bajnai et al., 2018, 2020).

Beside, a number of studies suggested that living temperature alone may not explain the variability of those different parameters within brachiopod shells. First, biogenic carbonate precipitation rate depends on environmental and biological factors (temperature and salinity tolerance, age, reproduction period, ontogeny, etc.), which causes biases in the paleoenvironment record extracted from the shell (Peck et al., 1997; Schöne, 2008). Brachiopod shell calcite is secreted by the outer epithelial cells of the mantle (Williams, 1968; Simonet Roda et al., 2019a, b, 2022). Calcite elements are in close relationship with organic matrices which control the mineral growth (Williams, 1968; Gaspard et al., 2008; Gaspard and Nouet, 2016; Simonet Roda et al., 2019a, b, 2022). Finally, nanoparticles are transported via ion transport through the cell membrane, highlighting calcite precipitation as an amorphous calcium carbonate precursor stage (Griesshaber et al., 2009; Schmahl et al., 2012). The chemistry of the calcifying environment of brachiopod shell calcite is poorly constrained (Immenhauser et al., 2016) Testing thermodynamic equilibrium at the mineralization site is therefore, to date, out of reach. Consequently, the chemistry of brachiopod shells is usually compared to that of the surrounding water in order to approach the processes involved in biomineralization from a chemical point of view. Numerous studies highlighted that brachiopod shell calcite does not precipitate sensu stricto in isotopic equilibrium with the surrounding water but, at best, mimics thermodynamic equilibrium (Carpenter and Lohmann, 1995; Auclair et al., 2003; Parkinson et al., 2005; Yamamoto et al., 2010a, b; Cusack et al., 2012; Takayanagi et al., 2012, 2013, 2015; Bajnai et al., 2018, 2020; Romanin et al., 2018; Rollion-Bard et al., 2019). Brand et al. (2019) characterize their oxygen isotope fractionation equation as representing a brachiopod based-equilibrium. This term is opposed to the thermodynamic equilibrium, which is best described by very-slow-growing calcites (Coplen, 2007; Daëron et al., 2019), thus ensuring isotopic equilibration of the dissolved inorganic carbon (DIC) species with water (Watkins et al., 2013, 2014). Trace element incorporation into calcite is also affected by kinetic effects. In synthetic carbonates, the relative abundance of elements such as Sr, Na and Li rises with higher growth rates (Lorens, 1981; Busenberg and Plummer, 1985; Tesoriero and Pankow, 1996; Gabitov et al., 2011, 2014). Such kinetic trends were evidenced at an intra-individual level among brachiopods by the negative correlation of element $/$ Ca ratios with stable isotope values (δ¹³C, δ¹⁸O; Ullmann et al., 2017; Rollion-Bard et al., 2019).

In this study, we assess the potential of various geochemical proxies (δ¹⁸O, Δ₄₇, $\mathrm{Sr}/\mathrm{Ca}$; $\mathrm{Mg}/\mathrm{Ca}$, $\mathrm{Li}/\mathrm{Ca}$, $\mathrm{Na}/\mathrm{Ca}$) as recorders of seawater temperature by analysing a new set of modern articulated brachiopod shells collected from various depths and latitudes during institutional oceanographic cruises and covering a broad range of water temperatures comprised between −1.9 and 25.5 ^∘C. We discuss the validity and robustness of our oxygen isotope fractionation equation established with measured δ¹⁸O_c of calcite and with seawater temperature and δ¹⁸O_w estimated from oceanographic data. The dependence of elemental ratios on seawater temperature is discussed in light of this new available dataset of modern brachiopods. Finally, we highlight kinetic effects as non-negligible sources of isotopic and trace element variability.

2 Material and methods

2.1 Sample collection and environmental parameters

The studied material consists of 37 articulated brachiopod shells that were collected in situ during institutional oceanographic cruises or other scientific missions (see Supplement for details). The sampling location was documented with depth and geographic coordinates for most samples. Temperature and salinity, which were not measured in situ for most samples, were estimated independently using the NOAA World Ocean Atlas 2018 (Locarnini et al., 2018; Zweng et al., 2018) or, whenever available, using long-term local records in the literature (Table 1). For both temperature and salinity, mean annual values (MAT), as well as higher and lower monthly averages, were considered. Similarly, δ¹⁸O_sw values were usually not measured in situ but were calculated from the salinity data with the appropriate regional δ¹⁸O_sw–S relationship from LeGrande and Schmidt (2006). The oceanographic parameters used in this study are listed in Table 1.

Table 1 Brachiopod samples with taxonomic identification, sampling location and environmental parameters. Most environmental parameters are derived from the World Ocean Atlas 2018 (Locarnini et al., 2018; Zweng et al., 2018). For a few samples, environmental parameters are derived from local literature, as listed below. Almost all δ¹⁸O_sw values are calculated using the regional δ¹⁸O_sw–salinity relationship published by LeGrande and Schmidt (2006), unless direct measurements covering seasonal variations were available in the literature.

¹ Goodwin and Cornelisen (2012); ² Woods et al. (2014); ³ Meredith et al. (2013); ⁴ Jacobson (1983).

2.2 Sample preparation

Encrusting organisms covering the shells were first mechanically removed using stainless steel dental tools, and each shell was then placed in diluted bleach (NaClO 5 %) for 5 to 10 min in an ultrasonic bath to remove organic matter and other contaminants. As preliminary tests have revealed that the more delicate shells (e.g. Macandrevia africana; Cooper, 1975) broke apart during the ultrasonic bath, immersion in diluted bleach for a few hours without the ultrasonic bath was subsequently preferred for a few specimens. In all cases, the specimens were rinsed with deionized water and oven dried at 50 ^∘C for a few hours.

The sampling of brachiopod shell calcite was performed following the recommendations of Romanin et al. (2018). The umbo, edges and muscle scar area were avoided because they record major kinetic or metabolic effects (Carpenter and Lohmann, 1995; Auclair et al., 2003; Parkinson et al., 2005; Yamamoto et al., 2010b; Ullmann et al., 2017; Romanin et al., 2018), and sampling was focused on the middle part of the shell. Considering the complex structure of brachiopod shells (Williams, 1968; Gaspard et al., 2018; Simonet Roda et al., 2022) relative to the precision of our sampling method (dental tools and engraving bit), we were unable to prepare pure samples for each layer of the brachiopod shell (acicular primary layer, fibrous secondary layer and columnar tertiary layer). In this study, we differentiate the outer layers, i.e. a mix of the primary acicular layer and some amount of the outer secondary fibrous layer, which record significant kinetic effects (Carpenter and Lohmann, 1995; Auclair et al., 2003; Parkinson et al., 2005; Bajnai et al., 2018; Romanin et al., 2018; Rollion-Bard et al., 2019), from the inner layers, i.e. the inner secondary fibrous layer and/or the tertiary columnar layer when present. Outer shell layers were removed manually with dental tools, and resulting powders were kept to compare their elemental ratios, δ¹⁸O and δ¹³C with those of the inner shell layers (Carpenter and Lohmann, 1995; Auclair et al., 2003; Pérez-Huerta et al., 2008). After removal of the outer layers, the inner layers were sampled using an engraving bit fitted to a DREMEL Micro™ drill adjusted to the lowest possible speed. For each sample, the powder was used for all selected geochemical analyses which correspond to δ¹⁸O, δ¹³C, Δ₄₇ and element concentration. As clumped isotope analyses require large amounts of material (>10 mg is necessary to operate multiple replicate analyses), a large area (a few cm²) was sampled in the middle part of the shell, allowing several mg of calcitic powder to be collected. It is worth noting that this large sampling area may correspond to several months or years of shell growth, depending on the species.

This sampling protocol was applied to most of our samples, except for the smallest and most fragile shells that quickly broke apart during sampling. Consequently, in the case of the small shells of Terebratulina latifrons, Frenulina sanguinolenta and some specimens of Fallax neocaledonensis, the outer layers were kept, and only the area covering the umbo up to the muscle scar was avoided for sampling. The weakness and thinness of M. africana precluded a total removal of the outer layers; however, the umbo, muscle scars and edges were removed from the bulk shell. Fragments of these more fragile shells were ground to a fine powder in an agate mortar.

2.3 Carbon and oxygen stable isotopes

Stable isotope compositions of the sampled powders were determined using a MultiPrep™ Autosampler coupled to a dual-inlet GV IsoPrime^® mass spectrometer. For each sample, an aliquot of about 400 µg of calcium carbonate was reacted with anhydrous oversaturated phosphoric acid at 90 ^∘C for 20 min. Oxygen isotope ratios of calcium carbonate were computed assuming an acid fractionation factor 1000ln α(CO₂–CaCO₃) of 8.1 between carbon dioxide and calcite (Swart et al., 1991). All sample measurements were duplicated and adjusted to the international references NIST NBS18 (δ¹⁸O_VPDB = −23.2 ‰; δ¹³C_VPDB = −5.01 ‰) and NBS19 (δ¹⁸O_VPDB = −2.20 ‰; δ¹³C_VPDB = +1.95 ‰) and an internal standard of Carrara Marble (δ¹⁸O_VPDB = −1.84 ‰; δ¹³C_VPDB = +2.03 ‰). Since 2019, reproducibility of the Carrara Marble in-house standard has been ±0.096 ‰ for δ¹⁸O (2σ, n=1062) and ±0.066 ‰ for δ¹³C (2σ, n=1062). These stable isotope data are completed with δ¹⁸O and δ¹³C data obtained from clumped isotope measurements, for which the analytical procedures are reported in detail in Sect. 2.5.

2.4 Elemental ratios

Elemental concentrations were obtained by dissolving 2 to 20 mg of carbonate powder in 10 mL HNO₃ (2 %). For all samples, pairs of aliquots were prepared with 10× and 100× dilutions, depending on the considered trace and major elements, along with a fixed amount of Sc and In that were added to correct concentrations from instrument drift. Solutions were analysed using an inductively coupled plasma-optical emission spectrometer (iCAP 7000 ICP-OES) and a quadrupole ICP-mass spectrometer (i-CAP-Q ICP-MS) for minor and trace elements, respectively. The 100× diluted aliquots were used to calculate Ca concentration, while the 10× diluted aliquots served for Mg, Sr, Na and Li concentrations. The calculated element $/$ Ca ratios are reported in mmol $/$ mol for $\mathrm{Mg}/\mathrm{Ca}$, $\mathrm{Sr}/\mathrm{Ca}$ and $\mathrm{Na}/\mathrm{Ca}$ and in µmol $/$ mol for $\mathrm{Li}/\mathrm{Ca}$. The reproducibility of measurements was assessed through the analysis of the carbonate standard CCH1 (Roelandts and Duchesne, 1988).

2.5 Clumped isotopes

Carbonate samples were converted to CO₂ by phosphoric acid reaction at 90 ^∘C in a common, stirred acid bath for 15 min. The initial phosphoric acid concentration was 103 % (1.91 g cm⁻³), and each batch of acid was used for 7 d. After cryogenic removal of water, the evolved CO₂ was helium-flushed at 25 mL min⁻¹ through a purification column packed with Porapak Q ($\mathrm{50}/\mathrm{80}$ mesh, 1 m length, 2.1 mm ID) and held at −20 ^∘C, then it was quantitatively recollected by cryogenic trapping and transferred into an IsoPrime 100™ dual-inlet mass spectrometer equipped with six Faraday collectors ($m/z$ 44–49). Each analysis took about 2.5 h, during which analyte gas and working reference gas were allowed to flow from matching 10 mL reservoirs into the source through deactivated fused-silica capillaries (65 cm length, 110 µm ID). Every 20 min, gas pressures were adjusted to achieve $m/z$ = 44 current of 80 nA, with differences between analyte gas and working gas generally being below 0.1 nA. Pressure-dependent background current corrections were measured 12 times for each analysis. All background measurements from a given session were then used to determine a mass-specific relationship linking background intensity (Z_m), total $m/z$ = 44 intensity (I₄₄) and time (t) as follows:

$$\begin{array}{}\text{(2)}& Z_m=a{I}_{\mathrm{44}}+P\left(t\right),\end{array}$$

with P being a polynomial of degree 2 to 4.

Background-corrected ion current ratios (δ₄₅ to δ₄₉) were converted to δ¹³C, δ¹⁸O and “raw” Δ₄₇ values as described by Daëron et al. (2016) using the IUPAC oxygen-17 correction parameters. The isotopic composition (δ¹³C, δ¹⁸O) of our working reference gas was computed based on the nominal isotopic composition of carbonate standard ETH-3 (Bernasconi et al., 2018) and an oxygen-18 acid fractionation factor of 1.00813 (Kim et al., 2007). Raw Δ₄₇ values were then converted to the I-CDES (Intercarb-Carbon Dioxide Equilibrium Scale) Δ₄₇ reference frame by comparison with four ETH carbonate standards (ETH 1–4; Bernasconi et al., 2021) using a pooled-regression approach (Daëron, 2021). Full analytical errors are derived from the external reproducibility of unknowns and standards (N_f = 89) and conservatively account for the uncertainties in raw Δ₄₇ measurements, as well as those associated with the conversion to the absolute Δ₄₇ reference frame.

3 Results

All geochemical results obtained from modern brachiopods (δ¹³C, δ¹⁸O, Δ₄₇, $\mathrm{Mg}/\mathrm{Ca}$, $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Na}/\mathrm{Ca}$ and $\mathrm{Li}/\mathrm{Ca}$) are reported in Table 2. Linear regression models of Δ¹⁸O_c-w, Δ₄₇, $\mathrm{Mg}/\mathrm{Ca}$, $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Na}/\mathrm{Ca}$ and $\mathrm{Li}/\mathrm{Ca}$, alongside mean annual temperature (MAT), are reported in Fig. 1 and Table 3.

Table 2 Geochemical data of modern-brachiopod samples. Stable isotope values and element $/$ Ca ratios.

^* Samples with the anterior part of the shell.

Table 3 Linear regression model parameters for the geochemical parameter–temperature relationship tested and displayed in Fig. 1.

Figure 1 Geochemical data and linear regression model for each geochemical parameter tested for its relationship with brachiopod living temperature. Parameters associated with each regression model are listed in Table 3.

3.1 Carbon and oxygen stable isotopes

The whole dataset of δ¹³C and δ¹⁸O values (n=73) from modern-brachiopod shells ranges from −2.2 ‰ to 3.2 ‰ and from −1.2 ‰ to 4.0 ‰, respectively. There is a robust negative linear correlation (R²>0.73, p slope < 0.001) between Δ¹⁸O_c-w (=δ¹⁸O_c − δ¹⁸O_w) and MAT (Fig. 1; Table 3). The correlation observed in bulk shell samples is indistinguishable from that in inner-layer samples (Fig. 1a). On the contrary, the correlation obtained from outer-layer samples shows a significant offset in relation to both inner-layer and bulk-layer samples. This is in line with the differences between the outer and inner layers observed within the same specimen, which are in most cases above analytical uncertainties and range from −3.7 ‰ to −0.2 ‰ (mean = −1.4 ± 0.4 ‰; 2σ; n=17) and −0.8 ‰ to 0.5 ‰ (mean = $-\mathrm{0.4}\pm \mathrm{0.3}$ ‰; 2σ; n=17) for δ¹³C and δ¹⁸O, respectively. Our new set of data from inner- and bulk-layers samples provides the following oxygen isotope fractionation equation:

$$\begin{array}{}\text{(3)}& T=-\mathrm{5.0}(\pm \mathrm{0.2})\left({\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{c}}{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{sw}}\right)+\mathrm{19.4}(\pm \mathrm{0.4}),\end{array}$$

with T being the temperature in ^∘C and Δ¹⁸O_c-w (= δ¹⁸O_c−δ¹⁸O_w) being the oxygen isotope fractionation between brachiopod calcite and seawater, excluding data from the outer-shell layers. As opposed to the ordinary least-squares regression used in Fig. 1 and Table 3, this equation is established using a York regression (York et al., 2004) using seasonal variation for temperature uncertainties and propagated uncertainties from δ¹⁸O_c measurements and δ¹⁸O_w estimates for Δ¹⁸O_c-w uncertainties.

3.2 Trace element $/$ Ca ratios

The $\mathrm{Mg}/\mathrm{Ca}$, $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Li}/\mathrm{Ca}$ and $\mathrm{Na}/\mathrm{Ca}$ ratios of modern-brachiopod shells range from 3.57 to 25.30 mmol $/$ mol, 0.45 to 2.05 mmol $/$ mol, 3.28 to 58.47 µmol $/$ mol and 1.86 to 18.28 mmol $/$ mol, respectively (N=47). $\mathrm{Sr}/\mathrm{Ca}$ ratios show significant but weak negative correlations ($\mathrm{0.28}<R^{\mathrm{2}}<\mathrm{0.44}$; p slope < 0.01) with temperature. $\mathrm{Li}/\mathrm{Ca}$ ratios show significant and good negative correlations ($\mathrm{0.61}<R^{\mathrm{2}}<\mathrm{0.86}$; p slope < 0.001) with temperature, with the strongest correlation (R² = 0.86) being reached in outer-layer samples. $\mathrm{Na}/\mathrm{Ca}$ ratios show significant negative correlations ($\mathrm{0.38}<R^{\mathrm{2}}<\mathrm{0.60}$; p slope < 0.001) with temperature. The concentrations of these three elements (Sr, Li and Na) in brachiopod shells are strongly positively correlated with each other ($\mathrm{0.66}<R^{\mathrm{2}}<\mathrm{0.93}$ for the whole dataset), especially when considering only inner-layers sample ($\mathrm{0.84}<R^{\mathrm{2}}<\mathrm{0.96}$). No significant correlation (p slope > 0.2) of $\mathrm{Mg}/\mathrm{Ca}$ with temperature is observed in the dataset (Fig. 1; Table 3). Measured element $/$ Ca ratios are systematically higher in the outer layers than in the inner layers of brachiopod shells (Student's t test: p < 0.01; Fig. 1c–e).

3.3 Clumped isotopes

Δ₄₇ values obtained from the selected brachiopod shells range from 0.59±0.01 ‰ to 0.69±0.01 ‰ I-CDES and are strongly correlated with ambient temperature (Table 3; Fig. 1b). Temperatures inferred from clumped isotope data were calculated using the INTERCARB calibration determined by Anderson et al. (2021), who reprocessed the data coming from various laboratories using the same carbonate standards for data correction. Clumped isotope temperatures are significantly lower than estimated MAT, with a mean deviation from environmental temperatures of $-\mathrm{2.4}\pm \mathrm{1.3}$ ^∘C (95 % CI (confidence interval); n=19).

4 Interpretation and discussion

4.1 Brachiopod shell paleothermometers

4.1.1 Validity and robustness of our ${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}\mathrm{O}$-based fractionation equation

The oxygen isotope fractionation equation derived from our data deviates substantially from that determined by Brand et al. (2019), which is also derived from articulated brachiopod shells (Fig. 2). The two equations follow similar trends and are only offset by ∼ 0.6 ‰ at tropical temperatures (20–30 ^∘C; Fig. 2a), with the equation from Brand et al. (2019) predicting a lower fractionation factor. At temperate and polar temperatures (20 to 0 ^∘C), our new equation highlights a higher sensitivity of oxygen fractionation to temperature than that of Brand et al. (2019; Fig. 2). The two equations overlap at around 10 ^∘C and deviate substantially around ∼ 0 ^∘C, with a Δ¹⁸O_c-w difference of up to ∼ 1 ‰ between the equation of Brand et al. (2019) and our equation (Fig. 2a), a pattern also observed in the geochemical datasets (Fig. 2b). Examining the geochemical dataset revealed that the large offset between the two equations at low temperatures partly results from the use of a second-order regression model by Brand et al. (2019), as the regression curve overlaps with the highest Δ¹⁸O_c-w data points at temperatures of ∼ 1 ^∘C. We preferred to use a linear regression that was more appropriate to data in the low-temperature (<50 ^∘C) linear domain of the hyperbolic curve ($\mathit{\alpha }=f(\mathrm{1}/T^{\mathrm{2}})$) predicted by equilibrium thermodynamics. Still, the main discrepancies between the equation proposed herein and that of Brand et al. (2019) are mainly driven by a difference in the geochemical dataset. The case of the polar brachiopods Liothyrella uva from the island of Rothera (Antarctica) highlights the likely reasons for this discrepancy at low temperatures. Indeed, δ¹⁸O values reported in both studies are the same (Brand et al., 2019), indicating that the discrepancy lies in the environmental parameters used to derive the equation. Brand et al. (2019) used published environmental parameters measured locally in late February for this site, while we approached this estimation using regional annual averages. This difference in integration time of the environmental parameters is important because the month of February records the annual maximum and minimum in seawater temperature and δ¹⁸O_sw values, respectively (Meredith et al., 2013). This particular brachiopod species has a higher growth rate during the Austral winter than during the Austral summer (Peck et al., 1997), further increasing the likely seasonal bias. These discrepancies are individually relatively small (± 1–2 ^∘C; ±0.5 ‰ VSMOW), but put together, they easily explain the 1 ‰ offset between both equations in the low-temperature range.

Figure 2 Plot comparison of (a) regression lines for different temperature–oxygen isotope fractionations for brachiopod calcite (this study; Brand et al., 2019), laboratory-precipitated calcite (Kim and O'Neil, 1997) and slow-growing cave calcite (Daëron et al., 2019); and (b) corresponding datasets for brachiopod calcite (modified after Brand et al., 2019). Uncertainty envelope of 95 % CI of the York least-squares regression (York et al., 2004).

We use a published independent dataset of modern-brachiopod δ¹⁸O values (Bajnai et al., 2018) to test our equation against other published equations in their ability to predict environmental temperatures when δ¹⁸O_sw is known (Fig. 3). This dataset has a MAT range of 0 to 29 ^∘C and uses two sets of δ¹⁸O_sw for most data, one derived from direct measurements and one derived from δ¹⁸O_sw–salinity relationships (LeGrande and Schmidt, 2006). Plus, $\mathrm{Mg}/\mathrm{Ca}$ values are available for all samples, allowing us to apply the equation of Brand et al. (2013) to the dataset. Isotopic temperatures were calculated and then normalized as the temperature deviation from environmental temperature (ΔT in ^∘C) with full propagation of uncertainties. Of all the equations tested, the equations derived from measurements of modern brachiopods (Brand et al., 2013, 2019; this study) yield the most accurate temperature predictions. Mean ΔT are statistically indistinguishable from 0, using salinity-based δ¹⁸O_sw estimate results for all four equations. Interestingly, the accuracy of isotopic temperature estimates is poorer when measured δ¹⁸O_sw values are used. The equations of Brand et al. (2013, 2019) now significantly underestimate growing temperatures (Fig. 3; $p_{(\mathrm{mean}\mathrm{\Delta }T=\mathrm{0})}<\mathrm{0.01}$), while our equation yields mean ΔT values that are statistically indistinguishable from 0 ($p_{(\mathrm{mean}\mathrm{\Delta }T=\mathrm{0})}=\mathrm{0.11}$). The equation of Kim and O'Neil (1997), derived from synthetic carbonates, and that of Epstein et al. (1953), derived from molluscs, both underestimate the growing temperature of brachiopod shells (Fig. 3; $p_{(\mathrm{mean}\mathrm{\Delta }T=\mathrm{0})}<\mathrm{0.001}$) in both δ¹⁸O_sw hypotheses. On the other end of the spectrum, the equation of Daëron et al. (2019), derived from slow-growing calcite, overestimates the growing temperature of brachiopod shells (Fig. 3, $p_{(\mathrm{mean}\mathrm{\Delta }T=\mathrm{0})}<\mathrm{0.001}$) in both δ¹⁸O_sw hypotheses. Although it lessens the accuracy of δ¹⁸O temperatures in this dataset, the use of locally measured δ¹⁸O_sw values over the salinity-based estimates improved their precision ($\mathrm{4.0}<\mathrm{2}\mathit{\sigma }<\mathrm{5.3}$ ^∘C, n=13 rather than $\mathrm{5.4}<\mathrm{2}\mathit{\sigma }<\mathrm{6.4}$ ^∘C, n=18).

Figure 3 Offsets between the δ¹⁸O temperatures and MAT from the modern-brachiopod dataset of Bajnai et al. (2018) using different fractionation equations. Circles: temperature offset for each sample; vertical transparent bars: temperature offset 95 % confidence interval for each sample; horizontal bar: temperature offset arithmetic mean and associated 95 % confidence interval of the mean. (a) Temperature offset for isotopic temperatures calculated with δ¹⁸O_sw values calculated from the LeGrande and Schmidt (2006) database; n=18. (b) Temperature offset for isotopic temperatures calculated with locally measured δ¹⁸O_sw values; n=13. OLS refers to ordinary least-squares regression; YLS refers to York least-squares regression (York et al., 2004).

Independent from the fractionation equation or the δ¹⁸O_sw value used, temperature estimates are associated with large scatter of more than 4 ^∘C (2σ) around mean ΔT. We suggest that this scatter is related primarily to the variability of brachiopod shell δ¹⁸O values observed within a population (Fig. 2; Brand et al., 2019) due to isotopic fractionation mechanisms other than temperature-dependent fractionation, namely kinetic effects, metabolic effects and pH effects. Using the MgCO₃-corrected equation of Brand et al. (2013) instead of that of Brand et al. (2019), both derived from very similar datasets, has the only effect of slightly improving the precision of the estimates, decreasing the scatter by ∼ 0.6 ^∘C around the mean deviation from environmental temperatures. This effect is about 1 order of magnitude below the remaining dispersal of data.

For the different equations derived from modern-brachiopod data, we push further the comparison of the deviation between isotopic temperature calculated with the different equations and environmental temperatures and look at the distribution of these offsets (ΔT) against seawater temperature (Fig. 4). This highlights a significant correlation between ΔT and seawater temperature when applying both the equations of Brand et al. (2013, 2019) to the comparative dataset of Bajnai et al. (2018; p_uncorrelated<0.001; Fig. 4). On contrary, there is no significant correlation between ΔT and seawater temperature when applying both our new equations to the comparative dataset of Bajnai et al. (2018; p_uncorrelated>0.3). This comparison highlights that the application of the equations of Brand et al. (2013, 2019) does not completely account for temperature as an explicative variable for brachiopod Δ¹⁸O_c-w in the dataset of Bajnai et al. (2018).

Figure 4 Distribution of isotopic temperature deviation from living temperature against living temperature using the equations of (a) this study, (b) Brand et al. (2013) and (c) Brand et al. (2019) in relation to the dataset of Bajnai et al. (2018).

In the following, we further assess the impact of trace element incorporation in brachiopod calcite on the fractionation factor between calcite and water using our new data. Brachiopod calcite can be viewed as a three-component molar mixture of CaCO₃, MgCO₃ and SrCO₃, and the fractionation factor of this mixture (Δ_{brach.calcite–water}) of divalent metal carbonates is equivalent to the sum of the mineral–water fractionation factors weighted from the mineral molar fractions. This bulk isotopic fractionation factor is thus expressed as follows:

$$\begin{array}{}\text{(4)}& {\mathrm{\Delta }}_{\text{(brach.calcite–water)}}=\sum _{i=\mathrm{1}}^n{X}_i\cdot {\mathrm{\Delta }}_{\text{(mineral–water)}},\end{array}$$

with the n (=3) minerals noted as i and X being the molar fraction of each mineral constituting the brachiopod calcite mixture. As, to our knowledge, a MgCO₃–water oxygen isotope fractionation equation has not yet been established from synthetic carbonates, we used the theoretical fractionation factors (α_(x−y)) of Chacko and Deines (2008) to compare the mixture and pure calcite fractionation factors. Δ_{(brach.calcite–water)} values, calculated for brachiopod inner calcite layers with the same range of CaCO₃, MgCO₃ and SrCO₃ contents as that recorded by our elemental data, do not exceed 0.18 ‰ compared to the calcite end member at environmental temperatures (mean = 0.08 ‰; n=24). This difference is mainly explained by Mg²⁺ incorporation, as (1) MgCO₃ is more abundant than SrCO₃ in brachiopod calcite; and (2) $\mathrm{1000}\ln {\mathit{\alpha }}_{\left({\mathrm{MgCO}}_{\mathrm{3}}\text{–water}\right)}$ is higher than $\mathrm{1000}\ln {\mathit{\alpha }}_{\left({\mathrm{CaCO}}_{\mathrm{3}}\text{–water}\right)}$ by about 13 ‰ in this temperature range, whereas $\mathrm{1000}\ln {\mathit{\alpha }}_{\left({\mathrm{SrCO}}_{\mathrm{3}}\text{–water}\right)}$ is lower than $\mathrm{1000}\ln {\mathit{\alpha }}_{\left({\mathrm{CaCO}}_{\mathrm{3}}\text{–water}\right)}$ by only ∼ 1 ‰ (Chacko and Deines, 2008). This basic calculation illustrates that Mg²⁺ incorporation into brachiopod shells with low Mg calcite can affect calcite δ¹⁸O by no more than 0.2 ‰, which is about twice the range of analytical uncertainties at 2σ. As a consequence, the incorporation of Mg²⁺ has a weak effect on the δ¹⁸O of brachiopod calcite. The MgCO₃ correction proposed by Brand et al. (2013) is not really needed for paleoclimatic and paleoceanographic studies, which most likely justifies why it was later abandoned (Brand et al., 2019).

4.1.2 Brachiopod shell clumped isotopes: an alternative paleothermometer?

Although there is a strong correlation between our new Δ₄₇ values from brachiopod shells and growing temperatures, Δ₄₇ values at a given temperature are generally higher than what is expected from the state-of-the-art equation of Anderson et al. (2021; Fig. 5). This observation is in line with the results of Bajnai et al. (2018) gathered in the CDES reference frame. To discuss modern-brachiopod Δ₄₇ values in relation to a larger dataset, the data from Bajnai et al. (2018) were adjusted to the CDES 90 reference frame for better comparison with our dataset in the I-CDES reference frame (Bernasconi et al., 2021). We acknowledge that the different standardization protocols used for the two datasets could hamper the strength of this comparison. However, this comparison is supported by the good agreement between sample 143 from Bajnai et al. (2018; ${\mathrm{\Delta }}_{\mathrm{47}}=\mathrm{0.670}\pm \mathrm{0.008}$ ‰ CDES 90) and its replicate sample Mv143b (${\mathrm{\Delta }}_{\mathrm{47}}=\mathrm{0.664}\pm \mathrm{0.006}$ ‰ I-CDES; Bajnai et al., 2020; Fiebig et al., 2021). The Δ₄₇–temperature relationship derived from this combined brachiopod dataset yields the following equation:

$$\begin{array}{}\text{(5)}& {\mathrm{\Delta }}_{\mathrm{47}\left({\text{I-CDES90}}^{\circ }\mathrm{C}\right)}=\mathrm{0.042}\pm \mathrm{0.002}\times {\displaystyle \frac{{\mathrm{10}}^{\mathrm{6}}}{T^{\mathrm{2}}}}+\mathrm{0.124}\pm \mathrm{0.016},\end{array}$$

with ${\mathrm{\Delta }}_{\mathrm{47}\left({\text{I-CDES90}}^{\circ }\mathrm{C}\right)}$ in ‰ and T being the temperature in K. This equation, while still nearly within the error of the equation of Anderson et al. (2021; Fig. 5b), highlights systematically higher brachiopod Δ₄₇ values relative to what is expected from the canonical clumped isotope equation.

Figure 5 Comparison of Δ₄₇ values from modern brachiopods with the state-of-the-art Δ₄₇ calibration of Anderson et al. (2021). Envelope and error bars are at 2σ. (a) Comparison of York least-squares regressions (York et al., 2004) derived from two modern-brachiopod datasets (Bajnai et al., 2018; this study) with the equation of Anderson et al. (2021). (b) Comparison of York least-squares regression derived from the combined modern-brachiopod dataset with the equation of Anderson et al. (2021). (c) Comparison of modern-brachiopod Δ₄₇ temperatures calculated using the equation of Anderson et al. (2021) with seawater temperatures.

A similar underestimation of Δ₄₇ temperatures using the equation of Anderson et al. (2021) was recently reported on culture-grown bivalves (de Winter et al., 2022). The authors argued that the Δ₄₇−1–T² relationship is non-linear, implying that equations including warm (>100 ^∘C) data points, such as that of Anderson et al. (2021), systematically underestimate temperatures in the 30 to −2 ^∘C range (de Winter et al., 2022). To test whether the equation of Anderson et al. (2021) is responsible for the cold bias reported here, we compare our data with three other equations established in the marine temperature range within the I-CDES reference frame (Meinicke et al., 2021; Huyghe et al., 2022; Peral et al., 2022). While applying the equation of Meinicke et al. (2021) slightly reduced the mean Δ₄₇–temperature offset relative to brachiopod living temperature ($-\mathrm{2.4}\pm \mathrm{1.2}$ ^∘C; 2 SE, n=37) compared to the equation of Anderson et al. (2021) ($-\mathrm{3.9}\pm \mathrm{1.2}$ ^∘C; 2 SE, n=37), applying of the equation of Peral et al. (2022; $-\mathrm{3.8}\pm \mathrm{1.2}$ ^∘C; 2 SE, n=37) or of Huyghe et al. (2022; $-\mathrm{3.0}\pm \mathrm{1.2}$ ^∘C; 2 SE, n=37) gave similar offsets (see File S3 in the Supplement for a detailed comparison). This comparison shows that equations derived from samples precipitated exclusively in the marine temperature range give very similar offsets and does not support the hypothesis that the offset is due to very warm data points (>30 ^∘C) influencing the slope of the linear regression of Anderson et al. (2021). We also argue that the difference between the equations of Peral et al. (2018, 2022) and Meinicke et al. (2020, 2021), both based on foraminiferal calcite, can be explained by different ways of estimating calcification temperature rather than by different observations, as discussed by Meinicke et al. (2020). At least a fourth of brachiopod calcite samples show statistically significant deviation from previously published equations (9 to 14 of the 37 samples considered here, depending on the equation). We conclude that the Δ₄₇ offsets observed in brachiopod calcite reflect primary deviations from previous observations that should be addressed to prevent temperature biases when using Δ₄₇ values of fossil brachiopods as paleotemperature proxies.

Our new data suggest that, at a given temperature, brachiopod shell Δ₄₇ values are higher than what has been observed among foraminifera (Breitenbach et al., 2018; Peral et al., 2018, 2022; Daëron et al., 2016; Meinicke et al., 2020, 2021), laboratory-precipitated calcites (Jautzy et al., 2020; Anderson et al., 2021; Fiebig et al., 2021), slow-growing cave calcites (Daëron et al., 2019; Anderson et al., 2021; Fiebig et al., 2021) and bivalve calcite (Huyghe et al., 2022). This observation is more relevant in the Bajnai et al. (2018) dataset, with a mean Δ₄₇ temperature deviation from seawater temperatures of $-\mathrm{5.6}\pm \mathrm{1.6}$ ^∘C (2 SE, n=18), than in our new dataset, with a mean Δ₄₇ temperature deviation from seawater temperatures of $-\mathrm{2.4}\pm \mathrm{1.3}$ ^∘C (2 SE, n=19; Fig. 5a). However, we argue that the difference between the two datasets of modern-brachiopod Δ₄₇ may largely result from distinct taxonomic assemblages and/or sampled locations in both datasets. Indeed, our new Δ₄₇ data from New Zealand brachiopods agree very well with previous data from the same brachiopod species and similar locations (Bajnai et al., 2018; see File S3 in the Supplement). The Bajnai et al. (2018) dataset is dominated by brachiopods from the suborder Terebratellidina ($\mathrm{12}/\mathrm{18}$), a group for which the combined dataset shows the most deviation from the Anderson et al. (2021) equation at low (−2 to 5 ^∘C) temperatures but most importantly in the middle (8–15 ^∘C) temperature range (Fig. 5c). This group is associated with a mean Δ₄₇ temperature deviation from seawater temperatures of $-\mathrm{5.6}\pm \mathrm{1.7}$ ^∘C (2 SE, n=19) when applying the Anderson et al. (2021) equation, while the other most represented group in the combined dataset, the suborder Terebratulidina, better fits the canonical equation (Fig. 5c) with a mean Δ₄₇ temperature deviation from seawater temperatures of $-\mathrm{1.8}\pm \mathrm{1.4}$ ^∘C (2 SE, n=13). This apparent taxonomic difference among brachiopod groups suggests that the deviation observed for brachiopod shell Δ₄₇ values relative to other calcifying organisms relates to processes other than ambient temperatures, controlling bound ordering within the calcite crystal of brachiopod shells (Bajnai et al., 2018). Such processes still need to be confidently identified (Sect. 4.2) but likely lie in the different ways the shell is secreted under different growth rate dynamics and/or different crystal arrangements. Available brachiopod Δ₄₇ values suggest that, while the equation of Anderson et al. (2021) could be confidently applied to brachiopods of the suborder Terebratulidina (Fig. 5), this would not be the case for species of the suborder Terebratellidina. For the other brachiopod groups, data are still too scarce to propose an enlightened conclusion.

4.1.3 $\mathrm{Mg}/\mathrm{Ca}$ thermometer for brachiopod shells

Our $\mathrm{Mg}/\mathrm{Ca}$ data obtained from the inner layers of brachiopod shells do not show any significant correlation with growth temperature (p=0.94; Fig. 1c; Table 3), at odds with the strong correlation between the $\mathrm{Mg}/\mathrm{Ca}$ molar ratio of the brachiopod shell secondary layer and the growth temperature found by Brand et al. (2013, 2019). However, these authors acknowledged that some species do not follow the trend they illustrated as the GBMgL and suggested that such deviations resulted from either lower growth rates (Brand et al., 2013) or higher growth rates (Brand et al., 2019) of the considered species relative to other brachiopods. This issue alone illustrates how, from the largest dataset of modern brachiopods available so far (i.e. Brand's database), one cannot firmly determine whether or not $\mathrm{Mg}/\mathrm{Ca}$ correlates with seawater temperatures. We emphasize that our dataset lacks brachiopod shell samples from equatorial surface marine waters (27 ^∘C < T < 32 ^∘C), which record the highest MgCO₃ molar content in Brand's database (Brand et al., 2013, 2019). Interestingly, the highest MgCO₃ contents reported by these authors were measured on brachiopods belonging to the order Thecideida. This taxonomic singularity is of prime importance because (1) Thecideide brachiopod shells are essentially made of a shell fabric unique to its group (Simonet Roda et al., 2022), and secondary shell fabric (fibrous layers, which form the majority of Rhynchonellida and Terebratulida shells) is restricted to small areas of the shell, such as the teeth and the inner sockets (Williams, 1968, 1973; Baker, 2006); (2) Thecideide shell is formed from high-Mg calcite (Ullmann et al., 2017). The exclusion of Thecideida data from the database published by Brand et al. (2013) produces a substantially weaker (R² = 0.23) yet significant (p<0.01) correlation between shell $\mathrm{Mg}/\mathrm{Ca}$ and temperature. In addition, magnesium has also been shown to be more concentrated in major growth bands (Gaspard et al., 2018; Müller et al., 2022) than in other parts of the secondary shell layer in articulated brachiopods, which could explain the high $\mathrm{Mg}/\mathrm{Ca}$ values observed from slow-growth and long-lifespan brachiopods (Brand et al., 2019), especially when considering samples obtained from a large sampling area covering multiple growth lines. Given our new data and all these considerations, the brachiopod $\mathrm{Mg}/\mathrm{Ca}$ should be regarded as an unreliable proxy for annual or seasonal temperature reconstructions, at least in the −2 to 25 ^∘C range (Fig. 1; Brand et al., 2013, 2019).

An alternative or complementary way to explain the scattering of temperature–$\mathrm{Mg}/\mathrm{Ca}$ paired data may arise in the variability of marine $\mathrm{Mg}/\mathrm{Ca}$. All the elements considered in this study (Mg, Sr, Li and Na) have much larger residence times (10, 3.5, 1.8 and 44 Myr, respectively) than water molecules in the oceans, meaning that, at first order, they are uniformly distributed in the main bodies of the oceans (Lécuyer, 2016). A recent study, however, has shown that sizable variations in element $/$ Ca may occur in the vicinity of continental margins under the influence of oceanic currents at the interface with groundwaters or the atmosphere or within water masses isolated from the global thermohaline circulation (e.g. Black Sea, Mediterranean Sea, Red Sea; Lebrato et al., 2020). For example, despite the large residence time of Sr, non-homogeneous marine strontium isotopic ratios (${}^{\mathrm{87}}\mathrm{Sr}{/}^{\mathrm{86}}\mathrm{Sr}$) have been measured in coastal environments (southern Okinawa Trough, South China Sea and Kaoping Canyon) as the result of seawater mixing with groundwaters (Huang et al., 2011; El Meknassi et al., 2020). Given the poor constraint on local coastal water chemistry in the deep time, brachiopod samples from offshore environments should be preferred for trace-element-based paleotemperature reconstructions, as was suggested for marine biogenic apatite by Balter and Lécuyer (2010).

4.1.4 Potential of other trace element $/$ Ca proxy

Our new brachiopod $\mathrm{Sr}/\mathrm{Ca}$ data show a significant but weak negative correlation with ambient temperature (p_slope<0.001; R²<0.40). This correlation is supported by a weak positive correlation between $\mathrm{Sr}/\mathrm{Ca}$ and Δ¹⁸O_c-w of the inner layers (p_slope<0.01; R²=0.31), as we also expect Δ¹⁸O to decrease with increasing temperatures (i.e. Sect. 4.1). This contradicts the pioneering findings of Lowenstam (1961), which suggested a negative correlation between brachiopod shell $\mathrm{Sr}/\mathrm{Ca}$ and Δ¹⁸O_c-w. The brachiopod database of Brand et al. (2013) shows no significant correlation between $\mathrm{Sr}/\mathrm{Ca}$ and temperatures (p_slope=0.20). Given the diverging conclusions from one dataset to another, $\mathrm{Sr}/\mathrm{Ca}$ values of brachiopod shells do not appear to be controlled by growing temperatures at the first order.

$\mathrm{Li}/\mathrm{Ca}$ in brachiopod shells is significantly correlated to temperature (Fig. 1e; Table 3), a covariation initially reported by Delaney et al. (1989) and confirmed by more recent studies (Dellinger et al., 2018; Rollion-Bard et al., 2019; Washington et al., 2020). Most of the data presented here follow the $\mathrm{Li}/\mathrm{Ca}$–temperature trend derived from literature data (Washington et al., 2020), with several samples (n=9), exclusively collected from the inner layer, showing significantly lower $\mathrm{Li}/\mathrm{Ca}$. Those low $\mathrm{Li}/\mathrm{Ca}$ values lie below the $\mathrm{Li}/\mathrm{Ca}$–temperature relationship derived from inorganic calcite (Marriott et al., 2004). Similarly, low $\mathrm{Li}/\mathrm{Ca}$ values from brachiopod shells were reported by previous studies (Delaney et al., 1989; Washington et al., 2020). Interestingly, most of those low $\mathrm{Li}/\mathrm{Ca}$ values belong to samples from the genus Tichosina (Delaney et al., 1989; Washington et al., 2020, Table 2), so one could consider some peculiar species-dependent trace element partitioning driven by metabolic or kinetic processes. Caution is warranted, however, because very low $\mathrm{Li}/\mathrm{Ca}$ values are also apparent in the other tropical genus studied here, Stenosarina Cooper 1977, as well as in one specimen of Aerothyris kerguelenensis among the nine specimens studied there.

Other environmental parameters or biological processes may control the element $/$ Ca incorporation into calcite. $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Li}/\mathrm{Ca}$ and $\mathrm{Na}/\mathrm{Ca}$ show significant negative correlations with the local salinity ($\mathrm{0.37}<R^{\mathrm{2}}<\mathrm{0.76}$; p slope < 0.01), with similar scattering as that for the temperature relationships. $\mathrm{Mg}/\mathrm{Ca}$ ratios, again, show no significant correlation with this parameter (p slope = 0.7). A wider look at the geochemical dataset highlights very strong correlations of $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Li}/\mathrm{Ca}$, and $\mathrm{Na}/\mathrm{Ca}$ ratios in brachiopod shells with one another ($\mathrm{0.68}<R^{\mathrm{2}}<\mathrm{0.92}$; p slope < 0.01) over the whole dataset and particularly from inner-layer samples ($\mathrm{0.87}<R^{\mathrm{2}}<\mathrm{0.97}$; p slope < 0.01). Correlation in our dataset is strongest between $\mathrm{Sr}/\mathrm{Ca}$ and $\mathrm{Na}/\mathrm{Ca}$ and weakest between $\mathrm{Sr}/\mathrm{Ca}$ and $\mathrm{Li}/\mathrm{Ca}$. Similar correlation trends also appear at the intra-specimen level (Romanin et al., 2018; Rollion-Bard et al., 2019). Such strong correlations, at both intra-individual and inter-specific scales, may suggest that the first-order variability of those different element concentrations within brachiopod shells may be explained by similar processes, either environmental or biological. Ullmann et al. (2017) reported a distinct geochemical signature between brachiopods species, especially within the $\mathrm{Mg}/\mathrm{Ca}$–$\mathrm{Sr}/\mathrm{Ca}$ space. Our inner-shell-layer dataset presents a similar distribution within the low-$\mathrm{Mg}/\mathrm{Ca}$ range (<10 mmol $/$ mol), where $\mathrm{Sr}/\mathrm{Ca}$ ratios of Terebratulidina shells are distinctly lower than those of Terebratellidina shells. Owing to the strong correlation between $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Li}/\mathrm{Ca}$ and $\mathrm{Na}/\mathrm{Ca}$, our dataset differentiates between Terebratulidina shells with low element $/$ Ca ratios ($\mathrm{0.4}<\mathrm{Sr}/\mathrm{Ca}<\mathrm{0.7}$ mmol $/$ mol; $\mathrm{0}<\mathrm{Li}/\mathrm{Ca}<\mathrm{10}$ µmol $/$ mol; $\mathrm{1}<\mathrm{Na}/\mathrm{Ca}<\mathrm{4}$ mmol $/$ mol) and Terebratellidina shells with high element $/$ Ca ratios ($\mathrm{0.9}<\mathrm{Sr}/\mathrm{Ca}<\mathrm{1.2}$ mmol $/$ mol; $\mathrm{35}<\mathrm{Li}/\mathrm{Ca}<\mathrm{60}$ µmol $/$ mol; $\mathrm{8}<\mathrm{Na}/\mathrm{Ca}<\mathrm{12}$ mmol $/$ mol), with the exception of the Antarctic brachiopod Liothyrella uva (Broderip, 1833; Terebratulidina) which shows the highest $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Na}/\mathrm{Ca}$ and $\mathrm{Li}/\mathrm{Ca}$ values (Fig. 6a–b). With this exception and that reported by Ullmann et al. (2017), a taxonomic relationship is not satisfactory, especially considering the significant yet noisy relationship between element $/$ Ca and environmental parameters such as temperature and salinity within our new dataset. Indeed, we also observe a strong relationship between brachiopod shell element $/$ Ca values and geographical location (Fig. 6c–d). Tropical shells (Guadeloupe and New Caledonia – only Terebratulidina) have low element $/$ Ca values, the middle latitudes (New Zealand and Crozet Islands – dominated by Terebratellidina) have intermediate to high element $/$ Ca values, and polar brachiopods (Antarctica and Norway – mixed assemblage) have high element $/$ Ca values. Despite this first-order distribution, it appears that, within the same environment, species difference appears to dominate, as is the case in Doubtful Sound (New Zealand), where Terebratella sanguinea (Leach, 1814; Terebratellidina) shells have higher element $/$ Ca values than Liothyrella neozelanica (Thomson, 1918; Terebratulidina) shells.

Figure 6 Distribution of trace element content of the inner shell layers among the different brachiopod taxonomic groups (same symbols as in Fig. 5c). The data show a strong correlation between $\mathrm{Sr}/\mathrm{Ca}$ and $\mathrm{Li}/\mathrm{Ca}$ (a) and $\mathrm{Sr}/\mathrm{Ca}$ and $\mathrm{Na}/\mathrm{Ca}$ (b), as well as strong latitudinal partitioning of $\mathrm{Li}/\mathrm{Ca}$ (c) and $\mathrm{Sr}/\mathrm{Ca}$ (d) values.

Our new data suggest that temperature is not the prime parameter that can describe trace element incorporation into articulated brachiopod shell calcite from a dataset of widely distributed and taxonomically diverse brachiopods. This conclusion does not exclude the possibility of any temperature control on trace element incorporation in brachiopod shell calcite, which may be relevant at the species level, as proposed by Butler et al. (2015) and observed in other calcifying organisms such as bivalves (Freitas et al., 2006) and foraminifera (Anand et al., 2003). However, the influence of other controls on trace element incorporation in brachiopod shell calcite (growth rate, environment fluid chemistry, etc.; Jurikova et al., 2020) is strong enough to blur the likely temperature control so that the trace element $/$ Ca ratios studied here are not suited to being used as paleothermometers.

4.2 Precipitation of brachiopod shell calcite out of equilibrium with seawater

Our new data confirm that brachiopod shell calcite does not precipitate in isotopic equilibrium with seawater, as highlighted by the significant deviation from thermodynamic isotopic equilibrium relationships constrained from very-slow-growing calcite (Daëron et al., 2019; Anderson et al., 2021, Figs. 2–5). It is worth noting that most brachiopod shell δ¹⁸O data (Brand et al., 2019; this study) plot in between thermodynamic equilibrium isotopic fractionation (Daëron et al., 2019) and the equation of Kim and O'Neil (1997; Fig. 2). The deviation of brachiopod shell calcite δ¹⁸O from thermodynamic equilibrium (Δ¹⁸O_eq = δ¹⁸O${}_{\mathrm{brach}.\mathrm{calcite}}-{\mathit{\delta }}^{\mathrm{18}}$O_equilibrium) is obvious when we compare our dataset to the equation of Daëron et al. (2019) in the clumped isotope–oxygen isotope fractionation space (Fig. 7a), with most of our dataset showing higher Δ₄₇ and lower δ¹⁸O than predicted (Fig. 7b), a feature previously highlighted by Bajnai et al. (2018). As a result, brachiopod calcite shows apparent Δ₄₇ and δ¹⁸O temperatures that are, respectively, lower and higher than those predicted by the thermodynamic equilibrium. This deviation from equilibrium values has a Δ₄₇–Δ¹⁸O_eq slope of 0.019±0.005 in our new brachiopod dataset, indistinguishable from the slope of 0.017±0.003 reported by Bajnai et al. (2018) using δ¹⁸O_sw values from the δ¹⁸O_sw database (LeGrande and Schmidt, 2006).

Figure 7 Non-equilibrium fractionation and kinetic effects in brachiopod shell calcite as evidenced by clumped isotopes and oxygen isotopes. (a) Plot of brachiopod inner shell layers in the Δ₄₇–1000ln α_{(calcite–water)} space compared with the equilibrium line (blue line) defined by very-slow-growing abiogenic calcite (Daëron et al., 2019; Anderson et al., 2021). The red highlight shows the results from A. kerguelenensis (Crozet Islands) and their expected geochemical values on the equilibrium line. Note the negative correlation between Δ₄₇ and 1000ln α_{(calcite–water)} evidenced in this subset. (b) δ¹⁸O_c and Δ₄₇ offset from their expected equilibrium values calculated from environmental parameters using the equations of Daëron et al. (2019) and Anderson et al. (2021), respectively. The trends followed by these offsets within our new dataset are similar to those reported from other brachiopod samples by Bajnai et al. (2018).

Brachiopod shell Δ¹⁸O data (Brand et al., 2003, 2013, 2019; this study) also show a large geochemical variability among brachiopods sampled from similar marine environments (T, S), with a Δ¹⁸O range often higher than 1 ‰ (Fig. 3; Brand et al., 2019). Isotopic fractionation equations (Brand et al., 2013, 2019; this study) usually rely on Δ¹⁸O_c-w values plotted against mean annual temperatures. Thus, it would be tempting to attribute such scattering of data to seasonal variations in both temperature and δ¹⁸O_sw, but this would only apply to specimens that lived above or in the upper part of the thermocline. Such isotopic scattering is indeed expected if we consider datasets based on multiple sampling along the growth axis of an individual (Yamamoto et al., 2011; Takayanagi et al., 2015), such as in the case of Brand et al. (2013, 2019). However, seasonal scattering in data resulting from intra-individual geochemical variability is unlikely for our measurements, which are performed on a large sampling area that should correspond to multiple months or years of growth.

For example, we note a large isotopic variability in both δ¹⁸O and Δ₄₇ among the nine specimens of Aerothyris kerguelenensis from the Crozet Islands even though the specimens were recovered from depths in the 105 to 355 m range with a narrow range of temperature (3–5 ^∘C) and salinity (33.8–34.2). These nine samples, all collected from the shell inner layers of different individuals, record a Δ¹⁸O variability of 2 ‰, which is equivalent to a temperature range of 8 to 10 ^∘C if interpreted as such, and a Δ₄₇ variability of 0.02 ‰, which is equivalent to a temperature range of about 5 ^∘C. This subsample shows a significant Δ₄₇–Δ¹⁸O_eq slope of 0.027±0.007 (2σ; p_slope<0.01; R²=0.95, n=5), as well as a strong positive correlation between Δ¹⁸O_c-w and δ¹³C (p_slope<0.01; R²=0.69, n=9). These covariations most likely result from kinetic effects as defined by McConnaughey (1989), here recorded at the intra-specific level.

Our new data replicate the observations made by Bajnai et al. (2018) concerning kinetic fractionation in brachiopod shells. In consequence, we follow their conclusions and consider it to be the case that different growth rates of the specimens are responsible for the second-order isotopic variability and the deviation from the stable isotope and clumped isotope equilibrium observed in brachiopod datasets, both at an inter-specific and intra-specific level.

Neither minor nor trace elements are incorporated into brachiopod shell calcite in equilibrium with seawater relative to abiogenic calcite (Rollion-Bard et al., 2019; Simonet Roda et al., 2022). Kinetic effects are known to influence trace element incorporation into calcite (Lorens, 1981; Busenberg and Plummer, 1985; Tesoriero and Pankow, 1996; Gabitov et al., 2011, 2014) and may be a process, among others (Jurikova et al., 2020), that explains trace element incorporation into brachiopod shell calcite. Owing to the very limited knowledge of brachiopod shell growth rates (for only 5 out of 17 species in our dataset), we are unable to directly test this hypothesis. However, if kinetic effects have such an influence on element incorporation into brachiopod shell calcite, we should observe strong negative correlations of $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Na}/\mathrm{Ca}$ and $\mathrm{Li}/\mathrm{Ca}$ ratios with Δ¹⁸O_eq. We observe such significant correlations in the new dataset in both the inner- and outer-layer subsets ($\mathrm{0.40}<R^{\mathrm{2}}<\mathrm{0.55}$, p slope < 0.01), with a significant slope change between the two subsets. The variation in element $/$ Ca related to a similar variation in Δ¹⁸O_eq, as suggested by linear regression models, is about twice as high in the inner-layer subset than in the outer-layer subset. This difference may be explained by the different pathways brachiopods use to form their different shell layers (acicular, fibrous and prismatic/columnar; Simonet Roda et al., 2022). However, while looking at Aerothyris kerguelenensis, which showed the most important kinetic effects according to stable isotope data, there is no significant correlation of element $/$ Ca ratios with Δ¹⁸O_eq. Element $/$ Ca ratios appear to be independent from isotope kinetic effects at the species level despite strong correlations of element $/$ Ca ratios and Δ¹⁸O_eq in the multi-specific dataset.

In line with the significant correlation between Δ¹⁸O_eq and element $/$ Ca values in brachiopod shell inner layers, the environmental and taxonomic distributions of isotopic offsets (Δ¹⁸O_eq, Δ₄₇ offset) in our dataset are similar to the distribution of element $/$ Ca values. We observe more deviation from the equilibrium in high-latitude environments than in tropical environments, and we generally observe more deviation from isotopic equilibrium in Terebratellidina than in Terebratulidina in the whole dataset, although the difference is not as striking as that documented for trace elements. As the geographical distribution of the two suborders in the new dataset is highly biased, this trend cannot be confidently interpreted as arising from different shell formation processes in the different groups. However, several examples of Terebratellidine and Terebratulidine living in the same environment show a similar pattern. In Doubtful Sound (New Zealand), Terebratella sanguinea (Terebratellidina) deviates more from equilibrium isotopic composition than Liothyrella neozelanica (Terebratulidina; this study; Bajnai et al., 2018). In Friday Harbor (Washington, USA), Terebratalia transversa (Sowerby, 1846; Terebratellidina) records lower and more variable δ¹⁸O_c than Terabratulina sp. (Terebratulidina; Ullmann et al., 2017). These authors also reported high isotopic variability associated with the strong correlation of δ¹⁸O with δ¹³C among Terebratellidine but not among Terebratulidine.

From the relationship between the different parameters studied here, we suggest that kinetic effects in brachiopod shells may be a major source of deviation from calcite precipitation in equilibrium with seawater. Kinetic effects can be, on first approximation, controlled by shell growth rates (McConnaughey, 1989; Bajnai et al., 2018). Therefore, both biological and environmental parameters should control the extent of the kinetic effects. At first order, one can expect species growing larger shells to record more important kinetic effects. Moreover, shell growth rates in calcifying organisms increase when environmental conditions are close to optimal living range and decrease when the conditions strongly derive from that optimum, a pattern that is dependent on taxon-related preferences (Schöne, 2008). However, the poor knowledge of brachiopod shell growth rate, relative to other marine calcifying organisms, prevents an in-depth appreciation of its exact role. Kinetic effects are still insufficient to describe all of the deviation from equilibrium observed here, but our dataset does not allow for a discussion of other potential factors (Jurikova et al., 2020).

5 Conclusions

Our multi-proxy (δ¹⁸O, δ¹³C, Δ₄₇, $\mathrm{Mg}/\mathrm{Ca}$, $\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Na}/\mathrm{Ca}$ and $\mathrm{Li}/\mathrm{Ca}$) geochemical analyses of 37 modern brachiopods from different locations, spanning living temperatures from −1.9 to 25.5 ^∘C, allowed us to investigate the temperature relationship of several geochemical markers and to evaluate or re-evaluate their potential as paleothermometers for brachiopod shell archives. δ¹⁸O and Δ₄₇ are the most reliable paleothermometers, although we propose a revised δ¹⁸O paleotemperature equation for brachiopod shell calcite, and Δ₄₇ values highlight a significant deviation from the canonical temperature relationship. The element $/$ Ca ratios investigated here should not be used to infer environmental temperatures from brachiopod calcite owing to either the absence ($\mathrm{Mg}/\mathrm{Ca}$) or a poor yet significant ($\mathrm{Sr}/\mathrm{Ca}$, $\mathrm{Li}/\mathrm{Ca}$ and $\mathrm{Na}/\mathrm{Ca}$) correlation with living temperature. Our new dataset confirms previous evidence for the interpretation of non-equilibrium isotopic fractionation interpreted as significant kinetic effects recorded in the brachiopod shell inner layers. In general, the precipitation of brachiopod shells that takes place out of isotopic equilibrium leads to (1) an underestimation of living temperature from Δ₄₇ and (2) an overestimation of living temperature from oxygen isotopes, which, once the proxies are combined, could lead to (3) an underestimation of living-water δ¹⁸O_sw relative to what would be expected from equilibrium precipitation. Kinetic effects and, by extension, shell growth rates could explain the second order of isotopic variability in brachiopod shells and can be a larger source of isotopic variability than environmental parameters (temperature, δ¹⁸O_sw) at the local scale. Thus, we advise that this variability be considered as a significant source of uncertainties when applying paleotemperature proxies. Although we reject their use as a paleothermometer, we do not exclude any temperature control on trace element (Mg, Sr, Li and Na) incorporation into brachiopod shell calcite. Rather, our findings indicate that other processes, such as kinetic effects (growth rate) or shell microstructure and their associated mineralization pathways, have a significant control over trace element incorporation into brachiopod shell calcite, resulting in the high scrambling observed in the element $/$ Ca–temperature relationship.