Exploring macroevolutionary links in multi-species planktonic foraminiferal Mg∕Ca and δ¹⁸O from 15 Ma to recent

The ratio of the trace element Mg over Ca ($\mathrm{Mg}/\mathrm{Ca}$) and the oxygen isotopic composition (δ¹⁸O) of foraminiferal calcite are widely employed for reconstructing past ocean temperatures, although geochemical signals are also influenced by several other factors that vary temporally and spatially. Here, we analyse a global dataset of $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O data of 59 middle Miocene to recent species of planktonic foraminifera from a wide range of depth habitats, many of which have never been analysed before for $\mathrm{Mg}/\mathrm{Ca}$. We investigate the extent to which $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O covary through time and space and identify several sources of mismatch between the two proxies. Once the data are adjusted for long-term non-thermal factors, $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O are overall positively correlated in a way consistent with temperature being the dominant controller through both space and time and across many different species, including deep dwellers. However, we identify several species with systematic offsets in $\mathrm{Mg}/\mathrm{Ca}$ values, to which multispecies calibrations should be applied with caution. We can track the appearance of such offsets through ancestor-descendent species over the last 15 Myr and propose that the emergence of these offsets may be the geochemical expression of evolutionary innovations. We find that virtually all of the $\mathrm{Mg}/\mathrm{Ca}$- and δ¹⁸O-derived temperatures from the commonly used genera Globigerinoides and Trilobatus are within uncertainty of each other, highlighting the utility of these species for paleoceanographic reconstructions. Our results highlight the potential of leveraging information from species lineages to improve sea surface temperature reconstruction from planktonic foraminifera over the Cenozoic.

1 Introduction

Geochemical analyses of foraminifera are commonly applied to reconstruct paleoceanographic conditions, such as marine temperatures, and therefore infer past climatic changes. In particular, the fossil tests of planktonic foraminifera (calcareous zooplankton) provide one of the most widely used paleoclimate archives. Here we focus on two of these parameters: δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$, both of which have been used widely as temperature proxies.

The oldest and possibly most widely utilized of these proxies is the ratio of oxygen isotopes in the calcite tests of these organisms, which, due to slight differences in reactivity of molecules containing the different isotopes, is temperature-dependent (Urey, 1947; see Pearson, 2012, for review). This effect has been quantified in experiments with inorganic calcite (e.g. Kim and O'Neil, 1997) and planktonic foraminifera in culture (e.g. Erez and Luz, 1983; Bemis et al., 1998). Tests of planktonic foraminifera calcifying in warmer waters are depleted in ¹⁸O relative to species living in cooler waters (Emiliani, 1954). A second, more recently established paleoclimate proxy is the ratio of magnesium to calcium in test calcite (Chave, 1954; Nürnberg et al., 1996). Inorganic precipitation experiments have shown that the $\mathrm{Mg}/\mathrm{Ca}$ ratio of calcite is found to be higher at greater temperatures (Mucci, 1987). This relationship led to the in-depth exploration of $\mathrm{Mg}/\mathrm{Ca}$ ratios in planktonic and benthic foraminifera and its potential application as a temperature proxy through culturing (Lea et al., 1999; von Langen et al., 2005), core top (Nürnberg, 1995; Elderfield and Gassen, 2000), and sediment trap studies (Anand et al., 2003).

As they represent two different chemical systems, the $\mathrm{Mg}/\mathrm{Ca}$ and oxygen stable isotope ratios in foraminifera are often used together as independent temperature proxies. For instance, δ¹⁸O-derived calcification temperatures have been combined with $\mathrm{Mg}/\mathrm{Ca}$ data to derive $\mathrm{Mg}/\mathrm{Ca}$ temperature calibrations (e.g. Anand et al., 2003; McConnell and Thunell, 2005; Mohtadi et al., 2009). Other studies have applied these two systems together to infer the influence of environmental parameters such as seawater salinity on $\mathrm{Mg}/\mathrm{Ca}$ (e.g. Mathien-Blard and Bassinot, 2009; Hönisch et al., 2013) and global ice volume (e.g. Lear et al., 2000; Katz et al., 2008). Works such as these rely on temperature being a primary driver of the two proxies for any given sample, which should be the case if both systems are impacted purely by calcification temperature. Nonetheless, there are known non-thermal effects influencing both $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O. For oxygen isotope values, these include the oxygen isotopic composition of seawater ($\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$) and, to a lesser degree, seawater pH or carbonate ion concentration (Spero et al., 1997; Zeebe, 1999). Seawater carbonate chemistry has also been shown to impact the $\mathrm{Mg}/\mathrm{Ca}$ proxy. Culture and sediment trap studies demonstrate that surface ocean seawater pH can influence $\mathrm{Mg}/\mathrm{Ca}$ in planktonic foraminifera (Lea et al., 1999; Evans et al., 2016b; Gray et al., 2018), with the sensitivity of $\mathrm{Mg}/\mathrm{Ca}$ to pH appearing to vary between species (Gray and Evans, 2019). $\mathrm{Mg}/\mathrm{Ca}$ values of foraminifera are also dependent on the $\mathrm{Mg}/\mathrm{Ca}$ of seawater (Evans et al., 2016a), and both oxygen isotope and $\mathrm{Mg}/\mathrm{Ca}$ values can be impacted by test recrystallization (Dekens et al., 2002). $\mathrm{Mg}/\mathrm{Ca}$ values are susceptible to the preferential loss of Mg during dissolution and are thus influenced by the calcite saturation state of bottom waters (Regenberg et al., 2014; Tierney et al., 2019). Seawater salinity has a minor secondary effect on $\mathrm{Mg}/\mathrm{Ca}$ values (Kisakürek et al., 2008; Hönisch et al., 2013), and whilst salinity has little direct effect on oxygen isotopes, a change in salinity is usually accompanied by a change in $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ because hydrological processes such as evaporation and precipitation drive local $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ as well as salinity (LeGrande and Schmidt, 2006). Lastly, so-called “vital effects”, which lump together a wide variety of species-specific processes such as metabolism (including the process of calcification and the incorporation of metabolic products), the position within the water column and life cycle depth migration, the presence of photosymbionts, and seasonality (see summary in Schiebel and Hemleben, 2017), also add complexity to the interpretation of both the oxygen isotope and $\mathrm{Mg}/\mathrm{Ca}$ proxies.

Here we use the dataset published in Boscolo-Galazzo et al. (2021) to examine covariance between $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O in planktonic foraminifera extracted from sediments across a wide range of geographic locations, time intervals, and species. The dataset is composed of δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ data measured on 59 species of planktonic foraminifera, of which 24 have never before been measured for $\mathrm{Mg}/\mathrm{Ca}$ (Tables S1 and S2 in the Supplement). These data are from different ocean basins and latitudes and a range of ages between the middle Miocene (∼15 Ma) and the Holocene. Paired $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O were measured on the same samples; hence this dataset is ideally suited to isolate potential ecological, environmental and preservational factors, which may imprint $\mathrm{Mg}/\mathrm{Ca}$ or δ¹⁸O or both and which are otherwise impossible to recognize in studies focusing on a limited number of species and/or a narrow study area or time interval. In particular, it provides the unique opportunity to simultaneously do the following:

1. compare coupled δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ data on a geographical and temporal scale that is broader than usual;

2. compare coupled δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ data across species of different ecologies;

3. evaluate $\mathrm{Mg}/\mathrm{Ca}$ data of extinct species against those of their modern descendants;

4. test whether temperature can still be recognized as predominantly driving covariance in the dataset when spatial, temporal, and ecological variables are simultaneously in play.

2 Material and analytical methods

2.1 Material

The dataset (Boscolo-Galazzo et al., 2021) was produced from a range of globally and latitudinally distributed DSDP (Deep Sea Drilling Program), ODP (Ocean Drilling Program), and IODP (Integrated Ocean Drilling Program/International Ocean Discovery Program) sites (Fig. 1) which have a high CaCO₃ content and composed of calcareous nannofossils and foraminiferal pelagic oozes, with some input of siliceous plankton and clay. Sites were selected based on the best available global and temporal coverage and preservation of foraminifera. Planktonic foraminiferal preservation ranges from excellent to very good (recrystallized but lacking overgrowth and infilling; Boscolo-Galazzo et al., 2021), with the exception of Sites U1490 and U1489, where there is some overgrowth and infilling in the middle Miocene (Fayolle and Wade, 2020; Boscolo-Galazzo et al., 2021). The target ages selected for sampling were 0, 2.5, 4.5, 7.5, 10, 12.5, and 15 Ma. Biostratigraphic analysis was used to assess age using the biochronology of Wade et al. (2011) calibrated to the timescale of Lourens et al. (2004) (Table S1).

Figure 1 Site map with present-day mean annual sea surface temperatures (°C) from the World Ocean Atlas 2013 (Locarnini et al., 2013).

2.2 Planktonic foraminifera

A total of 59 species of planktonic foraminifera were analysed for $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O. Planktonic foraminifera were picked from three constrained size fractions: 180–250 µm, 250–300 µm, and 300–355 µm. Planktonic foraminiferal geochemistry can change as a function of size (e.g. Birch et al., 2013), such that we used data from the size fraction 250–355 µm only, giving a total of 57 species in the dataset. For abundant species, up to 80 specimens were picked for geochemical analysis, with as many as possible picked in the case of less common species. Hence, our foraminiferal data represent an average from multiple specimens. Paleodepth habitat attributions follow Boscolo-Galazzo et al. (2021) and Boscolo-Galazzo et al. (2022). Planktonic foraminiferal taxonomy follows the concepts described in Boscolo-Galazzo et al. (2022).

2.3 Trace element and stable isotope analysis

Picked planktonic foraminifera were crushed between two glass slides to open all large chambers. When there was enough material, the crushed sample was split for stable isotope and trace element analysis. The trace element split was cleaned using a protocol to remove clays and organic matter (step A1.1–A1.3 of Barker et al., 2003). The samples did not undergo reductive cleaning due to their fragility and small sample size and because the reductive step may cause preferential removal of high $\mathrm{Mg}/\mathrm{Ca}$ calcite from the test (Yu et al., 2007). Samples were dissolved in trace metal pure 0.065 M HNO₃ and then diluted with trace metal pure 0.5 M HNO₃ and analysed at Cardiff University on a Thermo Fisher Scientific Element XR ICP-MS against standards with matched calcium concentration to reduce matrix effects (Lear et al., 2002). Long-term analytical precision determined from consistency standards (CS1 and CS2) with $\mathrm{Mg}/\mathrm{Ca}$ ratios of 1.24 and 7.15 mmol mol⁻¹ is ∼0.7 % and ∼0.8 % (relative standard deviation). $\mathrm{Mg}/\mathrm{Ca}$ was plotted against $\mathrm{Fe}/\mathrm{Ca}$ and $\mathrm{Mn}/\mathrm{Ca}$ to assess whether there was any relationship as a result of the presence of Fe–Mn oxyhydroxides affecting $\mathrm{Mg}/\mathrm{Ca}$, but there was no correlation between the contaminant indicators and $\mathrm{Mg}/\mathrm{Ca}$ (Fig. S1).

Stable isotopes were measured on a Delta V Advantage with a GasBench II mass spectrometer at the Cardiff University stable isotope facility. Stable isotope results were calibrated to the VPDB scale using an in-house carbonate standard (Carrara marble). Analytical precision was 0.05 ‰ for δ¹⁸O and 0.05 ‰ for δ¹³C.

2.4 Data analysis

Before performing the analysis, we screened the dataset for outliers and removed one anomalously high data point with a $\mathrm{Mg}/\mathrm{Ca}$ value $>\mathrm{9}\mathrm{mmol}{\mathrm{mol}}^{-\mathrm{1}}$, which we attributed to analytical error (Table S1).

2.4.1 Formulation of theoretical relationships between $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O

To test for covariation between $\mathrm{Mg}/\mathrm{Ca}$ and oxygen isotope data, we regressed these data against each other and compared the observed relationship with that expected from modern calibrations. We did so to initially explore the relationship we might expect between $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O and whether this manifests in the dataset, before applying corrections for the non-thermal influences on both proxies. Given the complexity of the sample set (e.g. multiple species, ages, locations, preservation), different expected relationships between $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O are possible, which depend on the following: (i) species-specific vital effects, (ii) non-thermal controls on $\mathrm{Mg}/\mathrm{Ca}$ (salinity, pH, $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$), and (iii) non-thermal controls on δ¹⁸O ($\mathrm{pH}/\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$), $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$, as well as how these factors change through time. To account for this, we calculated a number of possible expected theoretical relationships to give a sense of how much of the scatter in the raw data is likely to be explicable by these factors and inform our following data analysis accordingly. We stress that this exercise was conducted as an initial exploration of the dataset; no single relationship will be able to explain the dataset because it is influenced by multiple, often interlinked, variables.

Expected theoretical relationships were calculated starting with modern laboratory culture calibrations, onto which the key non-thermal long-term and spatial controls on these proxies were sequentially added to demonstrate how much each of these is expected to shift the slope of the expected $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O relationship (Fig. 2a). Specifically, we did the following:

• i.

  combined the $\mathrm{Mg}/\mathrm{Ca}$ temperature calibrations for Globigerinoides ruber and Trilobatus sacculifer of Gray and Evans (2019) with the δ¹⁸O–temperature equation of Erez and Luz (1983);

• ii.

  added the impact of a 0.15 unit whole-ocean pH change (approximating the magnitude of the Neogene whole-ocean change, e.g. Rae et al., 2021) using the pH–$\mathrm{Mg}/\mathrm{Ca}$ slope for G. ruber as an example (note that this is only applicable to species that show a pH sensitivity) (Evans et al., 2016b);

• iii.

  included the expected control of temperature on pH via the T-dependent dissociation of water (K_W), i.e. temperature-driven pH changes within a given time interval independent of whole-ocean pH shifts (Gray et al., 2018);

• iv.

  showed the impact of $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$ half of the modern ratio (Evans et al., 2016a);

• v.

  included the effect of pH or [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] on δ¹⁸O (Spero et al., 1997; Zeebe, 1999) given the covariance of temperature and pH described in point (iii) above using the multispecies average slope of Gaskell and Hull (2023); and, finally,

• vi.

  explored the likely impact of the covariance of $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ and temperature that is characteristic of the modern ocean and arises from the broad coupling of the hydrological cycle with surface temperatures.

This latter influence was calculated by combining sea surface temperature (SST) data from the 2013 World Ocean Atlas (Locarnini et al., 2013) and $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ from LeGrande and Schmidt (2006), taking all surface ocean data except those from polar meltwater regions, which demonstrates that, on average, in the modern ocean $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ increases by 0.0425 ‰ per degree (Celsius) of SST increase.

Figure 2 Raw δ¹⁸O plotted against $\mathrm{Mg}/\mathrm{Ca}$ for all samples presented here. (a) Several possible expected $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O slopes are shown for comparison, including the slope for G. ruber and T. sacculifer (at constant pH) in the modern ocean (solid and dashed black lines respectively). The additive impact of nonthermal controls is then explored using the G. ruber calibration as an example, specifically (a) the impact of a whole-ocean pH shift of 0.15 units (orange line), accounting for the covariation of pH and temperature (driven by the temperature-dependent dissociation of water, red line), seawater $\mathrm{Mg}/\mathrm{Ca}$ half of the present day value (thin blue line), the theoretical impact of pH on δ¹⁸O (blue line), and the covariance of temperature and $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ in the modern ocean (thick blue line). We stress that the magnitude of the pH, $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$, and $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ changes applied to these theoretical lines was chosen to illustrate the sensitivity of the two proxies to these factors and does not necessarily represent the degree to which these factors differed in any specific time interval compared to today. The length of each line depicts the expected $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O change across the same temperature range in each case (5–35 °C). All calculations assume $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}=\mathrm{0}\mathrm{‰}$. (b) As in panel (a), except with the long-term whole-ocean changes in pH, $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$, and $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ subtracted from the raw proxy values (see text, using the multispecies calibration of Gray and Evans, 2019, in the case of the $\mathrm{Mg}/\mathrm{Ca}$ corrections), i.e. accounting for the impact of these non-thermal $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O controls. Sample age is shown as a function of colour. Note that the black line is not a regression but rather shows one possible theoretical $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O relationship and is identical to the thick blue line from panel (a) including the combined effect of all factors explored there. The slope of this theoretical relationship is −2.08, which provides an approximate expected slope of Miocene–recent $\mathrm{Mg}/\mathrm{Ca}$ versus δ¹⁸O, although we note that no real-world dataset will conform exactly to this as pH, $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$, and $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ have not (perfectly) covaried through time.

Each of these factors was applied additively such that factor (iv) listed above ($\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$) in Fig. 2 includes factors (i) through (iii). The sum of the influence of these factors on the theoretical δ¹⁸O–$\mathrm{Mg}/\mathrm{Ca}$ relationship is represented by the thick blue line in Fig. 2a and the black line in Figs. 2b, 3, 5, 6, and 7, which has a slope of −2.08 in δ¹⁸O–$\ln \left(\mathrm{Mg}/\mathrm{Ca}\right)$ space.

The magnitude of some of these potential non-thermal controls on the two proxies over the time interval studied here is reasonably well constrained. Specifically, the long-term whole-ocean pH and $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$ changes are sufficiently well known (Rae et al., 2021; Zhou et al., 2021; Brennan et al., 2013) that they can be “subtracted” out of the raw proxy values, given that they are likely to apply to all or most species in the dataset. As such, we next explored the degree to which the $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O covariation improves once long-term whole-ocean pH and $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$ changes are removed. To avoid (possibly incorrect) a priori assumptions regarding, for example, which $\mathrm{Mg}/\mathrm{Ca}$–temperature calibration should be applied to each species in the dataset and the degree to which surface ocean $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ has varied at the study sites, we initially kept the $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O comparison in raw proxy space and did the following: (1) converted the raw $\mathrm{Mg}/\mathrm{Ca}$ values to temperature using the multispecies $\mathrm{Mg}/\mathrm{Ca}$–temperature calibration from Gray and Evans (2019), together with our best estimate of pH and $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$ (as described below (Sect. 2.4.2)), and (2) converted the temperatures back into $\mathrm{Mg}/\mathrm{Ca}$ using the same calibration but modern seawater $\mathrm{Mg}/\mathrm{Ca}$ and pH. Conceptually, this is equivalent to correcting the raw proxy values for these nonthermal controls. In addition, we subtracted out the long-term whole-ocean change in $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ related to continental ice growth using the sea level curve of Rohling et al. (2022) and a sea level–$\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ scaling factor of 1 ‰ per 67 m. This results in a raw proxy dataset in which the aforementioned long-term non-thermal factors are no longer present and which can be used to evaluate the occurrence of residual scatter independent of the long-term non-thermal controls on $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O (Fig. 2b).

2.4.2 Transformation of proxy values into paleotemperature

Measured foraminifera $\mathrm{Mg}/\mathrm{Ca}$ was transformed into paleotemperature using the MgCaRB tool (Gray and Evans, 2019; https://github.com/dbjevans/MgCaRB/releases/tag/v1.3, last access: 11 February 2025 (MATLAB); https://doi.org/10.5281/zenodo.14056982) which takes the following into account:

• Salinity. Although this has a minor effect on $\mathrm{Mg}/\mathrm{Ca}$ (Hönisch et al., 2013), whole-ocean changes are nonetheless accounted for using a salinity reconstruction derived from scaling the δ¹⁸O benthic stack (Westerhold et al., 2020) to the sea level record of Spratt and Lisiecki (2016) between 0.8 and 8 Ma, of Rohling et al. (2022) between 8 and 40 Ma, and assuming an ice-free world before this time. We applied the multispecies salinity sensitivity of Gray and Evans (2019) to all species (3.6 % per salinity unit).

• pH. Long-term whole-ocean changes were derived from a smoothing spline fit to the boron-isotope-derived pH data compiled by Rae et al. (2021). We applied species-specific pH–$\mathrm{Mg}/\mathrm{Ca}$ sensitivities of Gray and Evans (2019) where available for a given species/lineage (discussed in more detail below) and used the multispecies sensitivity in all other cases.

• $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$. This was derived by combining the [${\mathrm{Ca}}_{\text{sw}}^{\mathrm{2}+}$] record of Zhou et al. (2021) with a smoothing spline fit to the fluid inclusion [${\mathrm{Mg}}_{\text{sw}}^{\mathrm{2}+}$] data given in Brennan et al. (2013). Raw $\mathrm{Mg}/\mathrm{Ca}$ values were adjusted using the equation ${\mathrm{Mg}/\mathrm{Ca}}_{\text{corrected}}={\mathrm{Mg}/\mathrm{Ca}}_{\text{raw}}\times {\mathrm{Mg}/\mathrm{Ca}}_{\text{sw}}^H/{\mathrm{5.2}}^H$, where H=0.64 based on a data compilation of three foraminifera species and inorganic calcite (Holland et al., 2020; Evans et al., 2015, 2016b; Mucci and Morse, 1983).

We initially applied a multispecies equation because the dataset includes a mix of extant and extinct species, some of which have never been measured for trace elements before (Table S2) or lack an extant/well-calibrated modern relative. Specifically, we used the multispecies $\mathrm{Mg}/\mathrm{Ca}$–temperature equation of Gray and Evans (2019) and applied the multispecies pH, salinity, and temperature sensitivities, together with the Globigerinoides ruber pre-exponential coefficient (given a $\mathrm{Mg}/\mathrm{Ca}$–T relationship of the form $\mathrm{Mg}/\mathrm{Ca}=\text{BeAT}$ in its simplest form) as most of the species for which high-quality data exist are known to be characterized by a $\mathrm{Mg}/\mathrm{Ca}$–pH sensitivity (Lea et al., 1999; Kisakürek et al., 2008; Evans et al., 2016b).

We subsequently applied species-specific calibrations to selected lineages to explore the degree to which scatter in the dataset can be accounted for by taking into account phylogenetic relationships among ancestor-descendent species. Specifically, the Trilobatus sacculifer calibration was applied to the Trilobatus trilobus–Trilobatus sacculifer lineage, and the Orbulina universa calibration was applied to the Praeorbulina–Orbulina lineage, both from laboratory culture studies following Gray and Evans (2019). To Neogloboquadrina and its descendent lineage Pulleniatina we applied a Neogloboquadrina pachyderma calibration with the sensitivities (and their uncertainties) of Tierney et al. (2019) (implemented with a re-fit to the dataset following the MgCaRB approach; see Fig. S1). We then evaluate the improvement relative to the multispecies calibration in samples spanning the middle Miocene to modern. The attribution of phylogenetic relationships follows Aze et al. (2011), Spezzaferri et al. (2018), Leckie et al. (2018), and Fabbrini et al. (2021). All uncertainties were fully propagated via Monte Carlo simulation, including those related to analysis; calibration coefficients (the coefficients that relate $\mathrm{Mg}/\mathrm{Ca}$ to temperature, pH, and salinity); and the 95 % confidence intervals on the salinity, pH, and $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$ reconstructions. Specifically, 10⁴ random draws of each of these values within their respective uncertainty bounds were used to generate the reported values and 95 % CI (percentiles 2.5, 50, and 97.5 of the resulting dataset). The relative contributions of the various nonthermal factors impacting $\mathrm{Mg}/\mathrm{Ca}$-derived temperature are shown in Fig. S2.

The conversion of δ¹⁸O to paleotemperature followed Gaskell and Hull (2023) using the bayfox calibration (Malevich et al., 2019) and the global and local $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ of Rohling et al. (2022) and Gaskell and Hull (2023) respectively. The calculation was performed twice, both with and without a $\mathrm{pH}/\left[{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}\right]$ effect on δ¹⁸O (the former using the mean planktonic foraminiferal slope of Gaskell and Hull (2023), and the [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] record of Zeebe and Tyrrell, 2019). We also explored the impact of using a different sea level/deep-ocean temperature reconstruction (Miller et al., 2020) on our results, described in the Results Section and shown in Fig. S3.

When evaluating the paleotemperature reconstructions, we define whether or not the two proxy systems agree within uncertainty by determining if the root sum of squares of the two uncertainties is smaller than the temperature difference between the two proxies. We then proceed to identify possible drivers for the data deviating from the expected $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O relationship by evaluating the age of the sample, regional changes in $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$, depth ecology, and possible species-specific offsets.

We note that all of the above corrections assume surface ocean conditions, while the dataset contains a number of species that calcify at depth (Boscolo-Galazzo et al., 2021). Given the uncertainties surrounding past changes in vertical pH and $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ profiles, we do not attempt to account for this in our data analysis but note that this consideration should be borne in mind when interpreting data from deep-dwelling species.

3 Results

Our basic expectation is that higher $\mathrm{Mg}/\mathrm{Ca}$ should relate to more negative δ¹⁸O values for warmer temperatures and vice versa for colder temperatures. Despite the large number of variables included, the dataset as a whole shows a positive correlation (Fig. 2a; R²=0.37, RMSE=1.01, p≪0.01) between δ¹⁸O and $\ln \left(\mathrm{Mg}/\mathrm{Ca}\right)$. Hence, the δ¹⁸O–$\mathrm{Mg}/\mathrm{Ca}$ covariance can be considered a robust feature over the past 15 Myr for the majority of the species analysed and across the study sites. Nonetheless, there is a high degree of scatter in these data, which suggests that the temperature signal which should lead $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O data to change consistently in opposite directions is affected by other factors. Our exercise of generating theoretical $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O relationships (Fig. 2), exploring how the relationship between the two proxies might change through space and time, provides a qualitative indication as to whether the scatter can be attributed to long-term non-thermal factors generally corrected for when using the δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ proxies. The substantial differences between these expected relationships suggests that this is likely to be the case (Fig. 2a). In particular, Fig. 2a suggests that both a pH effect and $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$ changes through time may explain a substantial degree of the variability observed in the dataset compared to the modern relationships (compare the coloured and black lines; see also Fig. S2). Nonetheless, accounting for these long-term biases alone in the raw dataset does not remove all of the scatter (Fig. 2b), suggesting the importance of additional factors, such as vital effects and regional variations in $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$. Data converted into temperature, along with 95 % confidence intervals, are shown in Fig. 3. In this plot 62 % of data points fall within uncertainty, confirming that a high degree of variability in the raw data can be effectively explained and accounted for by correcting for the known spatially and temporally varying non-thermal effects influencing both proxies.

Figure 3 δ¹⁸O- versus $\mathrm{Mg}/\mathrm{Ca}$-derived paleotemperatures plotted as a function of age (a) and calcification depth (b), accounting for the impact of whole-ocean and regional changes in $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ following Gaskell et al. (2022) using the bayfox δ¹⁸O–temperature calibration (Malevich et al., 2019), whole-ocean changes in $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$ and pH on $\mathrm{Mg}/\mathrm{Ca}$ using MgCaRB (Gray and Evans, 2019), and including a pH correction on δ¹⁸O using the mean planktonic foraminifera slope (Gaskell and Hull, 2023). The grey envelopes in (a) and (b) are fully propagated uncertainties in both proxies incorporating analysis, calibration, pH, $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$, and $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$/salinity, while the dotted black lines are identical to the thick blue and black lines in Fig. 2a and b (see text for details). Site 516 data are shown with square symbols in panel (a).

Using the deep-ocean temperature and sea level records of Miller et al. (2020) rather than Rohling et al. (2022) (Fig. S3) would shift our δ¹⁸O temperature reconstructions of +0.5 to 3.3 °C depending on sample age, with the samples from the 7.5 Ma time slice being most affected. Whilst the choice to correct for global $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ exerts a significant control on δ¹⁸O-derived temperatures, it has a minor impact on both our long-term dataset overall (compare Fig. 3 with Fig. S3) and the proportion of proxy data within combined uncertainty, from 62 % to 64 %.

We find that deep-dwelling and surface-dwelling species fall within uncertainty in terms of the broad degree of agreement between $\mathrm{Mg}/\mathrm{Ca}$- and δ¹⁸O-derived temperatures; hence different depth-habitat ecologies or calcification environments do not represent a systematic source of offset per se (Fig. 3b). Indeed, all the species displaying a temperature offset between the two proxies are surface dwellers (0–200 m depth, Boscolo-Galazzo et al., 2022), with deep dwellers characterized by inter-proxy agreement in almost all cases (Fig. 3b). Additional explanations are required for those temperature reconstructions that, all non-thermal and spatial controls included, still plot outside their propagated error uncertainty. Data points outside the combined proxy uncertainties either display >5 °C warmer (colder) δ¹⁸O ($\mathrm{Mg}/\mathrm{Ca}$) temperatures (upper-left area of the plot) or >5 °C warmer (colder) $\mathrm{Mg}/\mathrm{Ca}$ (δ¹⁸O) temperatures (bottom-right area of the plot) (Fig. 3).

Figure 3a clearly shows that the large majority of the outliers in the upper-left part of the plot (warmer δ¹⁸O temperatures or cooler $\mathrm{Mg}/\mathrm{Ca}$ temperatures) are species of late Miocene to modern age, while the outliers in the bottom-right part of the plot (warmer $\mathrm{Mg}/\mathrm{Ca}$ temperatures) are mostly older species of middle Miocene age or from mid-latitude Site 516 (squares in Fig. 3a). We suggest and discuss possible main sources for these offsets between $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O temperatures in detail below.

4 Discussion

4.1 Diagenesis

Diagenesis is known to alter the test chemistry of foraminifera in three main ways: partial dissolution, overgrowth, and recrystallization (Edgar et al., 2015). The trace element and isotopic composition of tests react differently to these diagenetic processes. The trace element composition of foraminiferal calcite may be susceptible to partial dissolution because it is inhomogeneous (Fehrenbacher and Martin, 2014), which decreases trace element ratios in species with high- and low-Mg regions (Dekens et al., 2002; Edgar et al., 2015; Rongstad et al., 2017), resulting in lower temperature reconstructions. Overgrowth and recrystallization have been shown to add both low-Mg and high-Mg diagenetic calcite, potentially impacting the original signal in opposite directions (Branson et al., 2015), although $\mathrm{Mg}/\mathrm{Ca}$ is relatively robust to this type of diagenesis, at least in certain circumstances (Staudigel et al., 2022). The oxygen isotopic composition of planktonic foraminiferal tests is well known to be very sensitive to overgrowth and recrystallization (e.g. Sexton et al., 2006), whereby the addition of diagenetic calcite, or the replacement of the original calcite with diagenetic calcite precipitated at the seafloor, can significantly alter the original isotopic signal, shifting it to more positive values (Pearson, 2012; Edgar et al., 2015).

The Boscolo-Galazzo et al. (2021) dataset spans 15 Myr and includes sites with different average preservation of foraminiferal tests and oceanographic settings. When these data are regressed against each other (Fig. 2b), we find a total of 26 data points characterized by oxygen isotope values more positive than expected from their $\mathrm{Mg}/\mathrm{Ca}$ values (Fig. 2b; Table S1), resulting in δ¹⁸O temperatures >5 °C colder than $\mathrm{Mg}/\mathrm{Ca}$ temperatures (outside the error envelope) (Fig. 3; Table S1). Of these data points, 21 were from the older time slices (12.5 Ma and 15 Ma), 1 from the 7.5 Ma time slice and 4 from the core top of Site U1338 (Table S1).

As the majority of these data points were characterized by $\mathrm{Mg}/\mathrm{Ca}$ values of $\sim \mathrm{1.5}\text{–}\mathrm{2}\ln \left(\mathrm{Mg}/\mathrm{Ca}\right)$ (Fig. 2b), thus yielding more reasonable $\mathrm{Mg}/\mathrm{Ca}$ than δ¹⁸O temperatures for these sites/time intervals (Table S1), the observed offset is most likely best attributed to diagenetic overgrowth/recrystallization, shifting oxygen isotopes towards more positive values without affecting $\mathrm{Mg}/\mathrm{Ca}$ to the same extent. A recent study compared typical $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O from recrystallized planktonic foraminifera with chemical diffusive models simulating early diagenetic processes in calcite (Staudigel et al., 2022). According to that study, in a closed system, the bulk δ¹⁸O value will be altered faster than the $\mathrm{Mg}/\mathrm{Ca}$, regardless of what partitioning coefficient is used for Mg, leading to a progressive shift to more positive δ¹⁸O values leaving $\mathrm{Mg}/\mathrm{Ca}$ virtually unchanged (Staudigel et al., 2022).

Data points characterized by δ¹⁸O overprinted by overgrowth/recrystallization were distributed across most of the study sites (except for Sites 871/872), but they were more common at Site U1490 and U1489 ($\mathrm{15}/\mathrm{26}$) which are characterized by inferior preservation compared to the others (Boscolo-Galazzo et al., 2021). The core top samples at Site U1338 show clear signs of dissolution with highly fragmented tests. These data points presented the lowest $\mathrm{Mg}/\mathrm{Ca}$ values in the dataset (1.68 and 0.91 mmol mol⁻¹), with temperatures from $\mathrm{Mg}/\mathrm{Ca}$ lower than from δ¹⁸O. This suggests that partial dissolution and recrystallization affected both $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O in this sample but $\mathrm{Mg}/\mathrm{Ca}$ more so.

Overall, our scrutiny for diagenesis of the Boscolo-Galazzo et al. (2021) dataset is consistent with previous studies suggesting that δ¹⁸O values are more easily affected by recrystallization than $\mathrm{Mg}/\mathrm{Ca}$ (Sexton et al., 2006; Staudigel et al., 2022; John et al., 2023), similar to other trace element systems (Edgar et al., 2015).

Based on the considerations above we excluded the affected 26 data points from the subsequent analysis as being characterized by a stronger diagenetic overprint than the rest of the dataset (Table S1). Removing the affected data points in some but not all the cases equated to removing a whole sample (Table S1). This is because of variable diagenetically offset δ¹⁸O values from different species in a sample, as observed elsewhere (Sexton et al., 2006; Edgar et al., 2015). The approach used here to reconstruct δ¹⁸O temperature shows that for the rest of the study dataset a diagenetic offset would be comprised within the propagated error envelope (∼2–2.5 °C) (Fig. 3) and comparable or not distinguishable from an offset deriving from poorly constrained $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$.

4.2 Regional-scale spatial heterogeneity in seawater chemistry

Once converted into temperatures the dataset shows an overall good $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O agreement, consistent with temperature being a dominant controller of both proxies through time and across the broad geographical area investigated (Fig. 3). This suggests that, by and large, the seawater corrections applied for the local and global changes in ocean chemistry are adequate, although for one site this may not hold true. Site 516 is a mid-latitude site (south-west Atlantic, 30° S) characterized by a modern sea surface temperature around 20 °C (Fig. 1); most of these data points from this site have very positive δ¹⁸O values associated with high $\mathrm{Mg}/\mathrm{Ca}$ (Fig. 2b; Table S1). Once converted, this results in δ¹⁸O temperatures that are too cold (12–17 °C) compared to both modern (given long-term warming since the Miocene is not expected at any of these sites) and the equivalent $\mathrm{Mg}/\mathrm{Ca}$ temperatures from the same samples, which are around 21–25 °C (Fig. 3a; Table S1). We do not attribute this mismatch to diagenesis for a number of reasons. First, Site 516 is characterized by a very good test preservation, much better than at Site U1489 and U1490 for which diagenesis extensively affects the middle Miocene samples; second, this mismatch is observed in the entire dataset through samples spanning the middle Miocene to modern; and third, the mismatch is observed for surface-dwelling species only, with deep dwellers characterized by δ¹⁸O–$\mathrm{Mg}/\mathrm{Ca}$ in good agreement. Site 516 is situated in an area of complex surface hydrography, as it sits at the confluence between the warm Brazil western boundary current and the cold Falkland (Malvinas) Current spinning off from the Antarctic Circumpolar Current (e.g. Jonkers et al., 2021). Compared to the subtropical gyres, where many of the study samples come from, a large degree of spatial variability of surface water physical–chemical properties can be expected on a seasonal and multiannual scale. As such, we suggest that the mismatch between δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ observed in surface-dwelling species at Site 516 may result from changeable surface water properties from the mixing of two very different water masses creating deviations in pH, salinity, and $\mathit{\delta }{{}^{\mathrm{18}}\mathrm{O}}_{\text{sw}}$ beyond those that are typical in stratified open ocean environments and that are difficult to correct for. Sites with a changeable hydrography, such as Site 516, may hence not be ideal for the application of geochemical proxies affected by seawater chemistry, unless changes in seawater chemistry at the site can be reconstructed directly.

4.3 Species-specific offsets

The third and largest source of mismatch that we consider is the occurrence of species-specific offsets, particularly in $\mathrm{Mg}/\mathrm{Ca}$ given that, in general, the relative degree of inter-species $\mathrm{Mg}/\mathrm{Ca}$ variability is greater than for oxygen isotopes. For example, the range of inter-species δ¹⁸O offsets is around 1 ‰ (<4 °C), whereas there may be more than a factor of 2 difference in $\mathrm{Mg}/\mathrm{Ca}$ between species for a given temperature (e.g. compare Pearson, 2012; Gray and Evans, 2019; Regenberg et al., 2009; Däeron and Gray, 2023). When regressing these data against each other, we find several species characterized by systematically offset $\mathrm{Mg}/\mathrm{Ca}$ values (e.g. Fig. 2b; Table S1) leading to proxy–proxy disagreement (Fig. 3; Table S1). The occurrence of such offset $\mathrm{Mg}/\mathrm{Ca}$ values is not evenly distributed across species but is shared among related species in both spinose and non-spinose groups (Figs. 4 and 5; Table S1). Specifically, the spinose species Orbulina universa (2 out of 5 specimens, $\mathrm{2}/\mathrm{5}$), O. suturalis ($\mathrm{3}/\mathrm{6}$), and Praeorbulina glomerosa ($\mathrm{1}/\mathrm{1}$) present high $\mathrm{Mg}/\mathrm{Ca}$ ratios compared to their δ¹⁸O values and $\mathrm{Mg}/\mathrm{Ca}$ of other species (Fig. 5e). We also find that Globigerinella siphonifera ($\mathrm{6}/\mathrm{7}$), G. calida ($\mathrm{1}/\mathrm{1}$), G. praesiphonifera ($\mathrm{3}/\mathrm{4}$), and Globigerina bulloides ($\mathrm{5}/\mathrm{6}$) have offset $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O values, largely being characterized by higher-than-expected $\mathrm{Mg}/\mathrm{Ca}$, although three G. siphonifera data points show lower $\mathrm{Mg}/\mathrm{Ca}$ (Figs. 4 and 5k; Table S1). Among non-spinose species, Neogloboquadrina humerosa (11 out of 12 specimens, $\mathrm{11}/\mathrm{12}$), N. acostaensis ($\mathrm{3}/\mathrm{3}$), Pulleniatina obliquiloculata ($\mathrm{5}/\mathrm{6}$), P. praecursor ($\mathrm{3}/\mathrm{5}$), P. primalis ($\mathrm{4}/\mathrm{6}$), Sphaeroidinella dehiscens ($\mathrm{4}/\mathrm{4}$), and Sphaeroidinellopsis paenedehiscens ($\mathrm{5}/\mathrm{10}$) display $\mathrm{Mg}/\mathrm{Ca}$ values lower than expected for their oxygen isotope composition (Figs. 4 and 5g, i; Table S1). For the first time, these results highlight the occurrence of similarly offset $\mathrm{Mg}/\mathrm{Ca}$ values for ancestor-descendent species belonging to the same lineage as well as for whole new lineages (Figs. 4 and 5).

Figure 4 Phylogenetic relationships of offset spinose and non-spinose species. Shown here are the species discussed in the text and their most closely related species and ancestors. Coloured lines indicate species offset in $\mathrm{Mg}/\mathrm{Ca}$ in the study dataset relative to a multispecies calibration approach, and black lines indicate non-offset species. Red lines indicate offset species with higher $\mathrm{Mg}/\mathrm{Ca}$, blue lines offset species with lower $\mathrm{Mg}/\mathrm{Ca}$, and purple lines species displaying both types of offsets. Question marks indicate the lack of $\mathrm{Mg}/\mathrm{Ca}$ data for a given species in the dataset presented here. Phylogeny after Aze et al. (2011), Spezzaferri et al. (2018), Leckie et al. (2018), and Fabbrini et al. (2021). The phylogenetic chart was generated using Mikrotax (Huber et al., 2016; http://www.mikrotax.org/pforams, last access: 4 June 2024). Ages are from Lourens et al. (2004) for the Neogene and Pälike et al. (2006) for the Oligocene.

Figure 5 δ¹⁸O versus $\mathrm{Mg}/\mathrm{Ca}$ and proxy-derived paleotemperature estimates for the ancestor-descendent species Globigerinoides subquadratus–G. ruber (a, b); Trilobatus trilobus–T. sacculifer (c, d); Praeorbulina glomerosa–Orbulina suturalis–O. universa (e, f); Neogloboquadrina acostaensis–N. humerosa, N. pachyderma–N. incompta, and N. acostaensis–Pulleniatina primalis–P. praecursor–P. finalis (g, h); Sphaeroidinellopsis seminulina–S. paenedehiscens–Sphaeroidinella dehiscens and S. seminulina–S. kochi (i, j); and Globigerinella praesiphonifera–G. siphonifera–G. calida and Globigerina bulloides (k, l). Species are defined as offset relative to a multispecies calibration approach when their reconstructed temperature plots outside the grey error envelope (see text for details). Circles indicate data points from Site 516. Raw proxy values are given with the long-term non-thermal controls on $\mathrm{Mg}/\mathrm{Ca}$ subtracted out (as in Fig. 2), as well as an estimate of paleotemperature (as in Fig. 3). The continuous and dotted black lines depict one possible estimate of the expected slope between δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ (the thick blue and black lines from Fig. 2), adjusted to approximately match the location of these data by shifting them in the direction of δ¹⁸O. Data points which are considered strongly affected by diagenesis are not included in this plot. Note that one data point in panel (g) falls outside of the plot area.

Divergent $\mathrm{Mg}/\mathrm{Ca}$ values for G. siphonifera and O. universa have previously been reported (Opdyke and Pearson, 1995; Anand et al., 2003; Friedrich et al., 2012). In the case of O. universa, the offset may be related to pH change in the foraminiferal microenvironment due to symbiont photosynthetic activity (Eggins et al., 2004) or changes in seawater pH, with $\mathrm{Mg}/\mathrm{Ca}$ of the test increasing by as much as 6 %±3 % for each 0.1 unit decrease in pH (Lea et al., 1999; Russell et al., 2004). pH-related vital effects are reported for other spinose species of planktonic foraminifera such as Globigerina bulloides (Lea et al., 1999; Davis et al., 2017), which is related to the genus Globigerinella (Fig. 4).

Among the neogloboquadrinids, N. acostaensis and its descendant N. humerosa have the most clearly expressed offset with low $\mathrm{Mg}/\mathrm{Ca}$ values in the study dataset. In contrast, Neogloboquadrina pachyderma and N. incompta are not offset in the study dataset, perhaps simply because of the limited number of data (one data point each). More broadly, a $\mathrm{Mg}/\mathrm{Ca}$ offset compared to other species has been reported in the literature (Davis et al., 2017). Neogloboquadrina dutertrei, N. incompta, N. pachyderma, and Pulleniatina obliquiloculata have been shown to be characterized by much lower trace element concentrations (Mg–Ba–$\mathrm{Zn}/\mathrm{Ca}$) in the adult portions of their shells (crust and cortex), so that a greater amount of adult versus early ontogenetic calcite leads to low trace element values in bulk shell analysis (Jonkers et al., 2012; Davis et al., 2017; Fritz-Endres and Fehrenbacher, 2021). The low $\mathrm{Mg}/\mathrm{Ca}$ of crust and cortex has been found to be independent of ambient temperature in cultured Neogloboquadrina (Davis et al., 2017). It is found in specimens collected both in surface waters and at depth (Jonkers et al., 2021), indicating that the low $\mathrm{Mg}/\mathrm{Ca}$ is not acquired due to calcification in deeper, colder waters of the crust/cortex portion of the shell, although a greater incidence of crusts is reported for colder waters (Jonkers et al., 2021). Neogloboquadrina, Pulleniatina, and Sphaeroidinella/Sphaeroidinellopsis are all characterized by a thick crust or cortex, suggesting their $\mathrm{Mg}/\mathrm{Ca}$ is biased by low Mg adult calcite being quantitatively predominant, which is further corroborated by their $\mathrm{Mg}/\mathrm{Ca}$ being unrelated to temperature in our data and consistently falling outside of the δ¹⁸O-derived temperatures, even accounting for the combined uncertainty of the two proxies (Fig. 5g–h and i–j).

The majority of the data points from the offset spinose and non-spinose species result in temperature differences between the two proxies greater than 5 °C when using the multi-species calibration from Gray and Evans (2019) as described in Sect. 2.4.2, hence outside the calculated error envelope when taking all the non-thermal factors discussed above into account (Fig. 5). A similar temperature offset is not apparent for other lineages such Trilobatus trilobus–Trilobatus sacculifer and Globigerinoides subquadratus–G. ruber, to which the same treatment to the offset spinose and non-spinose species was applied (Sect. 2.4.2, Fig. 5a–d). Hence, we attribute the offset temperatures to the atypical $\mathrm{Mg}/\mathrm{Ca}$ signatures described above in the affected species, in turn resulting from biology-/ecology-dependent vital effects shared within a lineage and between related lineages (Figs. 4 and 5).

When a species-specific calibration for Neogloboquadrina pachyderma is applied to descendent species/lineages and sister taxa (Fig. 6) the offset is successfully corrected for all the Neogloboquadrina and Pulleniatina species which effectively no longer produce offset temperatures (Fig. 6). In contrast, we only observe a minor improvement when applying the Orbulina universa calibration to the Praeorbulina–Orbulina lineage, with most temperature data points remaining offset (Fig. 6). This may imply that the O. universa laboratory calibrations require revision for application to fossil samples. No large difference is observed when applying the Trilobatus sacculifer calibration to ancestor-descendent species in the genus Trilobatus (Figs. 5–7) although we recommend doing so, given that no $\mathrm{Mg}/\mathrm{Ca}$–pH effect is known for this genus, in contrast to (e.g.) G. ruber.

Figure 6 As in Fig. 5 extending the use of species-specific calibrations to all species in a lineage in the case of the Trilobatus trilobus–T. sacculifer (c, d), Praeorbulina glomerosa–Orbulina suturalis–O. universa (e, f), Neogloboquadrina acostaensis–N. humerosa, N. pachyderma–N. incompta, N. acostaensis and between related lineages in the case of Neogloboquadrina and the Pulleniatina primalis–P. praecursor–P. finalis lineage (g, h). When applying species-specific calibrations, the offset between $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O temperature is resolved for some species (e.g. all Neogloboquadrina, compared to Fig. 5) but not for others; i.e. their reconstructed temperature still plots outside the grey error envelope (e.g. most of Globigerina and Globigerinella).

Figure 7 As Fig. 3, with all offset lineages (Sect. 4.3; specifically, those that remain offset following the application of lineage-specific calibrations where possible) and diagenetically compromised (Sect. 4.1) samples removed. Removing these samples leaves 157 data points, of which 92 % fall within the combined uncertainty of $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O agreement (grey envelope).

Overall, this exercise demonstrates that the majority of data points characterized by proxy–proxy disagreement (Fig. 3) are from the lineages Praeorbulina–Orbulina, Globigerina–Globigerinella, Neogloboquadrina–Pulleniatina, and Sphaeroidinellopsis–Sphaeroidinella (Fig. 5). We find that using a “nearest descendant” approach in the choice of temperature calibration improves the agreement between δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ temperatures for the neogloboquadrinids and pulleniatinids (Fig. 6). At the same time, it enables us to identify “problematic” species and lineages which require further investigations before being used for temperature reconstructions (Fig. 6). Once all the potential sources of offset described above are taken into account, the species-specific calibration for Neogloboquadrina is applied, and data points which are still offset are removed from the dataset, the agreement between the two proxies increased from 62 % to 91 % of data points falling within the combined uncertainties of the proxies (Fig. 7).

4.4 Planktonic foraminiferal $\mathrm{Mg}/\mathrm{Ca}$ offsets as an expression of evolution

The analysis of the Boscolo-Galazzo et al. (2021) dataset performed here allows offsets through ancestor-descendent species to be tracked for the first time and the time of their appearance to be assessed. Specifically, we observe that offsets tend to appear with the origination of a new lineage and to continue to the modern representative(s) of that lineage (Fig. 8). In this way, an attempt can be made to interpret offsets as the geochemical expression of evolutionary new biochemical pathways or ecological strategies in emerging species.

Figure 8 Simplified phylogeny for the offset lineages discussed in the text with the occurrence of morphological changes associated with their evolution highlighted. Coloured lines indicate species offset in $\mathrm{Mg}/\mathrm{Ca}$ in the study dataset relative to a multispecies calibration approach, and black lines indicate non-offset species. Red lines indicate offset species with higher $\mathrm{Mg}/\mathrm{Ca}$, blue lines indicate offset species with lower $\mathrm{Mg}/\mathrm{Ca}$, and purple lines denote species displaying both types of offsets. Question marks indicate the lack of $\mathrm{Mg}/\mathrm{Ca}$ data for a given species in the dataset presented here. Phylogeny after Aze et al. (2011) and Spezzaferri et al. (2018). The phylogenetic chart was generated using Mikrotax (Huber et al., 2016; http://www.mikrotax.org/pforams, last access: 4 June 2024). Ages are from Lourens et al. (2004) for the Neogene and Pälike et al. (2006) for the Oligocene.

For spinose species, the observed high $\mathrm{Mg}/\mathrm{Ca}$ offset is shared by ancestor-descendent species such as Globigerinella praesiphonifera–G. siphonifera and Praeorbulina glomerosa–Orbulina suturalis–O. universa (Fig. 4). Globigerina bulloides shares the same type of offset as Globigerinella and a common ancestor (Globigerina archaeobulloides) in the earliest Oligocene (∼33.5 Ma) (Spezzaferri et al., 2018) (Fig. 4), suggesting that for this group the offset may go back to at least the early globigerinids of the Paleogene. The genus Praeorbulina originated from the genus Trilobatus at about 16 Ma (Fig. 4) (Pearson et al., 1997; Aze et al., 2011). Trilobatus trilobus is the last common ancestor between the Trilobatus and Praeorbulina–Orbulina lineages (Fig. 4) and does not present offset $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O values, similar to its modern descendants (Figs. 4 and 5). This suggests that the offset in Praeorbulina–Orbulina originated within the lineage and the morphological changes associated with it (Fig. 8) and carried on to the modern representative O. universa. Spinose Globigerinoides ruber is also reported to be sensitive to pH changes (Kisakürek et al., 2008; Evans et al., 2016a), G. ruber is not offset in the analysed dataset (both in the raw data and calculated temperatures), with ancestor-descendent G. subquadratus–G. ruber behaving similarly to the Trilobatus trilobus–T. sacculifer lineage through time (Figs. 5 and 6), suggesting that this non-thermal effect is adequately accounted for in this case (Fig. 6) (i.e. the laboratory calibrations are applicable downcore into deep time in correcting for this). The $\mathrm{Mg}/\mathrm{Ca}$ displayed by offset spinose Praeorbulina–Orbulina, Globigerinella species, and G. bulloides in the study dataset that is higher on average may suggest a lower pH environment, which we cannot directly account for (Fig. 2a). For Globigerina and Globigerinella, where a larger degree of scatter (Figs. 5 and 6k, l) is observed in association with a rather conservative morphology through time (Fig. 8), the offset may be linked to an opportunistic behaviour and capability to adapt to a broad range of environmental conditions with variable pH (Weiner et al., 2015). This may in turn be related to the complexity of genotypes association in both G. bulloides and Globigerinella species, with different genotypes having different ecologies but almost identical morphologies (e.g. Weiner et al., 2015; Morard et al., 2024).

We find that the low $\mathrm{Mg}/\mathrm{Ca}$ offset is shared between ancestor-descendent lineages Neogloboquadrina–Pulleniatina and Sphaeroidinellopsis–Sphaerodinella (Figs. 4 and 8). Neogloboquadrina evolved from Paragloborotalia continuosa in the late Miocene, at about 10 Ma (Fig. 4). While paired $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O measurements for P. continuosa are not available, paired $\mathrm{Mg}/\mathrm{Ca}$–δ¹⁸O measurements for Paragloborotalia siakensis, a species older than P. continuosa, do not show a low $\mathrm{Mg}/\mathrm{Ca}$ offset (Table S1). This may indicate that the occurrence of low $\mathrm{Mg}/\mathrm{Ca}$ crust/cortex started either with P. continuosa, the youngest representative of the genus Paragloborotalia in our study, or with the neogloboquadrinids. The genus Pulleniatina evolved from Neogloboquadrina acostaensis at about 6.5 Ma (Pearson et al., 2023; see Fig. 4), by modifying the chambers' arrangement and progressively developing a cortex (Fig. 8). Pulleniatina may have inherited the capability to thicken its test from Neogloboquadrina and modified it into a cortex (Fig. 8) (Pearson et al., 2023). The occurrence of a low $\mathrm{Mg}/\mathrm{Ca}$ cortex in Pulleniatina appears to start from the most ancient representative of this group, P. primalis (Fig. 4) and continues to the modern. Similar to Paragloborotalia, middle Miocene Sphaeroidinellopsis kochi and S. seminulina are not characterized by an offset to low $\mathrm{Mg}/\mathrm{Ca}$ values (Figs. 4–6). Nonetheless, an offset is observed in 5 out of 10 specimens for late Miocene–early Pliocene S. paenedehiscens and always occurs for its descendant Sphaeroidinella dehiscens (Figs. 4–6). Our analysis shows how the occurrence of a low $\mathrm{Mg}/\mathrm{Ca}$ offset in planktonic foraminifera becomes progressively rarer going back in time, in parallel with the rarity of crust/cortex features (Fig. 8). The occurrence of a crust/cortex is commonly observed in non-spinose modern planktonic foraminifera; however, only two early to middle Miocene genera are known to produce crusts (Globoconella and Paragloborotalia), and only one is known to produce a cortex (Sphaeroidinellopsis) (e.g. Fig. 8).

It is tempting to put the pattern of emergence of $\mathrm{Mg}/\mathrm{Ca}$ offsets in relationship with changes in ocean chemistry and global climate over the last 15 Myr. In particular, the offset spinose species are mostly tropical and evolved during a time when mean ocean pH was more than 0.1 pH unit lower than the preindustrial period (Rae et al., 2021). Further, the Praeorbulina–Orbulina plexus evolved at about 16 Ma (Fig. 8), in coincidence with a drop in ocean pH, likely linked to the global warmth of the Miocene Climatic Optimum (Rae et al., 2021). The particularly high $\mathrm{Mg}/\mathrm{Ca}$ signature of this group of species, along with their evolutionary timing, might testify their ability to withstand tropical surface waters more undersaturated than today thanks to changes in the biomineralization pathway as a consequence of their evolution during the Miocene Climatic Optimum.

The evolution of the offset non-spinose species happened several millions of years later during the long-term cooling trend of the last 10 Myr (Fig. 8). The offset species occur across tropical to high-latitude areas and mixed-layer (Sphaeroidinellopsis, Sphaeroidinella) to intermediate (Neogloboquadrina, Pulleniatina) depth habitats. The ability to develop a crust/cortex in species evolving over the last 10 Myr might have been of advantage as the global ocean was becoming progressively colder, in a similar way to the observed increases in shell mass across Pleistocene glacial cycles (e.g. Zarkogiannis et al., 2019).

The last 10 Myr was also characterized by decreasing concentration in Ca_sw ([Ca²⁺]) in step with global cooling, reaching concentrations approximately half those of the middle Miocene in the Pleistocene (Zhou et al., 2021). Partly as a consequence, $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$ approximately doubled over the past 15 Myr (Evans et al., 2016a; Rosenthal et al., 2022). With decreasing [Ca²⁺] and increasing [Mg²⁺] in seawater over the Neogene (Brennan et al., 2013; Evans et al., 2016a; Zhou et al., 2021), some species might have started to more actively control the $\mathrm{Mg}/\mathrm{Ca}$ ratio at their biomineralization site, e.g. by proportionally increasing the active transport of Mg²⁺, in order to buffer against the effects of the higher seawater $\mathrm{Mg}/\mathrm{Ca}$ and to keep the outer parts of their shell with low $\mathrm{Mg}/\mathrm{Ca}$ and thus more resistant to dissolution. Hence, the low $\mathrm{Mg}/\mathrm{Ca}$ offset observed in the modern and fossil non-spinose species above might be linked to their evolutionary emergence during a time of changing ocean physical–chemical properties, which may have promoted the evolution of thicker tests with a different elemental chemistry making them less buoyant and resistant to dissolution.

5 Conclusions

We analysed a multispecies planktonic foraminiferal δ¹⁸O and $\mathrm{Mg}/\mathrm{Ca}$ dataset spanning the last 15 Myr from multiple locations to test whether temperature is the main controller of both proxies and assess the major overprinting factors through time and space and for species with very distinct ecologies. Once diagenesis and possible regional hydrographic factors are taken into account, we find that species-specific offsets not accounted for in our calibration strategy remain a source of mismatch between the two proxies. Specifically, Globigerina, Globigerinella, Praeorbulina, and Orbulina species are consistently offset, with $\mathrm{Mg}/\mathrm{Ca}$ values higher than expected on average. Conversely, non-spinose Neogloboquadrina, Pulleniatina, Sphaeroidinellopsis, and Sphaeroidinella appear consistently offset with low $\mathrm{Mg}/\mathrm{Ca}$. The appearance of these geochemical offsets appears to be linked to the origination of new clades and is shared between ancestor-descendent species, such that we were able to track their evolutionary history. The variable offset in Globigerinella may go back to the early globigerinids of the Paleogene and could be related to the opportunistic behaviour of this group and emergence of multiple genotypes through geological time, leading to a wider range of habitat conditions. The high $\mathrm{Mg}/\mathrm{Ca}$ offset in Orbulina starts with Praeorbulina in the middle Miocene, while a low $\mathrm{Mg}/\mathrm{Ca}$ offset appears typical of groups evolving in the late Neogene characterized by a crust or cortex. This pattern suggests that the offsets observed in modern species may be a legacy of their parent groups originating millions of years ago, when ocean properties were different from today.

Overall, our study highlights the power of the multispecies and multi-time slice dataset presented here, enabling us to identify the evolutionary origin and timing of deviations in $\mathrm{Mg}/\mathrm{Ca}$–temperature/pH relationships. Furthermore, our study demonstrates the robustness of $\mathrm{Mg}/\mathrm{Ca}$ and δ¹⁸O proxies through geologic time when nonthermal factors (especially $\mathrm{Mg}/{\mathrm{Ca}}_{\text{sw}}$ and pH) are accounted for. For example, virtually all Globigerinoides and Trilobatus $\mathrm{Mg}/\mathrm{Ca}$- and δ¹⁸O-derived temperatures are within uncertainty of each other, highlighting the utility of these species for paleoceanographic reconstruction. In addition, our analysis enables us to identify species/lineages that should be treated with caution when interpreting $\mathrm{Mg}/\mathrm{Ca}$ data, at the very least demonstrating that care should be taken in selecting the calibration approach and highlighting the need for further work in understanding the nonthermal controls on Mg incorporation into the shells of these foraminifera.