A 15-million-year surface- and subsurface-integrated TEX₈₆ temperature record from the eastern equatorial Atlantic

TEX₈₆ is a paleothermometer based on Thaumarcheotal glycerol dialkyl glycerol tetraether (GDGT) lipids and is one of the most frequently used proxies for sea-surface temperature (SST) in warmer-than-present climates. However, GDGTs are not exclusively produced in and exported from the mixed layer, so sedimentary GDGTs may contain a depth-integrated signal that is also sensitive to local subsurface temperature variability. In addition, the correlation between TEX₈₆ and SST is not significantly stronger than that to depth-integrated mixed-layer to subsurface temperatures. The calibration of TEX₈₆ to SST is therefore controversial. Here we assess the influence of subsurface temperature variability on TEX₈₆ using a downcore approach. We present a 15 Myr TEX₈₆ record from Ocean Drilling Program Site 959 in the Gulf of Guinea and use additional proxies to elucidate the source of the recorded TEX₈₆ variability. Relatively high GDGT$[\mathrm{2}/\mathrm{3}]$ ratio values from 13.6 Ma indicate that sedimentary GDGTs were partly sourced from deeper (>200 m) waters. Moreover, late Pliocene TEX₈₆ variability is highly sensitive to glacial–interglacial cyclicity, as is also recorded by benthic δ¹⁸O, while the variability within dinoflagellate assemblages and surface/thermocline temperature records (U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ and $\mathrm{Mg}/\mathrm{Ca}$) is not primarily explained by glacial–interglacial cyclicity. Combined, these observations are best explained by TEX₈₆ sensitivity to sub-thermocline temperature variability. We conclude that TEX₈₆ represents a depth-integrated signal that incorporates a SST and a deeper component, which is compatible with the present-day depth distribution of Thaumarchaeota and with the GDGT$[\mathrm{2}/\mathrm{3}]$ distribution in core tops. The depth-integrated TEX₈₆ record can potentially be used to infer SST variability, because subsurface temperature variability is generally tightly linked to SST variability. Using a subsurface calibration with peak calibration weight between 100 and 350 m, we estimate that east equatorial Atlantic SST cooled by ∼5 ^∘C between the Late Miocene and Pleistocene. On shorter timescales, we use the TEX₈₆ record as a proxy for South Atlantic Central Water (SACW), which originates from surface waters in the South Atlantic Gyre and mixes at depth with Antarctic Intermediate Water (AAIW). Leads and lags around the Pliocene M2 glacial (∼3.3 Ma) in our record, combined with published information, suggest that the M2 glacial was marked by SACW cooling during an austral summer insolation minimum and that decreasing CO₂ levels were a feedback, not the initiator, of glacial expansion.

1 Introduction

Accurate tropical sea-surface temperature (SST) reconstructions are needed to assess oceanic heat distribution and ocean–atmosphere circulation in warmer-than-present climates, such as during the Pliocene. Tropical warming, even of a small magnitude, can lead to major changes in atmospheric circulation, with effects including intensification of monsoonal precipitation (Haywood et al., 2020; Zhang et al., 2016) and warming of the extratropics via teleconnections (Barreiro et al., 2006). At present, the most frequently used proxies for past SSTs are $\mathrm{Mg}/\mathrm{Ca}$ (foraminifer calcite), U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ (alkenones) and TEX₈₆ (glycerol dialkyl glycerol tetraethers; GDGTs). Each proxy has specific confounding factors and calibration issues. For example, $\mathrm{Mg}/\mathrm{Ca}$ paleothermometry requires corrections for the $\mathrm{Mg}/\mathrm{Ca}$ ratio of seawater and diagenesis (Dekens et al., 2002; Evans et al., 2016), and the U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ proxy is insensitive to temperatures above ∼28 ^∘C (Müller et al., 1998). TEX₈₆ is potentially one of the best proxies to reconstruct SSTs above 28 ^∘C, but one critical concern is that GDGTs are produced by marine Thaumarchaeota throughout the water column and often dominantly below the mixed layer (Ingalls et al., 2006; Kim et al., 2015; Lengger et al., 2014; Shah et al., 2008).

Although GDGTs may be dominantly produced at depths below the mixed layer, the efficiency of organic carbon export towards the sediment decreases exponentially with depth below the photic zone as shown with sediment trap data and simulations (e.g., Martin et al., 1987; Middelburg, 2019), which for GDGTs is supported by in situ measurements (e.g., Wüchter et al., 2005). However, as the production of GDGTs in the mixed layer is generally low relative to production within and below the thermocline, GDGTs produced below the mixed layer might still comprise a significant component of the sedimentary assemblage (e.g., Taylor et al., 2013; Ho and Laepple, 2017). The export depth of downcore GDGTs can potentially be traced with the fractional abundance of GDGT-2 to GDGT-3 (GDGT$[\mathrm{2}/\mathrm{3}]$ ratio; Hernández-Sánchez et al., 2014; Pearson et al., 2016; Taylor et al., 2013; Villanueva et al., 2015). GDGT$[\mathrm{2}/\mathrm{3}]$ values are typically low (<5) in suspended particulate matter (SPM) sampled from shallow (<100 m) water and sharply increase to >25 in deeper water (Hurley et al., 2018). This division is currently best explained by depth differentiation of Thaumarcheotal ecotypes with different cyclization patterns (Kim et al., 2016; Villanueva et al., 2015). Core tops do not have GDGT$[\mathrm{2}/\mathrm{3}]$ values indicative of a purely “deep” Thaumarcheotal GDGT source (Fig. 1). Instead, they show a gradual transition from low to moderately high values, which can be attributed to the integration of GDGTs from a range of shallow to intermediate depths. Likewise, TEX₈₆ (Eq. 1) may reflect a depth-integrated surface to subsurface signal (SubST) that is sensitive to temperature variability at and below the surface ocean (Kim et al., 2008, 2012; Schouten et al., 2002; Tierney and Tingley, 2015).

$$\begin{array}{}\text{(1)}& {\text{TEX}}_{\mathrm{86}}={\displaystyle \frac{[\text{GDGT-2}+\text{GDGT-3}+{\text{cren}}^{\prime }]}{[\text{GDGT-1}+\text{GDGT-2}+\text{GDGT-3}+{\text{cren}}^{\prime }]}}\end{array}$$

Indeed, several multi-proxy paleostudies suggest that dissimilarities between TEX₈₆-based SST reconstructions and other temperature records may be explained by a SubST signal recorded in TEX₈₆ variability (e.g., Lopes dos Santos et al., 2010; McClymont et al., 2012; Rommerskirchen et al., 2011; White and Ravelo, 2020). Together with the GDGT$[\mathrm{2}/\mathrm{3}]$ distribution in core tops (Fig. 1), these downcore records support the use of a SubST calibration. However, the targeted depth interval and weight distribution of the temperature integration remain the subject of discussion (Ho and Laepple, 2016, 2017; Tierney et al., 2014). The correlation of core top TEX₈₆ to SST is not significantly stronger than to SubST, as obtained from a wide range of depth intervals (Ho and Laepple, 2016; Schouten et al., 2002). Therefore, the sensitivity of sedimentary TEX₈₆ to temperature cannot be inferred from variability in the spatial dimension, but variability in the temporal dimension may provide solutions. Here, we present a new 0–15 Ma TEX₈₆ record from Ocean Drilling Program (ODP) Site 959 in the eastern equatorial Atlantic. We evaluate the TEX₈₆ and GDGT$[\mathrm{2}/\mathrm{3}]$ evolution on million-year and glacial–interglacial timescales during the late Pliocene. We complement these records with U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ SST estimates and dinoflagellate cyst (dinocyst) assemblages and compare them to the Site 959 benthic δ¹⁸O record (Norris, 1998a; van der Weijst et al., 2020) and $\mathrm{Mg}/\mathrm{Ca}$-based SST and thermocline temperature records (van der Weijst et al., 2022) to elucidate the local source of downcore TEX₈₆ index variability and, consequently, paleoclimatic change.

Figure 1 (a) Depth distribution of GDGT$[\mathrm{2}/\mathrm{3}]$ ratio in suspended particulate matter (SPM; solid lines) at several stations in the Atlantic Ocean (Hurley et al., 2018) and in core tops (circles) between 65^∘ N and 65^∘ S (Tierney and Tingley, 2015). (b) GDGT$[\mathrm{2}/\mathrm{3}]$ vs. TEX₈₆, with logarithmic fits for GDGT$[\mathrm{2}/\mathrm{3}]>\mathrm{5}$ (blue line) and GDGT$[\mathrm{2}/\mathrm{3}]<\mathrm{5}$ (orange line). (c) Depth-integrated (linear) GDGT$[\mathrm{2}/\mathrm{3}]$ from SPM data (station 23; Hurley et al., 2018). Integration start depths at 0, 50, 100 and 200 m and continuous end depth on the x axis. Range of typical Site 959 sedimentary GDGT$[\mathrm{2}/\mathrm{3}]$ values between orange dashed lines and corresponding range of potential integration end depths in the vertical. Note that the assumed linear depth weighing may not be realistic but that an integration start depth at 200 m is incompatible with the sedimentary GDGT$[\mathrm{2}/\mathrm{3}]$ values at Site 959, regardless of depth weighing.

2 Material and methods

2.1 Site, age model and sampling

ODP Site 959 was drilled during Leg 159 in the Gulf of Guinea, ∼160 km offshore Ghana and Côte D'Ivoire (3.62^∘ N, 2.73^∘ W; 2090 m depth; Mascle et al., 1996). It is presently located below the eastward-flowing Guinea Current, which originates from the North Equatorial Counter Current and the Canary Current (Norris, 1998a). SST currently varies between 25.3 and 28.6 ^∘C on a seasonal timescale (Locarnini et al., 2013), and surface salinity ranges between 34.6 and 35.0 (Zweng et al., 2013). The water column is characterized by a shallow thermocline, with the 20 ^∘C isotherm depth annually varying between ∼ 40 and 60 m in response to coastal upwelling (Locarnini et al., 2013; van der Weijst et al., 2022). A minor upwelling event occurs in boreal winter, and a longer and stronger event occurs in boreal summer (Djakouré et al., 2017; Verstraete, 1992; Wiafe and Nyadjro, 2015).

We studied the 0–160 m interval of Hole C and 162–194 m interval of Hole A. The lithology gradually changes from nannofossil/foraminifer chalk at the bottom to nannofossil/foraminifer ooze at the top (Mascle et al., 1996). The age model is based on nannofossil biostratigraphy (Shin et al., 1998) and benthic δ¹⁸O stratigraphy (Wagner, 1998) between 0 and 23 meters below seafloor (m b.s.f.), benthic δ¹⁸O stratigraphy between 33 and 46 mbsf (van der Weijst et al., 2020), astronomical tuning of high-resolution X-ray fluorescence data (Vallé et al., 2016) between 51 and 100 mbsf, and a spline function through planktic foraminifer and nannofossil biostratigraphy between 162 and 194 mbsf (Norris, 1998b; Shafik et al., 1998). We plot the data between 21.35 and 97.43 mbsf on the revised meters composite depth (rMCD) scale of Vallé et al. (2016). Because outside this interval a splice is unavailable, we constructed a revised meters below sea floor (rMBSF) scale in coherence with Vallé et al. (2016), assuming ∼5 % loss of sediment between consecutive cores (Supplement).

The core was sampled at a relatively low resolution between 0 and 15 Ma for biomarker analysis and at a higher resolution in the late Pliocene (2.7–3.5 Ma) interval for multi-proxy reconstructions. A total of 219 samples were extracted for TEX₈₆, of which 105 are closely spaced in the 2.7–3.5 Ma interval. U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ data were generated for 79 of these samples, of which 52 are from the 2.7–3.5 Ma interval. Additionally, we performed palynological assessments on 48 late Pliocene splits of samples that were also used for biomarker analysis.

2.2 Biomarker proxies

At Utrecht University, sediment samples were freeze-dried, and the outer surfaces were removed using a clean knife to prevent contamination. The sediment was finely ground using a mortar and pestle. In the process, a split of coarsely ground material was taken for palynology. Biomarkers were extracted at 100 ^∘C and 7.6×10⁶ Pa using a dichloromethane (DCM) : methanol (9:1, $v/v$) solvent mixture in an accelerated solvent extractor (ASE 350, Dionex). The total lipid extract was separated over an activated Al₂O₃ column into apolar, ketone and polar fractions using hexane : DCM (9:1, $v/v$), hexane : DCM (1:1, $v/v$) and DCM : methanol (1:1, $v/v$) as eluents, respectively. For GDGT analysis, a C46 GDGT (99 ng) standard was added to the polar fraction before it was filtered over a 0.45 µm PTFE filter and dissolved in hexane : isopropanol 99:1 ($v/v$). Measurements were performed on an ultrahigh-performance liquid chromatograph–mass spectrometer (UHPLC-MS) using the method of Hopmans et al. (2016). Samples with a branched to isoprenoid tetraether (BIT) index of >0.3 could be influenced by terrestrially produced GDGTs (Hopmans et al., 2004; Weijers et al., 2006) and were excluded from further analysis (n=3). All samples were within the reliable range of the methane index (MI; Zhang et al., 2011) and ring index (RI; Zhang et al., 2016), suggesting no appreciable inputs from methanogenic archaea or other non-thermal factors. We will present our results using various calibrations, including the conventional approach of calibrating TEX₈₆ to SST, but also using various models using the relation between TEX₈₆ and various depth integrations.

The ketone fraction was dissolved in hexane and analyzed on a gas chromatograph (GC) coupled to a flame ionization detector (GC-FID, Hewlett Packard 6890 series) equipped with a CP-Sil 5 fused silica capillary column (25 m×0.32 mm; film thickness 0.12 µm) and a 0.53 mm precolumn. Samples were injected on-column at 70 ^∘C with helium as a carrier gas and a flow rate of 2 mL min⁻¹. The oven program was as follows: 70 ^∘C for 1 min, then ramped to 130 ^∘C at 20 ^∘C min⁻¹, then to 320 ^∘C at 4 ^∘C min⁻¹, and then held isothermal for 10 min. U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ values were calculated from the fractional abundances of the C_37:2 and C_37:3 alkenones following Prahl and Wakeham (1987) and calibrated to SST using the calibration of Müller et al. (1998). We chose not to adopt the Bayspline approach (Tierney and Tingley, 2018) or other tropical calibrations (e.g., Sonzogni et al., 1997, DSR II) as they may amplify noise into signal because of difficulty in accurately determining the very small C_37:3 peak areas at the high range of SST (Herbert et al., 2020).

2.3 Palynology

Coarsely crushed freeze-dried samples were spiked with a known amount of Lycopodium clavatum spores to quantify absolute palynomorph (Stockmarr, 1972) and treated with 30 % HCl (2×) and 40 % cold HF (2×) to remove carbonates and silicates. The residue was sieved over nylon mesh, and from the 10–250 µm fraction, the lighter organic fraction was separated from the heavy mineral fraction by suspension. The remaining palynological residue was mounted on glass microscopic slides with glycerine jelly. Samples were analyzed using an optical microscope under 400× magnification and approximately 300 dinoflagellate cysts (dinocysts) were counted per sample. Taxonomy follows that of Williams et al. (2017). Capisocysta was only identified at the genus level, because known Capisocysta species are best identified by the number of antapical plates (Head, 1998), which is complicated by the common dissociation of the hypocystal plates. We calculate the relative abundances of a taxon based on the total dinocyst sum. Dinocysts with affinities for cold-temperate waters in the modern ocean (Boessenkool et al., 2001; Zonneveld et al., 2013) were grouped as “cool species” and used for estimating surface water temperature. This species are Ataxodinium choane, Corrudinium devernaliae, Corrudinium harlandii, Impagidinium pallidum, Nematosphaeropsis labyrinthus and Pyxidinopsis reticulata. We also use the $P/(P+G)$ ratio (Versteegh, 1994) as a paleoproductivity index. The $P/(P+G)$ ratio quantifies the number of cysts produced by heterotrophic dinoflagellates, which are characterized by peridinioid-type tabulation (P) over the total number of dinocysts, including those produced by autotrophic or mixotrophic species with a gonyaulacoid-type tabulation (G).

3 Results

3.1 TEX₈₆, U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ and GDGT$[\mathrm{2}/\mathrm{3}]$ evolution at Site 959

We explore our TEX₈₆ results using several calibration approaches in Fig. 2, including an SST calibration using the exponential TEX${}_{\mathrm{86}}^{\mathrm{H}}$ (Kim et al., 2010). It should be noted that this calibration yields significant residuals in tropical regions and suffers from regression dilution (Tierney and Tingley, 2014). We also use a linear calibration model using a Bayesian method, a model specific for the tropical Atlantic and models assuming contributions from GDGTs produced below the mixed layer. Below, we first explore our TEX₈₆ record in the context of other GDGT-based indicators and other data using the TEX${}_{\mathrm{86}}^{\mathrm{H}}$ calibration for SST (Figs. 3–7), while subsurface calibrations are discussed in Figs. 7 and 8.

Figure 2 Site 959 TEX₈₆ record, calibrated using TEX${}_{\mathrm{86}}^{\mathrm{H}}$ and alternative SST calibrations and with the HL16-SubST calibration (Ho and Laepple, 2015), compared to the U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ record. Late Pliocene glacial–interglacial variability in panel (a) and 15 Myr evolution in panel (b).

Between 15 and 11 Ma, TEX${}_{\mathrm{86}}^{\mathrm{H}}$-based SST values fluctuate around an average value of ∼29 ^∘C (Fig. 3), similar to the present-day local annual maximum surface temperature of 28.6 ^∘C. A long-term gradual cooling starts around 11 Ma and is punctuated by two major cooling steps of ∼2 ^∘C each, around 4.9 and 3.4 Ma, in which TEX${}_{\mathrm{86}}^{\mathrm{H}}$ drops well below the present-day local annual minimum surface temperature of 25.3 ^∘C. The cooling trend seems to end in the late Pleistocene, where TEX${}_{\mathrm{86}}^{\mathrm{H}}$ records a ∼4 ^∘C warming between 0.5 Ma and present, with a core top TEX${}_{\mathrm{86}}^{\mathrm{H}}$ of 24.0 ^∘C. U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ is at saturation (∼28 ^∘C) for most of the interval but registers some colder temperatures in the late Pliocene (27–28 ^∘C) and drops below 27 ^∘C after 1.8 Ma. The core top sample registers a temperature of 26.6 ^∘C, which is within calibration error (1.5 ^∘C; Müller et al., 1998) of the 27.5 ^∘C modern mean annual SST (Locarnini et al., 2013). The GDGT$[\mathrm{2}/\mathrm{3}]$ ratio is 7.0 in the core top. In the oldest part of the studied interval, GDGT$[\mathrm{2}/\mathrm{3}]$ values vary between 4 and 5.5. Following an abrupt increase at ∼13.6 Ma, GDGT$[\mathrm{2}/\mathrm{3}]$ varies largely between 5 and 8. After 1.7 Ma, GDGT$[\mathrm{2}/\mathrm{3}]$ persistently increases to >7.

3.2 Pliocene dinocyst assemblages

Palynological assemblages are dominated by dinocysts with only minor amounts of terrestrially derived pollen and spores. Dinocyst preservation is generally good. The assemblages are dominated by Brigantedinium spp. (19 %–60 %), closely followed by Spiniferites spp. (8 %–36 %) (Fig. 4; see the Supplement for a complete overview of the assemblages). The interval between 3.33 and 3.16 Ma was sampled at a higher resolution (∼5 kyr on average) and displays strong variability of the negatively correlated Brigantedinium spp. and Spiniferites spp. on timescales <10 kyr. Other relatively abundant and consistently present taxa are Operculodinium spp., Pentapharsodinium dalei and Impagidinium spp. Peak abundances of Lingulodinium spp. were recorded at ∼3.26 Ma (Lingulodinium machaerophorum) and at 2.77 Ma (Lingulodinium hemicystum). Capisocysta sp. is consistently present in the lower part of the interval, with peak abundances (>20 %) at 3.32 and 3.20 Ma, but has its last occurrence at 3.17 Ma. There is general coherence between TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SST variability and the relative abundance of the “cool species” group, which is dominated by N. labyrinthus. On average, the dinocysts with affinities for colder-temperate waters make up 6 % of the total assemblage, but peak abundances of 10 %–13 % occur during the M2 glacial. The $P/(P+G)$ paleoproductivity index is 0.48 on average, with minimum and maximum values of 0.23 and 0.70, respectively (Fig. 4). $P/(P+G)$ values are primarily driven by the abundance of Brigantedinium spp. The $P/(P+G)$ ratio shows no significant relationship with SST as derived from TEX${}_{\mathrm{86}}^{\mathrm{H}}$ (R²=0.05, p=0.14 in linear regression analysis).

Figure 3 New TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SST (Kim et al., 2010 calibration), U${}_{\mathrm{37}}^{\mathrm{k}\prime }$-SST (Müller et al., 1998 calibration) and GDGT$[\mathrm{2}/\mathrm{3}]$ records at Site 959, compared to benthic δ¹⁸O records from Norris (1998a) and van der Weijst et al. (2020). Shaded interval and arrow indicate annual SST range and mean annual SST in the Gulf of Guinea, respectively (Locarnini et al., 2013).

Figure 4 Relative abundances of the major groups of dinoflagellate cysts, the sum of dinocysts with cold affinities (“Cool species”) and $P/(P+G)$ paleoproductivity index compared to TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SST. Shaded bands indicate TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SST minima during the late Pliocene M2, KM2 and G20 glacials.

4 Discussion

4.1 Identifying the source depth of the TEX₈₆ signal

The 15-million-year TEX${}_{\mathrm{86}}^{\mathrm{H}}$ record at Site 959 shows a cooling from >30 ^∘C between 11 and 15 Ma to <20 ^∘C between 0 and 3.3 Ma (Fig. 3). The U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ record, on the other hand, is near or at proxy saturation (>27 ^∘C) for most of the interval, with some lower values after 1.8 Ma. Assuming that both proxies reflect SST, the relatively small offset between U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ and TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SST offset in the core top sample (2.7 ^∘C) could be explained by a combination of calibration errors and seasonality. However, the average downcore offset between TEX${}_{\mathrm{86}}^{\mathrm{H}}$ and U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ is 5.4 ^∘C and peaks at 7–8 ^∘C during the late Pliocene interglacials (Figs. 5 and 6), which demonstrates that TEX₈₆ underestimates past SST at this site, assuming that U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ is a more accurate representation of SST. This could potentially be resolved by using an alternative regional SST calibration. The BAYSPAR-SST calibration (Tierney and Tingley, 2014) accounts for spatial variability of the TEX₈₆-SST relationship and is the only available calibration that raises absolute SST estimates. However, similar to TEX${}_{\mathrm{86}}^{\mathrm{H}}$ (Kim et al., 2010) and low-latitude Atlantic (Zhang et al., 2018) calibrations, it likely overestimates SST variability (Fig. 2), indicating that the calibration slope is too steep. This is potentially a problem inherent to the calibration of TEX₈₆ to SST and could reflect downcore TEX₈₆ sensitivity to SubST instead of SST variability (Ho and Laepple, 2016). In the following sections, we explore multi-proxy data from Site 959 in search of indications for the water depth the downcore TEX₈₆ signal is reflecting.

Figure 5 Compilation of late Pliocene multi-proxy temperature data and benthic δ¹⁸O at Site 959. $\mathrm{Mg}/\mathrm{Ca}$_(G. ruber) (pink) and $\mathrm{Mg}/\mathrm{Ca}$_{(N. dutertrei)} (yellow) from van der Weijst et al. (2022), benthic δ¹⁸O from van der Weijst et al. (2020), and LR04 δ¹⁸O stack (Lisiecki and Raymo, 2005) for reference. Dashed lines are LOESS-smoothed trends, and vertical bands indicate the M2, KM2 and G20 glacial stages.

4.1.1 GDGT$[\mathrm{2}/\mathrm{3}]$ values

A substantial contribution of GDGTs from below the oceanic surface layer is supported by the relatively high GDGT$[\mathrm{2}/\mathrm{3}]$ ratio, which fluctuates between 5 and 8 in the majority of the record (Fig. 3). Such values are rarely observed in modern-day shallow water (<100 m) SPM and surface sediments (e.g., Besseling et al., 2019; Hernández-Sánchez et al., 2014; Hurley et al., 2018). These values are best explained as a mixed signal of GDGTs from the upper 200 m and intermediate waters (Fig. 1c), i.e., likely dominantly from below the mixed layer (∼50 m) at Site 959. A cross-plot of TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SST and GDGT$[\mathrm{2}/\mathrm{3}]$ (Fig. 7) shows that most data are clustered around a positive regression slope, but the oldest (>13.6 Ma) and youngest (<1.7 Ma) data plot outside this main cluster, potentially signaling systematic changes in GDGT production and export depth. Both 13.6 and 1.7 Ma mark the onset of a strong GDGT$[\mathrm{2}/\mathrm{3}]$ increase (Fig. 3), signaling systematic deepening of the source of the sedimentary GDGT signal.

The ratio of GDGT[2] to GDGT[3] in samples does not directly affect the TEX₈₆ index because they are both included in the denominator and divisor (Eq. 1). If the abundance of GDGT[2] increases relative to that of GDGT[3] with export depth, then samples with a higher GDGT$[\mathrm{2}/\mathrm{3}]$ ratio should be associated with lower TEX₈₆ values. The positive relationship between TEX₈₆ and GDGT$[\mathrm{2}/\mathrm{3}]$ (Fig. 7) is therefore an interesting feature. Taylor et al. (2013) found a similar relation in the calibration dataset and in the paleo-domain. Based on their analyses, they conclude that water depth is the dominant control of GDGT$[\mathrm{2}/\mathrm{3}]$ ratios in sediments, recently supported by SPM analyses (Hurley et al., 2018). They also conclude that the positive correlation between SST and GDGT$[\mathrm{2}/\mathrm{3}]$, notably the opposite relation to that expected based on homeoviscous adaptation, is an oceanographical artifact. Rather, the positive relation suggests that the long-term TEX₈₆ trends are not driven by GDGT export depth but reflect overall temperature changes instead. Alternatively, non-thermal factors such as water column oxygenation and nutrient supply may have influenced GDGT cyclization: high nutrient availability and ammonia oxidation rates have been linked to low TEX₈₆ values (Hurley et al., 2016; Park et al., 2018). However, following the GDGT$[\mathrm{2}/\mathrm{3}]$ shift at 13.6 Ma, TEX${}_{\mathrm{86}}^{\mathrm{H}}$ does not systematically change until ∼11 Ma (Fig. 3). Moreover, late Pliocene dinocyst assemblages (Fig. 4) are characterized by a highly variable $P/(P+G)$ ratio, indicating that upwelling and nutrient supply were highly variable on sub-Milankovitch timescales. In contrast to TEX${}_{\mathrm{86}}^{\mathrm{H}}$ (Fig. 5), the $P/(P+G)$ ratio shows no glacial–interglacial variability. It is therefore unlikely that TEX₈₆ variability was primarily driven by non-thermal factors.

Figure 6 Close-up of G20, KM2 and M2 glacial stages in 100 kyr windows. Multi-proxy temperature data and benthic δ¹⁸O as in Fig. 5. The resolution of the dinocyst record was too low for meaningful comparisons in KM2 and M2 glacials.

4.1.2 Surface, thermocline and sub-thermocline temperature variability

The late Pliocene TEX${}_{\mathrm{86}}^{\mathrm{H}}$ record follows a different evolution than SST records of U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ and the recent $\mathrm{Mg}/\mathrm{Ca}$ of the surface-dwelling Globigerinoides ruber ($\mathrm{Mg}/\mathrm{Ca}$_(G. ruber)) and thermocline-dwelling Neogloboquadrina dutertrei ($\mathrm{Mg}/\mathrm{Ca}$_{(N. dutertrei)}) generated on the same samples (Fig. 5). The $\mathrm{Mg}/\mathrm{Ca}$ records (van der Weijst et al., 2022) are based on very well-preserved foraminifera and apply a constant Pliocene correction for the $\mathrm{Mg}/\mathrm{Ca}$ of seawater, a species-specific core top calibration with a constant partitioning coefficient and embedded dissolution correction (Dekens et al., 2002). Whereas TEX${}_{\mathrm{86}}^{\mathrm{H}}$ decreases by ∼2 ^∘C, U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ and $\mathrm{Mg}/\mathrm{Ca}$-based SSTs show no significant net change in this interval. Moreover, although $\mathrm{Mg}/\mathrm{Ca}$ and U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ register some cooling during the late Pliocene glacials, their variability seems not primarily driven by glacial–interglacial cyclicity, as is best observed in the high-resolution interval from M2 to KM2 (Fig. 5). This is also true for the relative abundance of dinocysts with affinities for colder temperatures (Figs. 4 and 5). Furthermore, around the KM2 and M2 glacials, TEX${}_{\mathrm{86}}^{\mathrm{H}}$ cooling leads to increasing abundances of cool dinocysts, as well as $\mathrm{Mg}/\mathrm{Ca}$_(G. ruber) and U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ cooling by ∼ 5–10 ka (Fig. 6), which underscores the independent evolution of these records.

At sites where TEX₈₆ is suspected to be affected by GDGTs produced below the surface ocean, TEX${}_{\mathrm{86}}^{\mathrm{H}}$ has occasionally been interpreted to reflect thermocline temperature variability (e.g., Lopes dos Santos et al., 2010). Thermocline temperatures are lower and can be highly variable through time, which potentially explains both absolute SST underestimation and amplitude overestimation in specific TEX${}_{\mathrm{86}}^{\mathrm{H}}$ records. However, both could also be explained as a general artifact of an overestimated calibration slope that results from calibrating a depth-integrated GDGT signal to SST (Ho and Laepple, 2015, 2016). Our TEX${}_{\mathrm{86}}^{\mathrm{H}}$ record shows weak coherence with the Site 959 $\mathrm{Mg}/\mathrm{Ca}$_{(N. dutertrei)} thermocline temperature record (Fig. 5; van der Weijst et al., 2022). Whereas TEX₈₆ registers late Pliocene cooling, $\mathrm{Mg}/\mathrm{Ca}$_{(N. dutertrei)} registers warming in relation to thermocline deepening (van der Weijst et al., 2022). Moreover, while the range of $\mathrm{Mg}/\mathrm{Ca}$_{(N. dutertrei)} temperatures is similar to TEX${}_{\mathrm{86}}^{\mathrm{H}}$ on glacial–interglacial timescales, $\mathrm{Mg}/\mathrm{Ca}$_{(N. dutertrei)} does not follow the glacial–interglacial cyclicity as recorded in the TEX₈₆.

In contrast, the variability in the TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SST and benthic δ¹⁸O records at Site 959 is very similar (Figs. 3 and 5). Benthic δ¹⁸O is a faithful recorder of glacial–interglacial cyclicity because it registers a combined signal of deep ocean temperature and polar ice sheet volume. Could TEX₈₆ at Site 959 also be sensitive to glacial–interglacial cyclicity in deeper waters? At the onset of the M2 and KM2 glacials, the glacial expression of TEX${}_{\mathrm{86}}^{\mathrm{H}}$ leads δ¹⁸O by ∼5 kyr (Fig. 6), which indicates that, despite the similarities between TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SST and benthic δ¹⁸O, these proxies record variability in different water masses. Bottom waters at Site 959 are predominantly ventilated by North Atlantic Deep Water (NADW), whereas depths below the thermocline are occupied by South Atlantic Central Water (SACW; Fig. 8). SACW originates from surface waters in the South Atlantic gyre, which at depth mixes with Antarctic Intermediate Water (AAIW). It is therefore sensitive to climate variability in the mid-latitudes in the Southern Hemisphere and, like benthic δ¹⁸O, sensitive to higher-latitude climate change. Although our record samples SACW, AAIW dynamics have previously been suggested to be the main driver of TEX₈₆ variability in the Arabian Sea (Huguet et al., 2006) and in the southeast Atlantic (Rommerskirchen et al.,  2011). If TEX₈₆ at Site 959 records a mixed signal from surface and subsurface waters, as is supported by the relatively high GDGT$[\mathrm{2}/\mathrm{3}]$ ratio (Figs. 1 and 3), it is likely sensitive to SACW dynamics, which explains the strong coherence with the benthic δ¹⁸O record. Collectively, these observations suggest that the downcore TEX₈₆ record at Site 959 is substantially affected by temperature variability below the surface ocean and that calibration to SubST is more appropriate than to SST.

Figure 7 Cross-plot of TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SST and GDGT$[\mathrm{2}/\mathrm{3}]$. The oldest (>13.6 Ma; dark blue) and youngest (<1.7 Ma; dark red) samples largely stand out from the majority of the data. Linear regression lines determined from full dataset (stippled line) and samples dated between 13.6 and 1.7 Ma (solid line; 88 % samples, R²=0.37, p<0.001).

4.2 Calibration of the depth-integrated TEX₈₆ record

Several TEX₈₆ calibrations for subsurface temperature (TEX₈₆-SubST) are available, with each assuming, among other things, a different depth integration of the water column (Fig. 8). The BAYSPAR-SubST (Tierney and Tingley, 2014) and TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SubST (Kim et al., 2012) calibrations both target the upper 200 m of the water column, but the weight of the BAYSPAR-SubST calibration is centered at relatively shallow depths compared to the linearly weighted TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SubST calibration (Fig. 8). The HL16-SubST (Ho and Laepple, 2016) calibration targets the upper 950 m. The validity of this depth interval was questioned (Tierney et al., 2017), but with peak calibration weight at 100–350 m (Fig. 8), it is compatible with the documented vertical distribution of Thaumarcheotal cell counts and GDGT concentrations (e.g., Hernández-Sánchez et al., 2014; Schouten et al., 2012; Wuchter et al., 2005; Sintes et al., 2016), leaving export efficiency as the major uncertainty. Moreover, based on the large proportion of core tops with higher (>5) GDGT$[\mathrm{2}/\mathrm{3}]$ values in the global calibration set (Fig. 1), a slightly deeper calibration target (including depths below 200 m) is arguably preferable. This deeper calibration target lowers the slope of the calibration, which dampens reconstructed Miocene–Pliocene temperature variability at Site 959 (Fig. 2).

Figure 8 (a) Temperature and salinity distribution along a depth transect in the modern Gulf of Guinea (numbers: meters below sea surface), data from Locarnini et al. (2013) and Zweng et al. (2013). (b) Vertical weight distribution of the discussed SubST calibrations (silhouettes) and the water column composition in the modern Gulf of Guinea. (c) Site 959 TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SubST evolution according to different SubST calibrations, compared to global mean temperature evolution (Tierney et al., 2020). Δ Temperature estimated from linear trends (dashed lines) between the Middle Miocene and Pleistocene.

The ratio between temperature change in the surface and subsurface ocean is 1:1 when averaged across many sites and on longer timescales (Ho and Laepple, 2016). Under the assumption that the (integrated) depth of GDGT sourcing and the structure of the water column was stable, a depth-integrated TEX₈₆ record might still be of value to assess SST variability. Depending on the calibration, the magnitude of Late Miocene to Pleistocene cooling at Site 959 is 10 ^∘C (BAYSPAR-SubST), 7 ^∘C (TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SubST) or 5 ^∘C (HL16-SubST; Fig. 8). Assuming a 1:1 ratio, the long-term SST trend should be of a similar magnitude. Currently available tropical SST estimates on the studied timescale are either compromised by saturation of the U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ proxy or based on TEX₈₆ and are therefore not suitable for independent comparison. However, global mean surface temperature estimates based on benthic δ¹⁸O data suggest a ∼5 ^∘C cooling across the same interval (Hansen et al., 2013; Tierney et al., 2020). In the absence of major oceanographic changes, it is unlikely that the local cooling trend at Site 959 was larger, because temperature variability is generally lower in the tropics compared to high latitudes where deep-ocean waters derive from (i.e., polar amplification). This suggests that the BAYSPAR-SubST and TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SubST calibrations may overestimate the magnitude of long-term cooling (Fig. 8). From the currently available calibrations, HL16-SubST best approximates both expected glacial–interglacial variability and multi-million year trends (Figs. 2 and 8) and could potentially be used in combination with other SST records and/or modern water column data to correct for the local SST-SubST offset (Fig. 2). It is possible that the 1:1 ratio between SST and SubST did not persist across the entire record, e.g., due to changes in SACW production. To further improve tropical SST estimates in warmer climates, it should be further explored under which conditions a TEX₈₆ record that is driven by subsurface variability can be used to reconstruct SST variability and how to correct for the surface–subsurface temperature offset.

4.3 Late Neogene TEX₈₆ as an intermediate ocean signal at Site 959: new insights into M2 glacial inception

It is currently unclear what caused M2 glacial inception, but hypothesized mechanisms include declining atmospheric CO₂ levels (Berends et al., 2019; Dolan et al., 2015) and reduced latitudinal heat transport by the North Atlantic Current (De Schepper et al., 2009, 2013) in the Northern Hemisphere (NH) and by Indonesian throughflow (De Vleeschouwer et al., 2018) in the Southern Hemisphere (SH). Multi-proxy data from Site 999 in the Caribbean Sea show that declining CO₂ levels lagged the benthic δ¹⁸O glacial expression at the onset of M2 (de la Vega et al., 2020). This indicates that declining CO₂ levels did not initiate the M2 glacial, although they might have affected its intensity and duration.

Figure 9 Exploring the chronology around the M2 glacial with records aligned on peak glacial benthic δ¹⁸O values (original age models in Fig. S2). (a) Daily insolation at 60^∘ N/S during NH/SH summer solstice (Laskar, 1990). (b) East equatorial Atlantic SST at Site 959 (this study) and Site 662 (Herbert et al., 2010) and Caribbean Sea SST at Site 999 (de la Vega et al., 2020). (c) Site 959 TEX₈₆ HL16-SubST, indicative of intermediate ocean temperatures. (d) Atmospheric CO₂ reconstructions from Site 999 (de la Vega et al., 2020). (e)  Benthic δ¹⁸O records of Sites 959 (van der Weijst et al., 2020), 662 (Lisiecki and Raymo, 2005) and 999 (Haug et al., 2001; de la Vega et al., 2020), aligned to LR04 (Lisiecki and Raymo, 2005) on peak M2 values.

Because the TEX₈₆ record at Site 959 is affected by SACW, it can potentially be used to study the connection between our low-latitude site and higher latitudes. The TEX₈₆ and benthic δ¹⁸O records at Site 959 correspond closely but are out of sync during the M2 glacial (Fig. 6). Both records, generated on the same samples, display a double peak around M2, but TEX₈₆ leads δ¹⁸O by ∼5 kyr. Whereas the depth-integrated TEX₈₆ signal at Site 959 is sensitive to SH conditions through SACW, bottom waters at Site 959 consisted predominantly of NADW sourced from the NH during the late Pliocene (van der Weijst et al., 2020). The TEX₈₆ lead therefore indicates a lead of SH over NH cooling/glaciation at the onset of M2, in agreement with transient ice sheet simulation (Berends et al., 2019).

We further explore the relative chronology around M2 in Fig. 9. The age models were aligned with peak glacial δ¹⁸O values, and the same records are plotted in the original age models in Fig. S1 in the Supplement. Although it is difficult to discriminate between signal and noise at this resolution, some patterns seem to emerge from this compilation. Directly before M2, benthic δ¹⁸O, SubST and SST records from the eastern equatorial Atlantic (Sites 959 and 662) and Caribbean Sea (Site 999) indicate gradual, widespread cooling, while CO₂ levels were rising. A sharp drop in TEX${}_{\mathrm{86}}^{\mathrm{H}}$-SubST indicates intermediate ocean cooling at the onset of M2, leading pronounced glacial cooling (benthic δ¹⁸O) of the deep ocean and tropical surface ocean (best observed at Site 662). According to the chronology in Fig. 9, this intermediate ocean cooling occurred during a SH summer insolation minimum, suggesting that astronomical forcing may have amplified Southern Ocean cooling and tipped the system to glacial conditions. The sensitivity of the Antarctic ice sheet to astronomical forcing is also reflected by the 3.3 Ma onset of ∼20 kyr cyclicity in ice-rafted debris (IRD) at East Antarctic IODP Site U1361, with peak IRD mass accumulation rates at SH summer insolation minima (Patterson et al., 2014). Based on the chronology of Fig. 9, SH cooling led atmospheric CO₂ drawdown by ∼ 15–20 kyr. Lags of a similar magnitude were not unusual during Pleistocene glacial inceptions (Petit et al., 1999). Marginal Antarctic diatom assemblages and lithology suggest that considerable sea ice expansion occurred during M2, which could have inhibited CO₂ outgassing from the deep ocean (Ishino and Suto, 2020; McKay et al., 2012). Moreover, a slight increase in dust-mediated iron fertilization of the Southern Ocean may have facilitated export productivity (Martínez-Garcia et al., 2011). In combination with reduced mixing between Northern Component Water and Southern Component Water (van der Weijst et al., 2020), these processes would have led to increased carbon storage in the deep ocean in response to M2 cooling. Higher-resolution records and correlations would be required to fully test this hypothesis.

5 Conclusions

Several lines of evidence show that the TEX₈₆ record at Site 959 in the eastern equatorial Atlantic is best explained as a depth-integrated signal that is substantially affected by temperature variability below the thermocline. Relatively high GDGT$[\mathrm{2}/\mathrm{3}]$ ratios since ∼13.6 Ma indicate that the sedimentary GDGTs were partly sourced from deeper (>200 m) waters. Moreover, late Pliocene multi-proxy data show that TEX₈₆ variability is predominantly driven by glacial–interglacial variability, in contrast to other SST and thermocline records (U${}_{\mathrm{37}}^{\mathrm{k}\prime }$ and $\mathrm{Mg}/\mathrm{Ca}$_(G. ruber) SST and $\mathrm{Mg}/\mathrm{Ca}$_{(N. duterrei)}) and dinocyst abundances. The TEX₈₆ record strongly resembles the benthic δ¹⁸O record, indicating that it dominantly records temperature changes in a sub-thermocline water mass. At Site 959, intermediate depths are occupied by SACW, which derives from the mid-latitudes of the South Atlantic. According to the present-day water column composition, a substantial contribution of GDGTs from deeper (>200) waters is needed to obtain TEX₈₆ sensitivity to SACW. This favors a subsurface calibration that integrates across a wider range of water depths, such as HL16-SubST (Ho and Laepple, 2016), in which calibration weight peaks 100–350 m. This calibration target interval is compatible with pelagic Thaumarcheotal cell counts and GDGT concentrations and with core top GDGT$[\mathrm{2}/\mathrm{3}]$ values in the global calibration set.

Even if TEX₈₆ mainly reflects SubST variability, it may be used to reconstruct past SSTs if the temporal relationship between SST–SubST in a certain region is well-understood. Assuming a 1:1 relationship between long-term SubST and SST trends, our TEX₈₆ record suggests 5 ^∘C tropical SST cooling between the Late Miocene and Pleistocene when calibrated with the Ho and Laepple (2016) subsurface calibration. Additionally, TEX₈₆ is also highly informative as a SubST proxy, because the sensitivity of TEX₈₆ to SACW at Site 959 offers a chance to explore connections between high- and low-latitude climate variability. The TEX₈₆ and benthic δ¹⁸O records at Site 959 are highly similar, but TEX₈₆ is sensitive to high-latitude SH conditions through AAIW, whereas the benthic δ¹⁸O record is primarily forced by high-latitude NH conditions through NADW. A ∼5 kyr lead of TEX₈₆ relative to δ¹⁸O at the onset of the M2 glacial stage indicates a lead of SH over NH cooling. Albeit limited by the resolution of the datasets, multi-proxy data from Site 959, Site 999 (Caribbean Sea) and Site 662 (eastern equatorial Atlantic) aligned based on LR04 peak glacial δ¹⁸O values indicate that SH cooling also led tropical ocean cooling and CO₂ levels by 15–20 kyr, suggesting that CO₂ drawdown was a consequence of glacial conditions. Instead, glacial expansion and SACW cooling at SH insolation minima suggest that orbital forcing played a pivotal role in M2 glacial inception.