A Late Quaternary climate record based on long-chain diol proxies from the Chilean margin

In this study we have applied different indices based on long-chain diols, i.e., the long-chain diol index (LDI) as a proxy for past SST, the diol index as an indicator of past upwelling conditions, and the nutrient diol index (NDI) as a proxy for nitrate and phosphate concentrations in seawater. The proxies were analyzed in marine sediments recovered at ODP Site 1234, located within the Peru–Chile upwelling system, with a ∼2 kyr resolution covering the last 150 kyr. We also generated TEX${}_{\mathrm{86}}^{\mathrm{H}}$ and U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ temperature and planktonic δ¹⁸O records, as well as total organic carbon (TOC) and accumulation rates (ARs) of TOC and lipid biomarkers (i.e., C₃₇ alkenones, GDGTs, dinosterol, and loliolide) to reconstruct past phytoplankton production. The LDI-derived SST record covaries with TEX${}_{\mathrm{86}}^{\mathrm{H}}$- and U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$-derived SST records as well as with the planktonic δ¹⁸O record, implying that the LDI reflects past SST variations at this site. TOC and phytoplankton AR records indicate increased export production during the last interglacial (MIS 5), simultaneous with a peak in the abundance of preserved Chaetoceros diatoms, suggesting intensified upwelling during this period. The diol index is relatively low during the upwelling period, but peaks before and after this period, suggesting that Proboscia diatoms were more abundant before and after the period of upwelling. The NDI reveals the same trends as the diol index, suggesting that the input of nitrate and phosphate was minimal during upwelling, which is unrealistic. We suggest that the diol index and NDI should perhaps be considered as indicators for Proboscia productivity instead of upwelling conditions or nutrient concentrations.

1 Introduction

Paleoclimatic reconstructions typically rely on physical, biological, and geochemical proxies from sedimentary archives. Physical proxies include sediment composition, structure, grain size, density, and magnetic susceptibility, which may provide information on the paleodepositional environment. Biological proxies comprise preserved (micro)organisms such as diatoms, foraminifera, dinoflagellates, corals, and mollusks or remnants deriving from higher plants such as pollen and spores. Geochemical proxies are based on the chemical composition of either the sediment or fossilized organisms. Stable isotope and elemental concentrations are inorganic proxies typically measured in shells or skeletons of marine organisms (foraminifera, mollusks, corals) preserved in the sediment, providing insight into the chemistry of the seawater in which the organisms lived. Organic proxies are a relatively newer class of tools based on fossilized molecules that are unique for a specific organism or group of organisms, referred to as biomarkers. The ratios of specific biomarker molecules are often related to physical parameters of the environment in which the source organism grew, such as temperature, salinity, oxygen availability, or productivity, and therefore such organic proxies can also be applied to reconstruct past depositional environments (e.g., Brassell et al., 1986; Prahl and Wakeham, 1987; Schouten et al., 2002; Rampen et al., 2008, 2012; Willmott et al., 2010; Gal et al., 2018).

In the last decade long-chain diols (LCDs) have attracted attention as novel proxies to reconstruct past environmental conditions. The long-chain diol index (LDI) was proposed based on the distribution of C₂₈ and C₃₀ 1,13- and 1,15-diols in marine surface sediments, which shows a good correlation with mean annual sea surface temperature (SST; Rampen et al., 2012). These compounds have been detected in cultures of Eustigmatophyte algae (e.g., Volkman et al., 1992, 1999; Méjanelle et al., 2003; Rampen et al., 2014a), but the LCD distributions observed are different to those observed in the marine realm (e.g., Versteegh et al., 1997, 2000; Rampen et al., 2012, 2014a), and hence their role as LCD producers in the ocean remains uncertain. The diol index, which is an indicator for past upwelling and high nutrient conditions, is defined as the fractional abundance of 1,14-diols with respect to 1,13-diols (Willmott et al., 2010) or the C₃₀ 1,15-diol (Rampen et al., 2008), and 1,14-diols are biomarkers for Proboscia diatoms (Sinninghe Damsté et al., 2003; Rampen et al., 2014b). Proboscia grows during the early stages of upwelling since they need relatively little silica and they are able to migrate to deeper waters to obtain nutrients (Koning et al., 2001); therefore, it was proposed that the relative abundance of 1,14-diols can trace past upwelling conditions. C₂₈, C₃₀, and C₃₂ 1,14-diols have also been observed in the marine dictyochophyte Apedinella radians (Rampen et al., 2011) but its role as a 1,14-diol producer in the marine realm is unknown. Recently, a new index based on the saturated and monounsaturated C₂₈ 1,14-diol was proposed as a quantitative proxy for nitrate and phosphate concentrations (Gal et al., 2018), called the nutrient diol index (NDI). The authors found a strong positive correlation between the NDI and phosphate (R²=0.85) and nitrate (R²=0.80) concentrations for the marine sediment datasets of Rampen et al. (2014b) and de Bar et al. (2016). This suggests that the NDI might be a good indicator of past nutrient variations in surface waters.

Applications in sediment cores of these LCD proxies have shown that the LDI, as well as the diol index, is a promising paleotemperature and paleo-upwelling indicator (e.g., Pancost et al., 2009; Lopes dos Santos et al., 2012; Rampen et al., 2012; Seki et al., 2012; Rodrigo-Gámiz et al., 2014; Plancq et al., 2015; Jonas et al., 2017). However, a number of uncertainties still exist in the application of these biomarkers. For example, recent studies of surface sediments from coastal regions reveal different 1,13- and 1,15-diol distributions compared to open ocean sediments (de Bar et al., 2016; Lattaud et al., 2017a, b). Relatively high fractional abundances of the C₃₂ 1,15-diol along the coast were observed as a result of riverine input, significantly affecting the LDI signal, likely due to different 1,13- and/or 1,15-diol producers thriving in river outflow waters. Moreover, studies have related high Proboscia diatom abundances to stratified instead of upwelling conditions (e.g., Fernández and Bode, 1994). Similarly, Contreras et al. (2010) determined the relative abundance of the C₂₈ 1,14-diol at ODP Site 1229 located in the Peruvian upwelling system, and observed increased concentrations during the last interglacial (LIG; MIS 5e), and related this to enhanced water column stratification. Furthermore, Rodrigo-Gámiz et al. (2015) showed that for sediment traps and surface sediments around Iceland, which were characterized by high concentrations of 1,14-diols (> 75 % of all LCDs), the LDI did not correspond to SST. This makes the application of the LDI at sites with high input of diols from Proboscia diatoms (e.g., upwelling sites) uncertain. Finally, due to its recent development no studies have been performed yet to test the applicability of the NDI as a paleo-nutrient proxy.

To constrain the uncertainties in applying the LCD proxies (LDI, diol index, and NDI) at sites with upwelling and riverine input, we studied the Late Quaternary sedimentary record at ODP Site 1234. This site is located along the Chilean margin within the Chile–Peruvian upwelling system and near two major river mouths of the Andean river systems Río Biobío and Río Itata, both draining large basins (Muratli et al., 2010a). We sampled the last ∼150 kyr, including glacial and interglacial periods, and generated LDI-,TEX₈₆-, and U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$-based temperature records to constrain glacial–interglacial variations in SST throughout this entire interval. Additionally, we compared the diol index and the NDI record with other paleoproductivity indicators, including total organic carbon (TOC), organic matter stable carbon isotopes (δ¹³C) and phytoplanktonic lipid biomarkers that are characteristic for certain phytoplankton communities (C₃₇ alkenones, loliolide, and dinosterol).

2 Materials and methods

2.1 Study site

Marine sediments along the coast of Chile and Peru have been thoroughly studied, as it is a key region in the Southern Hemisphere for studying climate variability related to both atmospheric and oceanographic circulation (e.g., Lamy et al., 1998, 1999, 2002, 2004; Hebbeln et al., 2000, 2002; Mohtadi and Hebbeln, 2004; Stuut and Lamy, 2004; Heusser et al., 2006; Romero et al., 2006; Mohtadi et al., 2008; Muratli et al., 2010a, b; Verleye and Louwye, 2010; Chase et al., 2014). The main circulation patterns include the southern westerly winds (SWWs) and the Antarctic Circumpolar Current (ACC; Fig. 1). The ACC approaches the South American continent from the west and initiates both the Humboldt Current (or Peru–Chile Current) flowing northwards along the continental margin and the Cape Horn Current flowing southwards (e.g., Stuut et al., 2006; Fig. 1). The Humboldt Current flows along the Chilean coast and turns westwards as it approaches the Equator, forming the South Equatorial Current (SEC). Variations in the strength and location of the ACC and southern westerlies are thought to be the main climate controls in this region. Lamy et al. (2002) showed that long-term Holocene trends in temperature and productivity are linked to latitudinal shifts of the ACC and the westerlies. On a millennial to multidecennial scale, paleotemperature variations are strongly connected to millennial-scale climate variations in Antarctica. When the ACC migrates northward, potentially associated with an expansion of Antarctic sea ice, cold subantarctic waters are advected into the Southern Hemisphere midlatitudes, including the Chilean margin (e.g., Lamy et al., 2004; Kaiser et al., 2005). The subtropical high-pressure system is an important climate control north of 33^∘ S, which can enhance coastal upwelling (Mohtadi and Hebbeln, 2004).

Figure 1 Maps showing the location of ODP Site 1234. (a) A simplified illustration of the present-day atmospheric and oceanographic setting and the site location. SEC: South Equatorial Current; APF: Antarctic Polar Front; STF: subtropical front. (b) Present-day mean annual SST (^∘C) for the region. Data derived from the World Ocean Atlas 2013 (Locarnini et al., 2013); maps were drawn in Ocean Data View (Schlitzer, 2015) and modified manually.

ODP Site 1234 is located in the southeast Pacific Ocean (36^∘13.153^′ S, 73^∘40.902^′ W; 1015 m water depth), ca. 65 km offshore of Concepción, Chile. This site lies in the vicinity of two large Andean river systems: the Río Biobío and Río Itata, which drain basins of 24 000 and 11 200 km², respectively (Muratli et al., 2010a). The continental margin is under the influence of terrestrial input due to rainfall and land erosion (Lamy et al., 2004), and terrestrial material deposited along the coastal margin derives mainly from two mountain ranges, i.e., the Coast Range and the Andes (Lopez-Escobar et al., 1977; Martin et al., 1999). Moreover, ODP Site 1234 is located within the Peru–Chile upwelling system, which is one of the most important upwelling regions on Earth (Berger et al., 1987) due to the combination of upwelling-favorable winds (Daneri et al., 2000; Carr and Kearns, 2003) and an eastern boundary current (Humboldt Current). The upwelling regime stretches from 5 to 38^∘ S, corresponding to ca. 5000 km of coastline. South of 38^∘ S, onshore winds suppress upwelling (Strub et al., 1998), but high productivity is sustained by relatively high nutrient supply from the ACC (Hebbeln et al., 2000), which is typically rich in nitrate and limited in micronutrients such as iron (de Baar et al., 1995). Along the Chilean coast, iron availability (as well as other micronutrients) can increase due to fluvial supply as the precipitation of the southern westerlies onshore enhances river runoff (Mohtadi and Hebbeln, 2004).

2.2 Age model and sample strategy

The core was drilled as part of Ocean Drilling Program (ODP) Leg 202. Sedimentation rates were relatively high (∼80 cm kyr⁻¹ on average; Mix et al., 2003; Heusser et al., 2006). The age model used is from Heusser et al. (2006) and is based upon radiocarbon dates and the correlation of benthic δ¹⁸O with deep Atlantic core MD95204 (upper half) and the Vostok ice core chronology (lower half; see Heusser et al., 2006 and references therein). Heusser et al. (2006) switched between the two age models around 69 ka where they overlap, as indicated by the dashed line in Fig. 2a. However, this results in a substantial, and likely unrealistic, dip in the sedimentation rate between ca. 80 and 90 ka in the hiatus between the age models (Fig. 2b). Hence, we chose to use the upper age model (dark blue in Fig. 2), which extends up to ∼83 ka, and linearly interpolate between this age and the age of ∼90 ka constrained by the other age model for the lower half of the core. This depth–age range of linear interpolation between the two age models is highlighted in green in Fig. 2 and results in less abrupt changes in sedimentation rates. The core was sampled with a ∼2 kyr resolution covering the last ∼150 kyr. In total, 74 sediment samples were analyzed for bulk and organic geochemistry.

Figure 2 Age–depth relations for ODP 1234 derived from Heusser et al. (2006). (a) Age–depth constraints based on the correlation of benthic δ¹⁸O with Atlantic core MD95204 (dark blue; upper half) and the Vostok ice core chronology (orange; lower half). Heusser et al. (2006) switched between the two different age models at ∼69 ka and ∼53 mcd (dashed line). However, we used the full upper age model and then linearly interpolated between the age models from ∼60 to 64 mcd, as indicated by the connecting green line. (b) Sedimentation rates were calculated after linear interpolation between the age–depth tie points of the two age models. The sedimentation rates calculated after the linear interpolation between the age models are presented in green. In transparent red, the sedimentation rates are plotted when applying the dating strategy of Heusser et al. (2006).

2.3 Elemental analysis

All 74 sediment samples were freeze-dried and homogenized. Small aliquots (ca. 50–100 mg) were used for elemental analysis. For this purpose, all aliquots were acidified with 2 M hydrochloric acid (HCl) to remove all carbonates and rinsed with distilled water to remove salts. Subsequently, the decalcified sediment samples were analyzed for total organic carbon (TOC), total nitrogen (TN), and stable carbon isotope ratios (δ¹³C) by means of a Thermo Scientific Flash 2000 elemental analyzer coupled to a Thermo Scientific Delta V Advantage isotope ratio mass spectrometer. Total nitrogen (TN) was measured on aliquots that were not acidified. Results are expressed in standard δ notation relative to Vienna Pee Dee Belemnite (VPDB) for δ¹³C. The precision as determined using laboratory standards calibrated to certified international reference standards was in all cases < 0.2 ‰. TOC mass accumulation rates (ARs; MAR_TOC) were calculated by multiplication of the sedimentation rate (linear interpolation between depth points of the age–depth model; see above) with an estimated bulk density of 1.6 g cm⁻³ (Mix et al., 2003) and subsequent multiplication with the TOC percentage.

2.4 Foraminiferal stable isotope analysis

Benthic stable isotope data were previously published and were generated using standard techniques in laboratories at Oregon State University and the Woods Hole Oceanographic Institution (McManus et al., 1999, 2002, 2003; Heusser et al., 2006). Oxygen isotope ratios in Cibicidoides and Uvigerina were adjusted to each other by 0.64 ‰ (Shackleton, 1974). Additional data on the planktonic species Globigerina bulloides were generated at Lamont Doherty Earth Observatory (LDEO) of Columbia University using a Thermo Delta V Plus gas-source isotope ratio mass spectrometer (IRMS) equipped with a Kiel IV automated individual acid bath sample-preparation device. Although a recent study found no size-related influence on G. bulloides δ¹⁸O (Costa et al., 2017), specimens were generally picked from the 250–300 µm size range, with 8–12 individuals selected for analysis. Samples were replicated at ∼5 %–10 % frequency and measured isotope ratios were calibrated to the VPDB isotope scale with NBS-19 and NBS-18. Reproducibility of the in-house standard (1σ) is ±0.06 ‰ (δ¹⁸O) and ±0.04 ‰ (δ¹³C).

2.5 Lipid extraction and organic geochemical analysis

The sediment samples (ca. 15 g dry weight) were extracted using accelerated solvent extraction (ASE) using a DIONEX 200 at 100 ^∘C, a pressure of 7–8×10⁶ Pa, and a mixture of dichloromethane (DCM) and methanol (MeOH) (9:1; v:v). The total lipid extracts (TLEs) were dried under a stream of nitrogen (N₂) using a Caliper TurboVap LV. All TLEs were desulfurized using copper granules activated with 1 M HCL. The copper turnings were added to the TLEs, stirred overnight, subsequently dried over anhydrous sodium sulfate (Na₂SO₄) in order to remove precipitate and water, and dried down under N₂. For quantification purposes, three internal standards were added to the TLEs : 10-nonadecanone (C_19:0 ketone) for long-chain alkenones, C₂₂ 7,16-diol for LCDs (Rodrigo-Gámiz et al., 2015), and the C₄₆ GDGT for GDGTs (Huguet et al., 2006). The TLEs (aliquots of ∼4.5 mg) were separated into apolar, ketone (containing alkenones), and polar (containing LCDs and GDGTs) fractions by separation over activated (2 h at 150 ^∘C) Al₂O₃ and elution with hexane : DCM (9:1; v:v), hexane : DCM (1:1; v:v), and DCM : MeOH (1:1; v:v), respectively. Polar fractions were split for GDGT (25 %) and LCD (75 %) analysis.

2.5.1 GDGTs

Aliquots of the polar fractions were dissolved in hexane : isopropanol (99:1, v:v) to a concentration of ca. 2 mg mL⁻¹. The fractions were then filtered through 0.45 µm polytetrafluoroethylene (PTFE) filters. GDGTs were analyzed by means of ultrahigh-performance liquid chromatography mass spectrometry (UHPLC-MS) on an Agilent 1260 HPLC equipped with an automatic injector, which was coupled to a 6130 Agilent MSD, and HP Chemstation software according to Hopmans et al. (2016). The injection volume was 10 µL. Separation of the GDGTs was achieved in normal phase using two silica BEH HILIC columns in series (150 mm × 2.1 mm; 1.7 µm; Waters Acquity) at a temperature of 25 ^∘C. The mobile phases are hexane (A) and hexane : isopropanol (9:1, v:v) (B). Compounds were isocratically eluted for 25 min with 18 % B, followed by a linear gradient to 35 % B in 25 min and a linear gradient to 100 % B in the 30 min thereafter. The flow rate was kept constant (0.2 mL min⁻¹) during the analysis. The conditions for the APCI-MS were identical to Hopmans et al. (2016). GDGTs were detected in single ion monitoring (SIM) mode of the protonated molecules ($[\mathrm{M}+\mathrm{H}{]}^{+}$) of the various GDGTs. A standard mixture of C₄₆ GDGT (internal standard) and crenarchaeol was analyzed to determine the relative response factor (RFF) between these two compounds (see Huguet et al., 2006) and thereby quantify GDGTs in the sediments.

For reconstruction of past SST we used the TEX${}_{\mathrm{86}}^{\mathrm{H}}$ index as proposed by Kim et al. (2010), which is defined as the logarithmic function of the original TEX₈₆ (Schouten et al., 2002):

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{\mathrm{TEX}}_{\mathrm{86}}^{\mathrm{H}}=\log {\displaystyle \frac{\left[\text{GDGT-2}\right]+\left[\text{GDGT-3}\right]+\left[\mathrm{Cren}{}^{\prime }\right]}{\left\{\begin{array}{c}\left[\text{GDGT-1}\right]+\left[\text{GDGT-2}\right]\\ +\left[\text{GDGT-3}\right]+\left[\mathrm{Cren}{}^{\prime }\right]\end{array}\right\}}},\end{array}$$

where the numbers correspond to the amount of cyclopentane moieties in the isoprenoid GDGTs and where Cren^′ refers to the later eluting isomer of crenarchaeol (Sinninghe Damsté et al., 2002). We have discarded 13 samples for the TEX₈₆ calculation due to the partial co-elution of GDGT-2 with an unknown compound.

The TEX${}_{\mathrm{86}}^{\mathrm{H}}$ values were converted to SSTs by applying the global core-top calibration of Kim et al. (2010):

$$\begin{array}{}\text{(2)}& \mathrm{SST}=\mathrm{68.4}\times {\mathrm{TEX}}_{\mathrm{86}}^{\mathrm{H}}+\mathrm{38.6}.\end{array}$$

There is also a regional TEX${}_{\mathrm{86}}^{\mathrm{H}}$ calibration available based on Chilean surface sediments between 25 and 50^∘ S (Kaiser et al., 2015):

$$\begin{array}{}\text{(3)}& \mathrm{SST}=\mathrm{59.6}\times {\mathrm{TEX}}_{\mathrm{86}}^{\mathrm{H}}+\mathrm{33.0}.\end{array}$$

This calibration has a similar slope compared to the global core-top calibration of Kim et al. (2010) but a lower intercept, which results in SSTs that were ca. 2.4 ^∘C lower compared to the TEX${}_{\mathrm{86}}^{\mathrm{H}}$ SSTs calculated after Kim et al. (2010). Although the outcomes are relatively similar, the calibration of Kim et al. (2010) resulted in SSTs that agreed better with the U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ and LDI records, and therefore we have used this calibration. In Appendix A, the TEX${}_{\mathrm{86}}^{\mathrm{H}}$ temperatures based on the calibration of Kaiser et al. (2015) are plotted (Fig. A1 in Appendix). Additionally, we plotted TEX₈₆ temperatures calculated after the Bayesian calibration of Tierney and Tingley (2014, 2015) in Fig. A1 in the Appendix. These SSTs are on average 4 ^∘C lower compared to the TEX${}_{\mathrm{86}}^{\mathrm{H}}$-derived temperatures after Kim et al. (2010), which might be due to the relatively large number of high-latitude core tops on which the BAYSPAR calibration is based (Tierney and Tingley, 2014, 2015), whereas the TEX${}_{\mathrm{86}}^{\mathrm{H}}$ calibration of Kim et al. (2010) excludes (sub)polar core-top data. However, both reconstructions show the same trend.

To assess continental organic matter input into the marine realm, we calculated the branched isoprenoid tetraether (BIT) index, as proposed by Hopmans et al. (2004), including the 6-methyl brGDGTs as described by de Jonge et al. (2014, 2015):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\mathrm{BIT}={\displaystyle \frac{\left\{\begin{array}{c}\left[\mathrm{brGDGT}\mathrm{Ia}\right]+\left[\mathrm{brGDGT}\mathrm{IIa}+{\mathrm{IIa}}^{\prime }\right]\\ +\left[\mathrm{brGDGT}\mathrm{IIIa}+{\mathrm{IIIa}}^{\prime }\right]\end{array}\right\}}{\left\{\begin{array}{c}\left[\mathrm{brGDGT}\mathrm{Ia}\right]+\left[\mathrm{brGDGT}\mathrm{IIa}+{\mathrm{IIa}}^{\prime }\right]\\ +\left[\mathrm{brGDGT}\mathrm{IIIa}+{\mathrm{IIIa}}^{\prime }\right]+\left[\mathrm{Cren}\right]\end{array}\right\}}},\end{array}$$

where the numbers correspond to different branched GDGTs (Hopmans et al., 2004).

Additionally, we calculated the methane index (MI) as a proxy for the dissociation of marine gas hydrates conducive to the anaerobic oxidation of methane (AOM), as proposed by Zhang et al. (2011):

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\mathrm{MI}={\displaystyle \frac{\left[\text{GDGT-1}\right]+\left[\text{GDGT-2}\right]+\left[\text{GDGT-3}\right]}{\left\{\begin{array}{c}\left[\text{GDGT-1}\right]+\left[\text{GDGT-2}\right]+\left[\text{GDGT-3}\right]\\ +\left[\mathrm{Cren}\right]+\left[{\mathrm{Cren}}^{\prime }\right]\end{array}\right\}}}.\end{array}$$

The archaea living in hydrate-impacted environments mainly contain GDGT-1, GDGT-2, and GDGT-3 as their membrane lipids. Therefore, relatively high abundances of these GDGTs might indicate these types of environments in the past and potentially explain erroneous TEX${}_{\mathrm{86}}^{\mathrm{H}}$ results due to high GDGT-1, GDGT-2, and/or GDGT-3 abundances. For all samples, this index was < 0.3, except for one data point (MI = 0.6), which was therefore removed from further discussion. Additionally, to assess other potential influence on the TEX${}_{\mathrm{86}}^{\mathrm{H}}$, we determined the %GDGT-0 and ring index for all sediments. All values were below the advised thresholds, implying no substantial biases on the TEX₈₆ (Zhang et al., 2006, 2011; Sinninghe Damsté et al., 2012).

2.5.2 Long-chain alkenones

Ketone fractions were dissolved in ethyl acetate to a concentration of ∼1 mg mL⁻¹ and analyzed on an Agilent 6890N gas chromatograph (GC) with flame ionization detection (FID). Separation was achieved on a fused silica column with a length of 50 m and diameter of 0.32 mm, coated with a CP Sil-5 (thickness = 0.12 µm). Helium was used as a carrier gas. The flow mode was a constant pressure of 100 kPa. Alkenones were injected on-column at 70 ^∘C at the start of the analysis and increased by 20 ^∘C min⁻¹ to 200 ^∘C, followed by 3 ^∘C min⁻¹ until the final temperature of 320 ^∘C. This end temperature was held for 25 min. Quantification of the alkenones was achieved by means of the C_19:0 ketone internal standard. Identification of the long-chain alkenones was done on an Agilent 7890B GC system interfaced with an Agilent 5977A MS. Separation was achieved on a CP Sil-5 column with an identical diameter and film thickness as that of the GC-FID, but a length of 25 m. Helium was the carrier gas, maintaining a constant flow rate of 2 mL min⁻¹. The MS operated at 70 eV. For both systems, the injection volume was 1 µL. The long-chain alkenones were identified in full scan, scanning between m∕z 50 and 850, and comparison with the literature (de Leeuw et al., 1980; Volkman et al., 1980; Marlowe et al., 1984).

The U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ index was calculated according to Prahl and Wakeham (1987):

$$\begin{array}{}\text{(6)}& {\mathrm{U}}^{{\mathrm{K}}^{\prime }}\mathrm{37}={\displaystyle \frac{\left[{\mathrm{C}}_{\mathrm{37}:\mathrm{2}}\right]}{\left[{\mathrm{C}}_{\mathrm{37}:\mathrm{2}}\right]+\left[{\mathrm{C}}_{\mathrm{37}:\mathrm{3}}\right]}}.\end{array}$$

The U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ values were converted to SSTs using the calibration of Müller et al. (1998):

$$\begin{array}{}\text{(7)}& \mathrm{SST}={\displaystyle \frac{{\mathrm{U}}^{{\mathrm{K}}^{\prime }}\mathrm{37}-\mathrm{0.044}}{\mathrm{0.033}}}.\end{array}$$

There is also a Bayesian calibration available for the U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$, called the BAYSPLINE (Tierney and Tingley, 2018), but below ∼24 ^∘C temperature estimates are similar to those of the calibration of Müller et al. (1998). Since our temperatures are well below 24 ^∘C, we have applied the calibration of Müller et al. (1998).

2.5.3 LCDs

The LCDs were silylated prior to analysis. Polar fractions were dissolved in 25 µL of pyridine and 25 µL of N,O-bis(trimethylsilyl)trifluoroacetamide (BSTFA) and heated at 60 ^∘C for 20 min. Prior to injection, 450 µL ethyl acetate was added. GC-MS analysis was carried out on an Agilent 7890B gas chromatograph coupled to an Agilent 5977A mass spectrometer. Samples were injected at 70 ^∘C. The oven temperature was programmed to 130 ^∘C by 20 ^∘C min⁻¹ and subsequently to 320 ^∘C by 4 ^∘C min⁻¹; this final temperature was held for 25 min. The GC was equipped with an on-column injector and fused silica column (25 m × 0.32 mm) coated with CP Sil-5 (film thickness 0.12 µm). The carrier gas was helium at a constant flow of 2 mL min⁻¹. The mass spectrometer operated with an ionization energy of 70 eV. The injection volume was 1 µL. Identification of the LCDs was achieved in full scan, scanning between m∕z 50 to 850, and on the basis of their characteristic fragmentation (Versteegh et al., 1997). Quantification of the LCDs and proxy computation was done by analysis in SIM mode of the characteristic fragments (m∕z 299, 313, 327, and 341; Rampen et al., 2012; m∕z 187 for C₂₂ 7,16-diol internal standard). We applied the following correction factors for the relative contribution of the selected fragment ions during SIM to the total ion counts: 16 % for all saturated diols, 9 % for all unsaturated diols, and 25 % for the C₂₂ 7,16-diol internal standard.

Past SST was reconstructed by means of the LDI (Rampen et al., 2012):

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }\mathrm{LDI}={\displaystyle \frac{\left[{\mathrm{C}}_{\mathrm{30}}\mathrm{1},\mathrm{15}-\mathrm{diol}\right]}{\left\{\begin{array}{c}\left[{\mathrm{C}}_{\mathrm{28}}\mathrm{1},\mathrm{13}-\mathrm{diol}\right]\\ +\left[{\mathrm{C}}_{\mathrm{30}}\mathrm{1},\mathrm{13}-\mathrm{diol}\right]+\left[{\mathrm{C}}_{\mathrm{30}}\mathrm{1},\mathrm{15}-\mathrm{diol}\right]\end{array}\right\}}}.\end{array}$$

Subsequently, LDI values were converted to SST values via the following equation:

$$\begin{array}{}\text{(9)}& \mathrm{SST}={\displaystyle \frac{\mathrm{LDI}-\mathrm{0.095}}{\mathrm{0.033}}}.\end{array}$$

Three samples were discarded for the calculation of the LDI due to partial co-elution of the C₂₈ 1,13-diol with another unknown compound.

For the reconstruction of past upwelling conditions, the diol index as proposed by Willmott et al. (2010) was applied:

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }\mathrm{diol}\mathrm{index}={\displaystyle \frac{\left[{\mathrm{C}}_{\mathrm{28}}\mathrm{1},\mathrm{14}-\mathrm{diol}\right]+\left[{\mathrm{C}}_{\mathrm{30}}\mathrm{1},\mathrm{14}-\mathrm{diol}\right]}{\left\{\begin{array}{c}\left[{\mathrm{C}}_{\mathrm{28}}\mathrm{1},\mathrm{14}-\mathrm{diol}\right]+\left[{\mathrm{C}}_{\mathrm{30}}\mathrm{1},\mathrm{14}-\mathrm{diol}\right]\\ +\left[{\mathrm{C}}_{\mathrm{28}}\mathrm{1},\mathrm{13}-\mathrm{diol}\right]+\left[{\mathrm{C}}_{\mathrm{30}}\mathrm{1},\mathrm{13}-\mathrm{diol}\right]\end{array}\right\}}}.\end{array}$$

The diol index according to Rampen et al. (2008) was shown to also be affected by variations in SST due to the inclusion of the C₃₀ 1,15-diol in the ratio (Rampen et al., 2014b; Zhu et al., 2018), and since on a glacial–interglacial timescale we expect substantial SST differences, we consider the Willmott et al. (2010) ratio more appropriate.

Possible fluvial input of LCDs was assessed by the fractional abundance of the C₃₂ 1,15-diol, which is potentially derived from rivers (de Bar et al., 2016; Lattaud et al., 2017a, b):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }F{\mathrm{C}}_{\mathrm{32}}\mathrm{1},\mathrm{15}-\mathrm{diol}\\ \text{(11)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\left[{\mathrm{C}}_{\mathrm{32}}\mathrm{1},\mathrm{15}-\mathrm{diol}\right]}{\left\{\begin{array}{c}\left[{\mathrm{C}}_{\mathrm{28}}\mathrm{1},\mathrm{13}-\mathrm{diol}\right]+\left[{\mathrm{C}}_{\mathrm{30}}\mathrm{1},\mathrm{13}-\mathrm{diol}\right]\\ +\left[{\mathrm{C}}_{\mathrm{30}}\mathrm{1},\mathrm{15}-\mathrm{diol}\right]+\left[{\mathrm{C}}_{\mathrm{32}}\mathrm{1},\mathrm{15}-\mathrm{diol}\right]\end{array}\right\}}}.\end{array}$$

The NDI, a proposed proxy for ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations, was calculated following Gal et al. (2018):

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }\mathrm{NDI}={\displaystyle \frac{\left[{\mathrm{C}}_{\mathrm{28}}\mathrm{1},\mathrm{14}-\mathrm{diol}\right]+\left[{\mathrm{C}}_{\mathrm{28}:\mathrm{1}}\mathrm{1},\mathrm{14}-\mathrm{diol}\right]}{\left\{\begin{array}{c}\left[{\mathrm{C}}_{\mathrm{28}}\mathrm{1},\mathrm{14}-\mathrm{diol}\right]+\left[{\mathrm{C}}_{\mathrm{28}:\mathrm{1}}\mathrm{1},\mathrm{14}-\mathrm{diol}\right]\\ +\left[{\mathrm{C}}_{\mathrm{30}}\mathrm{1},\mathrm{14}-\mathrm{diol}\right]+\left[{\mathrm{C}}_{\mathrm{30}:\mathrm{1}}\mathrm{1},\mathrm{14}-\mathrm{diol}\right]\\ +\left[{\mathrm{C}}_{\mathrm{28}}\mathrm{1},\mathrm{13}-\mathrm{diol}\right]+\left[{\mathrm{C}}_{\mathrm{30}}\mathrm{1},\mathrm{13}-\mathrm{diol}\right]\\ +\left[{\mathrm{C}}_{\mathrm{30}}\mathrm{1},\mathrm{15}-\mathrm{diol}\right]\end{array}\right\}}}.\end{array}$$

The NDI was then translated to [${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$] and [${\mathrm{NO}}_{\mathrm{3}}^{-}$] concentrations by using the following equations (Gal et al., 2018).

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]={\displaystyle \frac{\mathrm{NDI}-\mathrm{0.015}}{\mathrm{0.413}}}\text{(14)}& {\displaystyle }& {\displaystyle }\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]={\displaystyle \frac{\mathrm{NDI}-\mathrm{0.075}}{\mathrm{0.026}}}\end{array}$$

Finally, loliolide, an indicator of diatom abundance (Klok et al., 1984; Repeta, 1989), and dinosterol, a biomarker for dinoflagellates (Boon et al., 1979; Volkman et al., 1998), were identified by GC-MS analysis (simultaneously with the LCDs) of the silylated polar fraction and quantified using the C₂₂ 7,16-diol standard in full scan, correcting for the molecular weights of the compounds.

3 Results

3.1 Bulk parameters and sedimentation rates

The TOC content varies between 0.4 % and 2.6 %. The TOC content is significantly higher during the interglacial periods (MIS 1, 3, and 5) compared to glacial periods (MIS 2, 4, and 6; two-tailed p< 0.001): the average TOC concentration is 1.4 % for the interglacial periods and 0.7 % for the glacial intervals (Fig. 3d), although TOC levels for MIS 3 are quite similar to those of MIS 1 and 2. Similarly, the average TN levels are 0.2 % and 0.1 % during the interglacial and glacial intervals, respectively (Fig. 3d). During Termination 2, both the TOC and TN contents increase rapidly (within < 2 kyr) towards interglacial values, and the highest TOC and TN values are observed around 78 ka (2.6 % and 0.3 %, respectively). The atomic C∕N ratio (Fig. 3b) varies from 4.3 to 10.5 and is on average higher during the interglacial periods (7.2) compared to glacial times (6.2). The organic matter δ¹³C record (δ¹³C_OM) also reveals a glacial–interglacial variation (Fig. 3c) corresponding to slightly ¹³C-enriched values during interglacial times (${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{average}}=-\mathrm{21.2}$ ‰) compared to the glacials (${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{average}}=-\mathrm{21.8}$ ‰). Although small, these changes are statistically significant (5 % significance level, two-tailed p< 0.001).

Figure 3 Geochemical bulk records for the studied interval of ODP 1234. (a) Sedimentation rate (blue) and MAR_TOC (orange), (b) atomic C∕N ratio, (c) bulk organic δ¹³C, (d) TOC (red) and TN (green) concentrations. The different color bands indicate different time periods: interglacial stages MIS 1, 3, and 5 in pink, including the MIS 5e in dark pink. The LGM is highlighted in blue. MIS ages are according to Lisiecki and Raymo (2005). The dashed lines indicate the onset ages of Terminations 1 (Pacific age of Stern and Lisiecki, 2014) and 2 (approximate age).

Using the age model modified from Heusser et al. (2006), we estimated sedimentation rates that varied between ca. 0.2 and 1.7 mm yr⁻¹ (Fig. 3a). Sedimentation rates were highest in MIS 4, during which the sedimentation rate reaches values of around 2 mm yr⁻¹. Sedimentation rates were lowest during the warmest periods (MIS 1 and 5e; ca. 0.3–0.4 mm yr⁻¹). The MAR_TOC varied between ca. 2 and 32 g m⁻² yr⁻¹ and shows a relatively similar pattern as the sedimentation rate, implying that the sedimentation rate strongly controls the AR of TOC (Fig. 3a). Around ca. 100 ka there is a pronounced maximum in the TOC AR, reaching ca. 32 g m⁻² yr⁻¹.

3.2 Lipid biomarker concentrations and ARs

In all sediments, long-chain alkenones, GDGTs, and LCDs were present above the quantification limit. The total concentration of crenarchaeol, a specific biomarker for marine Thaumarchaeota (Schouten et al., 2002, 2003), varies between ca. 2 and 16 µg g⁻¹, with highest values during MIS 5, after MIS 5e (Eemian Interglacial; Fig. 4a). Accumulation rates (ARs) range from 1 to 25 mg m⁻² yr⁻¹ (Fig. 4a). The summed concentration of the di- and tri-unsaturated C₃₇ alkenones, a productivity proxy for alkenone-producing haptophytes (e.g., Prahl et al., 1988; Prahl and Muehlhausen, 1989; Rostek et al., 1997; Villanueva et al., 1997), varies between ca. 0.5 and 10 µg g⁻¹. Highest concentrations are observed for MIS 5e (126–116 ka; Fig. 4b). Around the boundary of MIS 4 and 5, concentrations decrease from ca. 6 to 2 µg g⁻¹ sediment; however, the alkenone AR does not show this decrease. During the Late Holocene, C₃₇ alkenone concentrations increase again. The alkenone AR record resembles that of the MAR_TOC with a peak around 100 ka of ca. 10 mg m⁻² yr⁻¹ (Fig. 4b). Dinosterol shows the highest abundance during the last interglacial (MIS 5), peaking just after the Eemian Interglacial (between ca. 116 and 110 ka; Fig. 4c) with concentrations of ca. 4 µg g⁻¹. This peak is followed by a gradual decrease towards the LGM with values near 1 µg g⁻¹ sediment. The AR record of dinosterol highly resembles the concentration record but peaks somewhat later, i.e., between ca. 105 and 98 ka. The loliolide concentration shows one maximum during the Late Holocene of ca. 2.3 µg g⁻¹ sediment (Fig. 4d). Its AR varies between ca. 0 and 1 mg m⁻² yr⁻¹ throughout the record, with one pronounced peak around 100 ka of around 2 mg m⁻² yr⁻¹.

Figure 4 Biomarker proxy and accumulation records for ODP 1234. (a) Crenarchaeol concentrations and MARs. (b) C₃₇ alkenone concentrations and MARs. (c) Dinosterol concentrations and MARs. (d) Loliolide concentrations and MARs. (e) 1,13- and 1,15-diol MARs. (f) 1,14-diol concentration and MARs. (g) Diol index (blue) and biogenic opal concentrations (green; Muratli et al., 2010b) for the last ∼28 kyr. (h) Phosphate and nitrate concentrations calculated after the NDI. (i) Total diatom counts (purple) and Chaetoceros diatom counts (blue) (Mix et al., 2003). The thick smoothed lines reflect 2-point running averages of the diatom records. The counts are descriptive. F: few; C: common; A: abundant; VA: very abundant. (j) Mass AR of TOC. For panels (a)–(d) and (f): concentrations are in grey and MARs are in black. The different color bands indicate different time periods: interglacial stages MIS 1, 3, and 5 in pink, including the MIS 5e in dark pink. The LGM is highlighted in blue. The purple band indicates the period of presumed enhanced upwelling. MIS ages are according to Lisiecki and Raymo (2005). The dashed lines indicate the onset ages of Terminations 1 (Pacific age of Stern and Lisiecki, 2014) and 2 (approximate age).

The main LCDs detected throughout the core are the C₂₈ and C₃₀ 1,13-diols, the C₂₈, C₃₀, and C_30:1 1,14-diols, and the C₃₀ and C₃₂ 1,15-diols. We did not detect the C₃₂ 1,14-diol (characteristic for Apedinella radians), but we did identify the C_30:1 1,14-diol and C₂₉ 12-OH fatty acid, which are typical biomarkers for Proboscia diatoms (Sinninghe Damsté et al., 2003), implying that Proboscia diatoms are most likely the source of the 1,14-diols detected in the sediments. The profile of the summed concentration of 1,14-diols shows four distinct sharp peaks during the MIS 5 (Fig 4f), reaching concentrations up to 1.2 µg g⁻¹. The 1,14-diol AR shows one distinct peak at 94 ka of ca. 1.7 mg m⁻² yr⁻¹, whereas throughout the rest of the record the AR varies between 0 and 0.6 mg m⁻² yr⁻¹. Figure 4h shows the ARs of the summed 1,13- and 1,15-diols. The AR of the 1,13-diols ranges between ca. 0.1 and 0.8 mg m⁻² yr⁻¹ throughout MIS 6 and 5 and peaks (ca. 1.2 mg m⁻² yr⁻¹) at the end of MIS 5, revealing the highest ARs during MIS 4, followed by a subsequent gradual decrease towards Holocene values of around 0.2 mg m⁻² yr⁻¹. The AR record of the 1,15-diols is highly similar to that of the 1,13-diols for MIS 1 to 4, as well as for MIS 6, but shows higher values during the second half of MIS 5 with maxima around 94, 104, and 122 ka of between ca. 1 and 1.3 mg m⁻² yr⁻¹.

Figure 5 Foraminiferal oxygen isotope and organic temperature proxy records for ODP 1234. (a) The benthic stable oxygen isotope records for ODP Site 1234 of Heusser et al. (2006), in which the blue data reflect the benthic δ¹⁸O dated by correlation with Atlantic core MD95204, and the orange data represent the δ¹⁸O data correlated with the Vostok ice core chronology. (b) The organic geochemical seawater surface temperature records: U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ (green), TEX${}_{\mathrm{86}}^{\mathrm{H}}$ (red), and LDI (blue). The thick lines reflect 3-point running averages. (c) Planktonic foraminiferal δ¹⁸O of the species Globigerina bulloides with the 3-point running average. (d) Global compilation of benthic stable oxygen isotopes (Lisiecki and Raymo, 2005). The pink colored bands reflect the interglacial periods (MIS 1, 3, and 5). The pink lines within MIS 5 reflect the different substage δ¹⁸O minima and maxima. MIS (substage) ages are according to Lisiecki and Raymo (2005). The blue band reflects the Last Glacial Maximum, and the star symbol represents the present-day SST (World Ocean Atlas 2013 version 2). The dashed lines indicate the onset ages of Terminations 1 (Pacific age of Stern and Lisiecki, 2014) and 2 (approximate age).

3.3 Foraminiferal stable isotopic composition

We generated a δ¹⁸O record of G. bulloides for MIS 5 from ca. 70 to 142 ka, as well as for a short interval within MIS 3 from ca. 36 to 41 ka, which is plotted in Fig. 5b. The δ¹⁸O values of G. bulloides varied between ca. 0.3 ‰ and 4.2 ‰ and are most depleted during the Eemian Interglacial around 123 ka (δ¹⁸O = 0.3 ‰); enriched values are observed for MIS 6 and 3.

3.4 Organic proxy records

The diol index record shows several maxima, the most evident of which is during MIS 5 (Fig. 4g). During MIS 5e, the diol index reveals two peaks around 126 and 116 ka (up to 0.7), followed by a sharp drop and a subsequent increase around 110 ka, reaching a maximum value at ∼86 ka (up to 0.8). After 86 ka, the index gradually decreases. During MIS 2–4 the diol index varies between ca. 0.1 and 0.5. During MIS 1, the diol index shows a gradual increase from ca. 0.15 to 0.35.

We did not detect the monounsaturated C₂₈ 1,14-diol (C_28:1 1,14-diol), and therefore this diol could not be included in the calculation of the NDI. In ∼64 % of the surface sediments of the datasets of Rampen et al. (2014b) and de Bar et al. (2016), on which the NDI calibration of Gal et al. (2018) is based, the C_28:1 1,14 was also not detected. The NDI record reveals similar variations as the diol index, although more pronounced. During MIS 6 the NDI is close to zero, followed by a small peak at the end of the Eemian Interglacial (NDI = 0.4), coincident with the peak in the diol index at 116 ka. NDI-derived [${\mathrm{NO}}_{\mathrm{3}}^{-}$] and [${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$] concentrations around 116 ka are 11 and 0.8 µmol L⁻¹, respectively (Fig. 4h). Similar to the diol index, the ratio then decreases after which the NDI shows a broad peak (∼100–80 ka) with a maximum at 88 ka (NDI = 0.6), which gives concentrations of 20 and 1.4 µmol L⁻¹ when translated to [${\mathrm{NO}}_{\mathrm{3}}^{-}$] and [PO${}_{\mathrm{4}}^{\mathrm{3}-}$], respectively (Fig. 4h). Between 80 and 0 ka (i.e., MIS 4–1), the NDI does not reveal any distinct maxima, and the ratio varies between ca. 0 and 0.2.

Figure 6 Cross-correlations of sea surface temperature estimates based on the TEX${}_{\mathrm{86}}^{\mathrm{H}}$ index, U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ index, and LDI. Linear regressions are indicated by the black solid lines, together with the 95 % confidence intervals. All correlations are significant (p< 0.001). The orange dashed line represents the 1:1 line.

Overall the organic SST proxy records broadly follow the same trend, clearly revealing glacial–interglacial temperature variability (Fig. 5b). Indeed, when we cross-correlate the three proxies, all combinations display significant positive correlations (p< 0.001; Fig. 6), albeit with some scatter.

Additionally, the SST trends show good correspondence with the planktonic δ¹⁸O record of G. bulloides (Fig. 5c), as well as with the benthic oxygen isotope record for ODP Site 1234 (Heusser et al., 2006; Fig. 5a) and the global stack of benthic δ¹⁸O values (e.g., Lisiecki and Raymo, 2005; Fig. 5d). Termination 2 (around 130 ka) is clearly expressed in all three organic proxy SST records, showing a rise in temperature of approximately 4 ^∘C within ca. 2 kyr. During MIS 5e, SSTs were between ca. 16 and 18 ^∘C based on the three temperature proxies. Then, during MIS 5, the interglacial substages representing alternating cold and warm periods are clearly reflected in all three proxy records. The three records show temperature drops between ca. 4 and 6 ^∘C during MIS 5a (ca. 70 ka). During MIS 4, 3, and 2, all SST records show a gradual decrease towards the LGM, during which U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ and TEX${}_{\mathrm{86}}^{\mathrm{H}}$ reveal temperatures between 9 and 10 ^∘C and an LDI between 7 and 8 ^∘C. Around 25 ka, U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ and TEX${}_{\mathrm{86}}^{\mathrm{H}}$ temperatures show a similar steady rise to Holocene temperatures of around 16 and 18 ^∘C. The LDI reveals an approximate 7 ^∘C warming over Termination 1, starting slightly later at around 22 ka. During the last ∼ 20 kyr, LDI-derived SSTs are ∼3–4 ^∘C lower compared to the U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$- and TEX${}_{\mathrm{86}}^{\mathrm{H}}$-derived SSTs.

4 Discussion

4.1 Productivity

The average bulk organic carbon δ¹³C over the last 150 kyr is −21.4 (±0.6, SD) ‰, and the atomic C∕N ratio displays an average value of 6.8 (±1.3), both implying a dominant marine organic carbon source (Bordovskiy, 1965; Emerson and Hedges, 1988; Meyers, 1997). Today, the amount of sedimentary organic carbon is a good indicator of export production along the Chilean margin (Hebbeln et al., 2000). Hence, higher TOC levels during the interglacial intervals in the ODP 1234 record may suggest enhanced marine productivity during these intervals for this site. However, enhanced interglacial productivity is in contrast with previous studies that state that coastal productivity on the Peru margin was highest during the LGM and diminished during the Holocene (e.g., Thomas et al., 1994; Marchant et al., 1999; Thomas, 1999; Hebbeln et al., 2000, 2002; Lamy et al., 2002, 2004; Romero and Hebbeln, 2003; Mohtadi and Hebbeln, 2004). These studies suggest that during the LGM the ACC migrated northward, supplying nutrients (particularly nitrate and phosphate; Levitus et al., 1994), together with a northward shift of the southern westerly belt as the main precipitation source onshore, resulting in enhanced micronutrient supply via continental runoff. This combined effect would have stimulated productivity along the Chilean coast during the LGM. Upon deglaciation the climate zones propagated southward, resulting in a lowering of productivity, as indicated by pollen, sedimentological, and continental studies (e.g., Heusser, 1990; Lamy et al., 1998, 1999, 2001, 2004; Brathauer and Abelmann, 1999; Haberle and Bennett, 2004; Stuut and Lamy, 2004; Kaiser et al., 2005; Heusser et al., 2006). These studies suggest a 5 to 6^∘ northward movement of the southern westerly belt during the LGM.

In contrast to the hypothesis of stimulated productivity during the LGM and lower productivity during the Holocene, two studies suggested low paleoproductivity during the LGM for core sites at the same latitude as that of our core (35–36^∘ S) (Romero et al., 2006; Mohtadi et al., 2008). Mohtadi et al. (2008) hypothesized that if the climate zones shifted 5 to 6^∘ northward during the LGM, the SWW would be just above the core site (35–36^∘ S), blowing directly onshore and thereby preventing coastal upwelling. In turn, during the Early Holocene the subtropical high pressure would have become the dominant atmospheric player according to Romero et al. (2006), favoring upwelling. However, Muratli et al. (2010b), who reconstructed paleoproductivity over the last ∼30 kyr for ODP 1234, suggested that the rise in TOC after the LGM is probably not the result of increased productivity, but of a lower oxygen availability due to decreased Antarctic Intermediate Water ventilation and thus increased preservation. The authors based these suggestions on the fact that neither opal concentrations nor opal or TOC MARs increased simultaneously with TOC concentrations. Chase et al. (2014) support this hypothesis based on the Th-normalized organic carbon fluxes. Our (unnormalized) MAR_TOC supports these findings as well since although TOC levels show a steady rise after the LGM, we do not observe this increase in the accumulation of organic carbon, and the MAR_TOC is higher during the LGM compared to the Holocene, suggesting higher glacial productivity. In summary, our findings seem consistent with previous studies suggesting that higher TOC levels during the Holocene likely do not result from increased productivity, but rather from increased preservation.

Interestingly, during MIS 3, 4, and 5 we observe even higher MAR_TOC values, reaching a pronounced maximum around 100 ka during MIS 5, i.e., during an interglacial period. Hence, this record suggests that the general productivity over the last ∼150 kyr was highest around 100 ka. This interpretation is supported by most of the AR records of the individual lipid biomarkers peaking around 100 ka. The C₃₇ alkenone, crenarchaeol, and dinosterol ARs are all at their maximum between ca. 105 and 95 ka (Fig. 4b, a, and c, respectively), simultaneous with maximal MAR_TOC (Fig. 4j). The AR of loliolide (indicating diatom abundance) peaks at 100 ka (Fig. 4d), coincident with a maximum in total diatoms and Chaetoceros diatom counts (Mix et al., 2003; Fig. 4c) and concurrent with the peak in MAR_TOC. Chaetoceros diatoms are generally associated with upwelling (e.g., Abrantes, 1988; Abrantes and Moita, 1999), and also for the Peru–Chile upwelling system, the diatom genus Chaetoceros is thought to dominate the diatom community during upwelling conditions (e.g., Anabalón et al., 2007; Schrader and Sorknes, 1991; Romero et al., 2001; Vargas et al., 2004; Abrantes et al., 2007; González et al., 2007; Sanchez et al., 2012). Therefore, the high abundance of Chaetoceros diatoms at ∼100 ka suggests that there was a peak in upwelling intensity around this time. This upwelling would have introduced macronutrients from deeper colder waters into the euphotic zone, stimulating several phytoplankton communities, including haptophytes, dinoflagellates, and diatoms. During MIS 4 we also observe a peak in MAR_TOC (Fig. 3a) but there is no corresponding peak in Chaetoceros abundance (Fig. 4i). However, the high MAR_TOC during MIS 4 is linked to high sedimentation rates (Fig. 3c), whereas the peak in MAR_TOC during MIS 5 is not, confirming that the MAR_TOC maximum during MIS 5 indeed likely resulted from increased primary productivity, whereas during MIS 4 the high MAR_TOC did not.

4.2 1,14-diols as past upwelling and nutrient indicators

Interestingly, peaks in the diol index and 1,14-diol concentrations and accumulation rates occur before and after (around 116 and 88 ka; Fig. 4a and j), but not during, the time interval with enhanced upwelling (i.e., around 100 ka; purple band in Fig. 3). Contreras et al. (2010) also reported high abundances of the C₂₈ 1,14-diol around 120 ka near the coast of Peru (11^∘ S). This suggests that Proboscia diatoms were more abundant before and after this period of intense upwelling and is in agreement with observations that in the present-day Chile–Peru upwelling region, Proboscia alata is more dominant when upwelling is less intense (Tarazona et al., 2003; Herrera and Escribano, 2006). Although the recent timescale is quite different from our long time record, which represents an integrated signal of several hundreds to thousands of years, it suggests that over the time period 120–80 kyr, on average, upwelling became stronger, reaching a maximum at 100 ka before subsequently decreasing. Shortly before and after the period of maximum upwelling, the conditions (averaged over multiple years) were apparently optimal for Proboscia diatoms. As previously suggested by Rampen et al. (2014b), this indicates that the diol index should perhaps be considered as a specific indicator for Proboscia productivity, rather than upwelling strength generally, as the environmental conditions determining Proboscia abundance likely differ from region to region.

Recently, the NDI was introduced as a quantitative paleo-nutrient proxy. The diol index, the 1,14-diol AR, and the NDI are low during the time of most intense upwelling (around 100 ka), suggesting a minimum in annual mean nitrate and phosphate concentrations, which is highly unlikely. The NDI is based on the saturated and monounsaturated C₂₈ 1,14-diol relative to other diols. Whereas we observe a strong correlation between the ARs of the monounsaturated and saturated C₃₀ 1,14-diol (R²=0.70), the correlations between the C₂₈ 1,14 and C₃₀ 1,14-diol (R²=0.32) and the C₂₈ 1,14 and C_30:1 1,14-diol (R²=0.36) are weak, indicating different source organisms for the C₂₈ and C₃₀ 1,14-diols, in agreement with previous studies (Rampen et al., 2014b; de Bar et al., 2016; Gal et al., 2018). In fact, Sinninghe Damsté et al. (2003) found that 98 % of the diols in P. alata consisted of the saturated and monounsaturated C₂₈ 1,14-diol, whereas P. indica produces similar amounts of the C₃₀ and C₂₈ 1,14-diol and P. inermis mainly produced the C₂₈ 1,14-diol (Rampen et al., 2007). P. alata is the Proboscia species detected along the modern coast of Chile and Peru (Tarazona et al., 2003; Herrera and Escribano, 2006), and therefore it is likely that P. alata is the predominant producer of the C₂₈ 1,1,4-diol detected in the ODP 1234 sediments. Consequently, the NDI record likely reflects variations in the abundance of P. alata over the last 150 kyr in the Chilean margin, whereas the diol index (which also includes the C₃₀ 1,14-diol) more likely reflects the abundance of multiple species of Proboscia.

Figure 7 (a) Pacific U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ SST reconstructions off South America for the last ∼160 kyr (3-point running averages). (b) Sea surface temperature (^∘C) map for the region. The colors of the records correspond to the colors of the site symbols in the map. Alkenone-based SST data is from Caniupán et al. (2011) (MD07-3128; 53^∘ S), Ho et al. (2012) (GeoB 3327-5; 43^∘ S), Kaiser et al. (2005) (ODP1233; 41^∘ S), this study (ODP 1234; 36^∘ S), Kim et al. (2002) (GeoB 3302-1 and GIK 17748-2; 33^∘ S), Bin Shaari (2013) (ODP 1237; 16^∘ S), Bin Shaari et al. (2013) (ODP 1239; 1^∘ S) and Bin Shaari et al. (2014) (ODP 1241; 6^∘ N). Surface temperature data derived from the World Ocean Atlas 2013 (Locarnini et al., 2013); map is drawn in Ocean Data View (Schlitzer, 2015), and modified manually.

Previous studies showed that P. alata preferentially blooms over other diatoms when nutrients are high but silica concentrations are low (as P. alata needs little Si to build its frustule) (Haake et al., 1993; Sakka et al., 1999; Koning et al., 2001; Smith, 2001). In fact, although P. alata is often observed in high nutrient and/or upwelling regions (Hernández-Becerril, 1995; Lange et al., 1998; Koning et al., 2001; Smith, 2001), the conditions in which the species is abundant are often described as post-bloom, stratification, early upwelling season and/or the oceanic side of the upwelling front (e.g., Hart, 1942; Takahashi et al., 1994; Katsuki et al., 2003; Moita et al., 2003; Tarazona et al., 2003; Herrera and Escribano, 2006; Sukhanova et al., 2006; see references in Table 1 of Rampen et al., 2014b). Moreover, sediment trap studies from the Arabian Sea showed that the maximum flux of Proboscia lipids was at the start of the upwelling season (Prahl et al., 2000; Wakeham et al., 2002; Sinninghe Damsté et al., 2003; Rampen et al., 2007). When silicate concentrations increase (i.e., during upwelling), P. alata is likely outcompeted by more heavily silicified diatoms (such as Chaetoceros; e.g., Riegman et al., 1996). This suggests that the C₂₈ 1,14-diol likely reflects early or post-upwelling nutrient conditions. Hence, the dip in reconstructed [${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$] and [${\mathrm{NO}}_{\mathrm{3}}^{-}$] concentrations at ∼100 ka is not realistic and likely due to the low abundance of P. alata at this time of intense upwelling. In other regions, the C₂₈ 1,14-diol might be produced by other Proboscia species (e.g., P. inermis) and the NDI might thus not reflect solely P. alata productivity. In summary, we suggest that the NDI at the Chilean margin likely reflects Proboscia productivity and may therefore not be suitable as paleo-nutrient tracer.

4.3 Sea surface temperature evolution

The three organic proxy-based SST records broadly follow the trend of the planktonic oxygen isotope record for MIS 5, as well as the benthic oxygen isotope record for ODP Site 1234 (Heusser et al., 2006; Fig 5), and the global stack of benthic δ¹⁸O (e.g., Lisiecki and Raymo, 2005) for the entire record, suggesting that the evolution of SST at our study location generally follows global climate patterns. The resolution of our record is not high enough to recognize possible millennial-scale SST variations related to climatic events such as the Antarctic Cold Reversal. Nevertheless, the overall SST patterns are similar to other Southern Hemisphere records (e.g., Kaiser et al., 2005; Kaiser and Lamy, 2010; Caniupán et al., 2011; Lopes dos Santos et al., 2013) and Antarctic ice core stable isotope records (e.g., Blunier and Brook, 2001). The three temperature proxies are all significantly positively correlated and present a coherent view of regional climate variability (Fig. 6).

In principle, we would expect similar reconstructed temperatures from the LDI and U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ index, since both proxies are based on biomarkers produced by photosynthetic algae (Prahl and Wakeham, 1987; Rampen et al., 2012). Indeed, Kim et al. (2002) showed that for surface sediments off Chile, U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$-derived temperatures strongly correlated with annual mean temperatures of the sea surface mixed layer. The linear regression for the U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ vs. the LDI is close to the 1:1 line, suggesting the LDI is reflecting SST. Interestingly, the LDI reveals a greater amplitude over the record compared to the U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ and TEX${}_{\mathrm{86}}^{\mathrm{H}}$, with a maximum temperature difference between the coldest (LGM) and warmest (Eemian Interglacial) temperature of ca. 10 ^∘C, whereas this is ca. 7 and 8 ^∘C for the U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ and TEX${}_{\mathrm{86}}^{\mathrm{H}}$, respectively. Other LDI and U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ applications on glacial–interglacial timescales also show somewhat greater amplitudes for the LDI compared to the U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ (Rampen et al., 2012; Lopes dos Santos et al., 2013; Rodrigo-Gámiz et al., 2014; Jonas et al., 2017). This might potentially suggest that the surface sediment calibration of the LDI requires some modification, such as a lower slope. Further analysis of surface sediments should reveal this.

For the TEX₈₆, it has been shown to potentially reflect subsurface rather than surface water temperatures (Huguet et al., 2007; Kim et al., 2010, 2015; Schouten et al., 2013; Chen et al., 2014) due to the production of isoprenoid GDGTs below the surface mixed layer. Overall, the TEX${}_{\mathrm{86}}^{\mathrm{H}}$ record agrees reasonably well with the U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ record for ODP 1234, suggesting that it mainly reflects SST. Also, Kaiser et al. (2015), who established a regional TEX${}_{\mathrm{86}}^{\mathrm{H}}$ calibration, suggested that this proxy mainly reflects SST. However, all cross-correlations still reveal relatively large scatter, reflecting the multiple different constraints on the respective proxies, as each one is affected by different parameters and can reflect different seasonal temperatures. Nevertheless, for the largest part of the record, the absolute temperature differences between the three temperature proxies is smaller than the maximal possible discrepancy that can be explained by the combined calibration errors (calibration errors of U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$, TEX${}_{\mathrm{86}}^{\mathrm{H}}$, and LDI are 1.5, 2.5, and 2.0 ^∘C, respectively; Müller et al., 1998; Kim et al., 2010; Rampen et al., 2012, respectively).

The only period during which the offsets are larger than the calibration errors is the last ∼25 kyr, i.e., from the LGM to the Holocene, and the relatively brief interval between 52 and 56 ka. During these time spans, the LDI-derived temperatures are ca. 3 to 6 ^∘C lower than those derived from U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ or TEX${}_{\mathrm{86}}^{\mathrm{H}}$. This offset could potentially be related to terrestrial input; however, we have assessed the relative contribution of terrestrially derived organic carbon by means of the BIT index (Hopmans et al., 2004) and the fractional abundance of the C₃₂ 1,15-diol (de Bar et al., 2016; Lattaud et al., 2017a, b), which were both always < 0.2, implying that there was no substantial riverine input from terrigenous organic matter. Alternatively, the LDI might be compromised by a high input of 1,14-diols (Rodrigo-Gámiz et al., 2015); however, we do not observe a correlation between the LDI (or the SST offset between the LDI and U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ or TEX${}_{\mathrm{86}}^{\mathrm{H}}$) and the fractional abundance of the summed 1,14-diols. Potentially, the production of the 1,13- and 1,15-diols shifted to a colder season due to regional climatic change associated with the southward migration of the ACC and SWW upon the deglacial warming. Finally, we observe an early onset of the deglacial warming for the LDI compared to the U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ and TEX${}_{\mathrm{86}}^{\mathrm{H}}$ (ca. 6 kyr; Fig. 5b), the reasons for which are elusive.

We can compare our SST records of ODP Site 1234 with those of other SST records along the coast of South America, in particular those generated by the U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$, which has been used the most often (Fig. 7). Our U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ temperature record agrees well with other records in the vicinity of ODP 1234, both in terms of absolute temperature and in glacial–interglacial temperature amplitude. Our U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ SST record shows an approximate 7 ^∘C warming over Termination 1, which is in agreement with other past SST records for the central and southern Chilean margin (Kim et al., 2002; Lamy et al., 2002, 2004; Kaiser et al., 2005; Romero et al., 2006; Kaiser and Lamy, 2010). Although the timing of deglacial warming associated with Termination 2 is comparable to other records, the last deglaciation (Termination 1) reveals an early start (around 26 ka) at Site 1234 compared to neighboring sites. The TEX${}_{\mathrm{86}}^{\mathrm{H}}$ record reveals this same timing, whereas the LDI suggests a timing more comparable to other records in the vicinity (around 20 ka). Regional differences in the timing of deglacial warming have previously been related to the southward shift of the SWW, which directly and indirectly influences local paleoproductivity and upwelling intensity, thereby potentially leaving site-dependent, unique signatures in the SST records (Mohtadi et al., 2008). This may partially explain the discrepancies that we observe between our U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ record and the other records. For instance, site GeoB 3327-5 is located at the northern extent of the ACC and thus largely influenced by its latitudinal movement, which is closely related to the migration of the SWW (Ho et al., 2012), whereas for ODP 1234 this influence might be less. Furthermore, the combined regional records show an increase in the magnitude of U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ temperature variations over glacial–interglacial timescales towards the south (Fig. 7), which agrees with the idea of greater glacial–interglacial temperature differences at high latitudes compared to low latitudes (e.g., CLIMAP, 1976; Mohtadi and Hebbeln, 2004; Mohtadi et al., 2008).

5 Conclusions

We have applied long-chain diols as tracers for past SST (using the LDI), upwelling and nutrient conditions (diol index), and [${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$] and [${\mathrm{NO}}_{\mathrm{3}}^{-}$] (NDI) at ODP 1234, located within the Peru–Chile upwelling system. The LDI agrees with the TEX${}_{\mathrm{86}}^{\mathrm{H}}$ and U${}_{\mathrm{37}}^{{\mathrm{K}}^{\prime }}$ sea surface temperature records, as well as with the planktonic δ¹⁸O, suggesting that the LDI reflects past SSTs. During the last interglacial (MIS 5), increased accumulation of TOC and phytoplankton lipid biomarkers centered around 100 ka indicate enhanced primary productivity. Concurrently, there is a peak abundance in preserved Chaetoceros diatoms, suggesting a peak in upwelling intensity. The diol index peaks before and after this peak in upwelling, agreeing with present-day diatom distributions along the Chile–Peru margin, with Chaetoceros diatoms being dominant during upwelling and Proboscia alata thriving in more stable waters. The NDI (based primarily on the C₂₈ 1,14-diol) shows the same trend as the diol index (based on both the C₂₈ and C₃₀ 1,14-diols), i.e., also showing a dip around 100 ka, suggesting low mean annual nitrate and phosphate concentrations during this upwelling interval, which is not realistic. Likely P. alata is the main C₂₈ 1,14-diol producer at ODP 1234, suggesting that the NDI likely reflects P. alata productivity at the core site and is therefore not suitable as a paleo-nutrient tracer, since the species is generally outcompeted when Si concentrations increase, i.e., during upwelling (potentially explaining the minimum in NDI around 100 ka). Overall, these data suggest that the NDI and the diol index should perhaps be considered indicators for Proboscia productivity, rather than nutrient concentrations or general upwelling conditions.

Appendix A

Figure A1 TEX₈₆ records for ODP 1234 with different calibrations, i.e., the global core-top calibration of Kim et al. (2010) in red (TEX${}_{\mathrm{86}}^{\mathrm{H}}$), the local core-top calibration of Kaiser et al. (2015) in purple, and the Bayesian model of Tierney and Tingley (2014, 2015) in green applying the BAYSPAR tool (http://bayspar.geo.arizona.edu, last access: April 2018). For the BAYSPAR calculation we have assumed a prior mean 13 ^∘C and a search tolerance of 0.2 TEX₈₆ units. The calculated values are based on TEX₈₆ values from 62 latitude–longitude grid boxes that are 20^∘ by 20^∘.