Late Paleocene–early Eocene Arctic Ocean sea surface temperatures: reassessing biomarker paleothermometry at Lomonosov Ridge

A series of papers published shortly after the Integrated Ocean Drilling Program Arctic Coring Expedition (ACEX, 2004) on Lomonosov Ridge indicated remarkably high early Eocene sea surface temperatures (SSTs; ca. 23 to 27 ^∘C) and land air temperatures (ca. 17 to 25 ^∘C) based on the distribution of isoprenoid and branched glycerol dialkyl glycerol tetraether (isoGDGT and brGDGT) lipids, respectively. Here, we revisit these results using recent analytical developments – which have led to improved temperature calibrations and the discovery of new temperature-sensitive glycerol monoalkyl glycerol tetraethers (GMGTs) – and currently available proxy constraints.

The isoGDGT assemblages support temperature as the dominant variable controlling TEX₈₆ values for most samples. However, contributions of isoGDGTs from land, which we characterize in detail, complicate TEX₈₆ paleothermometry in the late Paleocene and part of the interval between the Paleocene–Eocene Thermal Maximum (PETM; ∼ 56 Ma) and the Eocene Thermal Maximum 2 (ETM2; ∼ 54 Ma). Background early Eocene SSTs generally exceeded 20 ^∘C, with peak warmth during the PETM (∼ 26 ^∘C) and ETM2 (∼ 27 ^∘C). We find abundant branched GMGTs, likely dominantly marine in origin, and their distribution responds to environmental change. Further modern work is required to test to what extent temperature and other environmental factors determine their distribution.

Published Arctic vegetation reconstructions indicate coldest-month mean continental air temperatures of 6–13 ^∘C, which reinforces the question of whether TEX₈₆-derived SSTs in the Paleogene Arctic are skewed towards the summer season. The exact meaning of TEX₈₆ in the Paleogene Arctic thus remains a fundamental issue, and it is one that limits our assessment of the performance of fully coupled climate models under greenhouse conditions.

1 Introduction

The Eocene epoch (56 to 34 million years ago; Ma) has long been characterized by warm climates. The earliest signs of a balmy Eocene Arctic region – fossil leaves of numerous plant species – were documented 150 years ago (Heer, 1869). Subsequent findings identified palms, baobab and mangroves, indicating the growth of temperate rainforests and year-round frost-free conditions in the Eocene Arctic region (Schweitzer, 1980; Greenwood and Wing, 1995; Suan et al., 2017; Willard et al., 2019). Fossils of animals, including varanid lizards, tortoises and alligators, also indicate warm Arctic climates (Dawson et al., 1976; Estes and Hutchinson, 1980). These earliest findings sparked interest in the climatological mechanisms allowing for such polar warmth about a century ago (Berry, 1922). Ever since, paleobotanists have focused on Arctic plant fossils and significantly refined their paleoclimatological interpretation towards estimates of precipitation as well as seasonal and mean annual temperature (e.g., Uhl et al., 2007; Greenwood et al., 2010; Eberle and Greenwood, 2012; Suan et al., 2017; Willard et al., 2019).

Novel insights in Paleogene Arctic paleoclimate research were made in the years following the Arctic Coring Expedition 302 (ACEX; Integrated Ocean Drilling Program – IODP, 2004; Fig. 1). This expedition recovered upper Paleocene and lower Eocene siliciclastic sediments deposited in a shallow marine environment, in Hole 4A (87^∘52.00^′ N; 136^∘10.64^′ E; 1288 m water depth), on the Lomonosov Ridge in the central Arctic Ocean (Backman et al., 2006). The succession was deposited at a paleolatitude of ∼ 78^∘ N based on a geological reconstruction (Seton et al., 2012) projected using a paleomagnetic reference frame (Torsvik et al., 2012) (see http://paleolatitude.org/, model version 2.1; Van Hinsbergen et al., 2015). The sediments are devoid of biogenic calcium carbonate but rich in immature organic matter, including terrestrial and marine microfossil assemblages and molecular fossils (e.g., Pagani et al., 2006; Sluijs et al., 2006; Stein et al., 2006).

Figure 1 Location of ACEX Hole 4A within a paleogeographic reconstruction of the Arctic region at the time of the PETM. Reconstruction made using gplates (Müller et al., 2018), with the tectonic reconstruction of Seton et al. (2012), the paleomagnetic reference frame of Torsvik et al. (2012) and modern coastlines. The red shape is Lomonosov Ridge in this reconstruction, and grey lines are structural features including spreading ridges.

As the upper Paleocene and lower Eocene sediments of the ACEX core lack biogenic calcium carbonate and alkenones, sea surface temperature (SST) reconstructions are based on the biomarker-based paleothermometer TEX₈₆. This proxy is based on membrane lipids (isoprenoid glycerol dibiphytanyl glycerol tetraethers; isoGDGTs) of Thaumarchaeota, which adapt the fluidity of their membrane according to the surrounding temperature by increasing the number of cyclopentane rings at higher temperatures (De Rosa et al., 1980; Wuchter et al., 2004; Schouten et al., 2013, and references therein). The proxy was introduced in 2002 by Schouten et al. (2002) and was calibrated to mean annual SST using modern marine surface sediments.

Initial papers suggested that Arctic SST increased significantly during two episodes of transient global warming. Maximum values of ∼ 23 and ∼ 27 ^∘C occurred during the Paleocene–Eocene Thermal Maximum (PETM-56 Myr ago; Sluijs et al., 2006) and the Eocene Thermal Maximum 2 (ETM2-54 Myr ago; Sluijs et al., 2009), respectively. Lower SSTs, generally exceeding 20 ^∘C, characterized the remainder of the early Eocene (Sluijs et al., 2008b). Such temperatures were immediately recognized to be remarkably high and could not be explained using fully coupled climate model simulations (Sluijs et al., 2006). Even the current-generation of IPCC-class (IPCC: Intergovernmental Panel on Climate Change) models are unable to match early Eocene Arctic mean annual SSTs, although reconstructions of tropical and midlatitude SSTs and deep ocean temperatures are consistent with some newer simulations (Frieling et al., 2017; Cramwinckel et al., 2018; Evans et al., 2018; Zhu et al., 2019).

Since the publication of the ACEX SST records, constraints on the applicability of the TEX₈₆ proxy have tremendously improved (see review by Schouten et al., 2013, and subsequent work by Taylor et al., 2013; Elling et al., 2014; Qin et al., 2014; Elling et al., 2015; Kim et al., 2015; Qin et al., 2015; Hurley et al., 2016; Zhang et al., 2016). This work has delivered new constraints on the ecology of Thaumarchaeota, the dominant depth at which they reside in the ocean and the depth from which their isoGDGTs are exported towards the seafloor. It also identified potential confounding factors such as variation in dominant isoGDGT export depth (e.g., Taylor et al., 2013; Kim et al., 2015), the input of non-thaumarchaeotal-derived isoGDGTs (e.g., Weijers et al., 2011; Zhang et al., 2011), growth phase (Elling et al., 2014), and environmental ammonium and oxygen concentrations (Qin et al., 2015; Hurley et al., 2016). Moreover, several indicators to detect such anomalies have been developed. Improvements in the chromatography method used for GDGT analysis now allow for better separation of previously co-eluting compounds, leading to enhanced analytical precision and sensitivity (Hopmans et al., 2016). Finally, recent work has described new GDGTs from oceans and sediments, notably branched glycerol monoalkyl glycerol tetraethers (brGMGTs or “H-shaped” brGDGTs) (e.g., Schouten et al., 2008; Liu et al., 2012), characterized by a covalent carbon–carbon bond that links the two alkyl chains. Their presence and distribution in peats and lake sediments have been linked to land air temperatures (LAT) (e.g., Naafs et al., 2018a; Baxter et al., 2019). However, these compounds have not yet been reported from ancient marine sediments.

Considering these developments and the paleoclimatological importance of the ACEX dataset, we reanalyzed the original lipid extracts for the PETM, ETM2 and the interval spanning these events (Sluijs et al., 2006; Sluijs et al., 2009) according to the latest chromatography protocols. We also compile published and generate new GDGT data from modern and Paleogene terrestrial deposits and use these to better assess the potential confounding influence of isoGDGTs from terrestrial sources, which was already recognized as a potential problem in the early work (Sluijs et al., 2006).

2 GDGT-based SST indices, calibration and confounding factors

2.1 TEX₈₆ and its calibration to SST

TEX₈₆ is based on the relative abundance of four different GDGTs (Fig. 2), following (Schouten et al., 2002)

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

where a higher relative abundance of cyclopentane moieties implies higher SSTs.

Figure 2 Molecular structures of the relevant isoGDGTs, brGDGTs and brGMGTs as well as their terminology as described in this study. Crenarchaeol isomer (not shown) differs from crenarchaeol in the stereochemistry of the cyclopentane moiety adjacent to the cyclohexyl moiety (Sinninghe Damsté et al., 2018b). For the terminology of the brGMGTs, for which the exact chemical structure is still unclear, we follow Baxter et al. (2019), since we identify the same isomers (see Fig. S2 for a chromatogram).

A number of models are used to calibrate TEX₈₆ to SST (Schouten et al., 2002; Schouten et al., 2003; Schouten et al., 2007; Kim et al., 2008; Liu et al., 2009; Kim et al., 2010; Tierney and Tingley, 2014; O'Brien et al., 2017), all based on a modern ocean surface sediment database. The currently available culture and mesocosm experiments as well as surface sediment data suggest that the relation between SST and TEX₈₆ is close to linear for a large portion of the modern ocean (Kim et al., 2010; Ho et al., 2014; Tierney and Tingley, 2014; O'Brien et al., 2017). In polar regions, the TEX₈₆ response to temperature diminishes (e.g., Kim et al., 2010; Tierney and Tingley, 2014). The response of TEX₈₆ to SST at the high-temperature end remains the subject of discussion (e.g., Cramwinckel et al., 2018; Hollis et al., 2019). Several authors prefer a linear relation (e.g., Tierney and Tingley, 2014; O'Brien et al., 2017). However, physiological considerations and multiple temperature-dependent GDGT indices might also imply a nonlinear relation at the high-temperature end, as can be observed at the high end of the modern ocean dataset and beyond the reach of the modern ocean in paleoclimate data (Cramwinckel et al., 2018). At higher temperatures, membrane adaptation may increasingly be established using isoGDGTs not included in the TEX₈₆ ratio, leading to a diminished TEX₈₆ response at very high temperatures (Cramwinckel et al., 2018). A nonlinear response has thus been proposed in other calibrations (Liu et al., 2009; Kim et al., 2010). The most recent nonlinear calibration, ${\mathrm{TEX}}_{\mathrm{86}}^{\mathrm{H}}$ (Kim et al., 2010), represents an exponential relation between SST and TEX₈₆ (Hollis et al., 2019). Unfortunately, ${\mathrm{TEX}}_{\mathrm{86}}^{\mathrm{H}}$ is mathematically problematic and has systematic residuals in the modern ocean (Tierney and Tingley, 2014).

Tierney and Tingley (2014) introduced a spatially varying Bayesian method to convert TEX₈₆ to SST and assume a linear relationship (BAYSPAR). BAYSPAR extracts TEX₈₆ values from the modern core-top dataset that are similar to the measured TEX₈₆ value from the geological sample based on a tolerance defined by the user, and it subsequently calculates regressions based on these core-top data. The uncertainty in SST reflects spatial differences in the correlation coefficient, intercept and error variance of the regression model.

Currently, it is generally encouraged to present results using both a linear and a nonlinear function (Hollis et al., 2019). The assumption of a linear or nonlinear relation between SST and TEX₈₆ leads to very different SST reconstructions for geological samples when TEX₈₆ values are > 0.70 (Kim et al., 2010; Tierney and Tingley, 2014; Frieling et al., 2017; O'Brien et al., 2017; Cramwinckel et al., 2018; Hollis et al., 2019). However, TEX₈₆ values for the early Eocene ACEX samples (0.5–0.7, Sluijs et al., 2006; Sluijs et al., 2008b; Sluijs et al., 2009) are below this value and well above most values observed in the polar regions (Kim et al., 2010; Tierney and Tingley, 2014; O'Brien et al., 2017), indicating that all calibrations will yield similar absolute SST values.

2.2 Caveats and confounding factors

Several confounding factors and caveats have been identified that could potentially bias TEX₈₆ data relative to mean annual SST. These notably relate to additions of isoGDGTs that were not produced in the upper water column by Thaumarchaeota, seasonal biases and choices that are made in the calibration between SST and TEX₈₆. Below we summarize methods that have been developed to assess if isoGDGT distributions might have been biased by confounding factors.

2.2.1 isoGDGTs of terrestrial origin

Previous work (Sluijs et al., 2006; Sluijs et al., 2008b; Sluijs et al., 2009) recognized that high contributions of terrestrially derived isoGDGTs could compromise the TEX₈₆ signal for portions of the upper Paleocene to lower Eocene interval of the ACEX core. This contribution can be tracked using the branched and isoprenoid tetraether (BIT) index, a ratio of mostly soil-derived branched GDGTs (brGDGTs; Fig. 2) and crenarchaeol, which is dominantly marine-derived (Hopmans et al., 2004; Schouten et al., 2013; adjusted to include 6-methyl brGDGT isomers):

$$\begin{array}{}\text{(2)}& \mathrm{BIT}\mathrm{index}={\displaystyle \frac{\left({\displaystyle \begin{array}{c}\left[\mathrm{brGDGT}-\mathrm{Ia}\right]+\left[\mathrm{brGDGT}-\mathrm{IIa}\right]+\left[\mathrm{brGDGT}-{\mathrm{IIa}}^{\prime }\right]\\ +\left[\mathrm{brGDGT}-\mathrm{IIIa}\right]+\left[\mathrm{brGDGT}-{\mathrm{IIIa}}^{\prime }\right]\end{array}}\right)}{\left(\begin{array}{c}\left[\mathrm{brGDGT}-\mathrm{Ia}\right]+\left[\mathrm{brGDGT}-\mathrm{IIa}\right]+\left[\mathrm{brGDGT}-{\mathrm{IIa}}^{\prime }\right]\\ +\left[\mathrm{brGDGT}-\mathrm{IIIa}\right]+\left[\mathrm{brGDGT}-{\mathrm{IIIa}}^{\prime }\right]+\left[\mathrm{crenarchaeol}\right]\end{array}\right).}}\end{array}$$

Most studies define a BIT value (typically 0.3 or 0.4) above which TEX₈₆-derived SSTs are unreliable (e.g., Weijers et al., 2006). However, the threshold of 0.4 is conservative in some settings, and the impact of terrigenous GDGTs on reconstructed SST will depend on the nature and temperature of the source catchment (Inglis et al., 2015). In addition, a cutoff value based on BIT values is difficult given the relatively large differences in BIT between labs, which originate from methodological differences (Schouten et al., 2009). A strong linear relationship between BIT and TEX₈₆ values is often taken as an indication of a bias in TEX₈₆ through land-derived isoGDGTs to the marine TEX₈₆ signature (e.g., Douglas et al., 2014). An earlier study used a somewhat subjective threshold of 0.3 for an interval spanning ETM2 in the ACEX core (Sluijs et al., 2009).

2.2.2 isoGDGTs of deepwater origin

Thaumarchaeota, the source of most isoGDGTs in marine waters (Zeng et al., 2019; Besseling et al., 2020), are ammonium oxidizers (Könneke et al., 2005; Wuchter et al., 2006a), making them independent of light. Although they occur throughout the water column, maximum abundances occur at depths < 200 m, generally around NO₂ maxima (e.g., Karner et al., 2001; Pitcher et al., 2011a). In most oceans, sedimentary GDGTs dominantly derive from the upper few hundred meters based on analyses of suspended particular organic matter and sediment traps (Wuchter et al., 2005; Wuchter et al., 2006b; Yamamoto et al., 2012; Richey and Tierney, 2016). A deeper contribution has also been inferred based on ¹⁴C analysis (Shah et al., 2008), implying possible contributions of isoGDGTs from the thermocline. Moreover, contributions of isoGDGTs produced in the deep sea have regionally been identified (e.g., Kim et al., 2015). Taylor et al. (2013) also found that deep-dwelling (> 1000 m) archaea might contribute to the sedimentary isoGDGT assemblage. They indicate that such deep contributions can be tracked using the GDGT-2 ∕ GDGT-3 ratio; high values of > 5 indicate contributions of archaea living deeper in the water column. Given that upper Paleocene and lower Eocene ACEX sediments were deposited in a shallow shelf environment (Sluijs et al., 2008b), a significant contribution of deep ocean archaeal lipids is not expected.

2.2.3 isoGDGTs of methanotrophic and methanogenic archaea

Contributions of isoGDGTs to the sedimentary pool might also derive from anaerobic methanotrophs and/or methanogens. Several indices have been developed to track such contributions based on relatively high contributions of particular isoGDGTs of these groups of archaea. The methane index (MI) was developed to detect the relative contribution of anaerobic methanotrophic Euryarchaeota, notably represented by GDGT-1, GDGT-2 and GDGT-3 (Pancost et al., 2001; Zhang et al., 2011), and is therefore defined as

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

MI values greater than 0.5 indicate a significant contribution of anaerobic methanotrophy. Such values may yield unreliable TEX₈₆ values. Another tracer for contributions of anaerobic methanotrophic archaea is the analogous GDGT-2 ∕ crenarchaeol ratio (Weijers et al., 2011).

Methanogenic archaea can synthesize GDGT-0 and smaller quantities of GDGT-1, GDGT-2 and GDGT-3. The ratio GDGT-0 ∕ crenarchaeol is indicative of contributions of methanogenic archaea to the isoGDGT pool (Blaga et al., 2009); values > 2 indicate a substantial contribution of methanogenic archaea. Up to now, high index values have often been observed near methane seeps or anoxic basins (e.g., Jaeschke et al., 2012) but rarely in open marine waters in the modern domain and paleodomains (Inglis et al., 2015; Zhang et al., 2016). Given the reducing conditions in the sediment and water column at the study site across the late Paleocene and early Eocene (Sluijs et al., 2006; Stein et al., 2006; Sluijs et al., 2008b; März et al., 2010), an influence of methane cycling might be expected.

2.2.4 isoGDGTs of the “Red Sea type”

Sedimentary isoGDGT distributions from the Red Sea are anomalous to other marine settings and are characterized by the low abundance of GDGT-0 and the high abundance of the crenarchaeol isomer. Presumably, this is due to an endemic thaumarchaeotal assemblage. The Red Sea isoGDGT distribution yields a different relationship between SST and TEX₈₆ (Trommer et al., 2009; Kim et al., 2015). Inglis et al. (2015) attempted to quantify a “Red Sea type” GDGT distribution in geological samples using the following index:

$$\begin{array}{}\text{(4)}& \mathit{\%}\mathrm{GDGTrs}={\displaystyle \frac{\left[\mathrm{crenarchaeol}\mathrm{isomer}\right]}{(\left[\mathrm{GDGT}-\mathrm{0}\right]+\left[\mathrm{crenarchaeol}\mathrm{isomer}\right])}}\times \mathrm{100}.\end{array}$$

However, as noted by Inglis et al. (2015) this ratio is also strongly SST-dependent such that the Red Sea type GDGT assemblage cannot be discerned from GDGT distributions that occur at high temperatures in normal open marine settings.

2.2.5 Seasonal bias

TEX₈₆ is calibrated to mean annual SST. However, particularly in midlatitude and high-latitude areas where production and export production are highly seasonal, the sedimentary GDGT distribution might not represent annual mean conditions (Wuchter et al., 2006b; Pitcher et al., 2011b; Mollenhauer et al., 2015; Richey and Tierney, 2016; Park et al., 2019). This issue should partly be reflected in the calibration uncertainty of the modern ocean database (several degrees Celsius, depending on the calibration and method; see section 2.7). Sluijs et al. (2006, 2008b, 2009) originally argued that the TEX₈₆ results from the ACEX core could be biased towards summer temperature because the export of organic matter from the surface ocean towards the sediment likely peaked during the season of highest production, i.e., the summer. However, we also note that the TEX₈₆–temperature relationship is not improved when using seasonal mean ocean temperatures (Kim et al., 2010; Tierney and Tingley, 2014) and modern observations indicate homogenization of the seasonal cycle at depth (Wuchter et al., 2006b; Yamamoto et al., 2012; Richey and Tierney, 2016), implying that seasonality has a relatively limited effect on modern sedimentary TEX₈₆ values.

2.2.6 Additional isoGDGT-based temperature indicators

The underlying mechanism of TEX₈₆ is that isoGDGTs produced at higher SSTs contain more rings than those produced at low SSTs. Although the combination of compounds included in TEX₈₆ seems to yield the strongest relation with temperature in the modern ocean (Kim et al., 2010), it implies that isoGDGT ratios other than TEX₈₆ also provide insights into SST. One alternative temperature-sensitive isoGDGT index is the ring index (RI), which represents the weighed number of cyclopentane rings of isoGDGTs 0–3, crenarchaeol and the crenarchaeol isomer (Zhang et al., 2016), defined as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}\mathrm{RI}& =\mathrm{0}\times \left[\mathit{\%}\mathrm{GDGT}-\mathrm{0}\right]+\mathrm{1}\times \left[\mathit{\%}\mathrm{GDGT}-\mathrm{1}\right]\\ & +\mathrm{2}\times \left[\mathit{\%}\mathrm{GDGT}-\mathrm{2}\right]+\mathrm{3}\times \left[\mathit{\%}\mathrm{GDGT}-\mathrm{3}\right]\\ & +\mathrm{4}\times \left[\mathit{\%}\mathrm{crenarchaeol}+\mathit{\%}\mathrm{crenarchaeol}\mathrm{isomer}\right].\end{array}\end{array}$$

Note that the abundance of GDGT-0 is important for determining the percentage of the other GDGTs in the total isoGDGT pool.

The close relation between TEX₈₆ and RI can also be used to detect aberrant distributions, including those produced by methanogenic, methanotrophic and terrestrial sources, as these sources typically contribute disproportionate amounts of specific lipids. RI_TEX, calculated from TEX using the polynomial fit of Zhang et al. (2016), is subtracted from the RI to arrive at the ΔRI. Cutoff values for sample deviation from the modern ocean calibration dataset are defined as 95 % confidence limits of the TEX₈₆–RI relation, or above $\left|\mathrm{0.3}\right|\mathrm{\Delta }$RI units.

2.3 H-shaped branched GDGTs; brGMGTs

BrGMGTs (Fig. 2) were first identified by Liu et al. (2012) in marine sediments, who identified a single acyclic tetramethylated brGMGT (m∕z 1020). This compound was later detected within the marine water column and appeared to be abundant within the oxygen minimum zone (Xie et al., 2014). Naafs et al. (2018a) identified a larger suite of brGMGTs (including m∕z 1048 and 1034) in a quasi-global compilation of modern peat samples. They argued that these compounds were preferentially produced at depth within the anoxic catotelm. Analogous to the continental paleothermometer based on bacterial brGDGTs produced in surface soils, termed MBT'_5me (Weijers et al., 2007b; De Jonge et al., 2014), they showed that the degree of methylation of brGMGTs in peats relates to mean annual air temperature. They calculated the degree of methylation of brGDGTs without cyclopentane moieties, designed for comparison to the methylation of brGMGTs, defined by H-MBT_acyclic:

$$\begin{array}{}\text{(6)}& {\mathrm{MBT}}_{\mathrm{acyclic}}={\displaystyle \frac{\mathrm{brGDGT}-\mathrm{Ia}}{\left(\begin{array}{c}\mathrm{brGDGT}-\mathrm{Ia}+\mathrm{brGDGT}-\mathrm{IIa}+\mathrm{brGDGT}-{\mathrm{IIa}}^{\prime }\\ +\mathrm{brGDGT}-\mathrm{IIIa}+\mathrm{brGDGT}-{\mathrm{IIIa}}^{\prime }\end{array}\right)}},\end{array}$$

$$\begin{array}{}\text{(7)}& \mathrm{H}-{\mathrm{MBT}}_{\mathrm{acyclic}}={\displaystyle \frac{\mathrm{H}\mathrm{1020}}{\left(\begin{array}{c}\mathrm{H}\mathrm{1020}\\ \mathrm{H}\mathrm{1034}+\mathrm{H}\mathrm{1048}\end{array}\right)}}.\end{array}$$

Based on the strong relation between MBT_acyclic and H-MBT_acyclic in their peat samples, Naafs et al. (2018a) suggested that the brGMGTs have the same origin as the brGDGTs, presumably Acidobacteria (Sinninghe Damsté et al., 2011; Sinninghe Damsté et al., 2018a). In addition, they showed that the abundance of brGMGTs (relative to the total quantity of brGMGTs and brGDGTs) positively correlates with mean annual air temperature, suggesting that the covalent bond in the brGMGTs is used to maintain membrane stability at higher temperature (Naafs et al., 2018a).

Baxter et al. (2019) identified a total of seven different brGMGTs from a suite of African lake sediments (Fig. 2) and found their relative distribution to correlate with mean annual air temperature. Accordingly, they proposed a proxy for mean annual air temperature termed brGMGT-I (see Fig. 2 for the molecular structures referred to here):

$$\begin{array}{}\text{(8)}& \mathrm{brGMGT}-\mathrm{I}={\displaystyle \frac{\left[\mathrm{H}\mathrm{1020}\mathrm{c}\right]+\left[\mathrm{H}\mathrm{1034}\mathrm{a}\right]+\left[\mathrm{H}\mathrm{1034}\mathrm{c}\right]}{\begin{array}{c}\left[\mathrm{H}\mathrm{1020}\mathrm{b}\right]+\left[\mathrm{H}\mathrm{1020}\mathrm{c}\right]+\left[\mathrm{H}\mathrm{1034}\mathrm{a}\right]\\ +\left[\mathrm{H}\mathrm{1034}\mathrm{c}\right]+\left[\mathrm{H}\mathrm{1048}\right]\end{array}}}.\end{array}$$

3 Material and methods

We used the polar fractions previously analyzed by Sluijs et al. (2006, 2009) from the PETM through ETM2 interval at IODP Expedition 302 Hole 4A. These fractions originate from a total lipid extract produced using a Dionex accelerated solvent extractor and fraction separations by Al₂O₃ column chromatography using hexane:dichloromethane (DCM) (9:1, v∕v) and DCM:methanol (1:1; v∕v) to yield the apolar and polar fractions, respectively. Polar fractions were redissolved in hexane:isopropanol (99:1, v∕v) and passed through a 0.45 µm polytetrafluoroethylene filter. This fraction was then analyzed by high-performance liquid chromatography (HPLC) and atmospheric pressure chemical ionization–mass spectrometry using an Agilent 1260 Infinity series HPLC system coupled to an Agilent 6130 single-quadrupole mass spectrometer at Utrecht University following Hopmans et al. (2016) to measure the abundance of GDGTs. Based on long-term observation of the in-house standard, the analytical precision for TEX₈₆ calculates to ±0.3 ^∘C in the SST domain.

To gain further insights into the potential impact of terrestrial isoGDGT input on TEX₈₆ values, we compiled isoGDGT and brGDGT distributions from modern peats (n=473; Naafs et al., 2017) and early Paleogene lignites (n=58; Naafs et al., 2018b). Note that the fractional abundance of crenarchaeol isomer was not reported in the early Paleogene dataset of Naafs et al. (2018b). We therefore revisited the original chromatograms from Naafs et al. (2018b) and integrated the crenarchaeol isomer (m∕z 1292).

4 Results

The new GDGT distributions (Table S1 in the Supplement) are consistent with the TEX₈₆ and BIT index data generated over a decade ago using the older analytical HPLC setup (Hopmans et al., 2000; Hopmans et al., 2016) (Fig. 3). TEX₈₆ exhibits some scatter and the slope of the regression is slightly off the 1:1 line, indicating that the new data have somewhat higher TEX₈₆ values. Less scatter is apparent in the BIT record, but the original BIT index values were slightly higher than recorded here, as indicated by the regression (Fig. 3). This result is consistent with previous analyses with the new analytical setup (Hopmans et al., 2016). This does not impact previous qualitative interpretations of this record (Sluijs et al., 2006; Sluijs et al., 2008b; Sluijs et al., 2009). In the Discussion section, we assess indicators of potential confounding factors (Sect. 2.2), including the influx of terrestrially derived isoGDGTs to the sediments (Figs. 4, 5 and S1) and several indices related to methane and depth of production (Fig. 6).

Figure 3 Comparison of the original GDGT dataset for the upper Paleocene and lower Eocene of ACEX Hole 4A (Sluijs et al., 2006; Sluijs et al., 2009) and the new data generated according to the latest chromatography protocols.

Although we did not detect significant quantities of isoprenoid GMGTs, high abundances of various brGMGTs are present in the ACEX samples, in total between 10 % and 45 % of the total brGDGT assemblage (Fig. 7). We consistently identify at least five brGMGTs across the three different mass-to-charge ratios (m∕z 1020, 1034 and 1048). Based on their (relative) retention times and overall distribution we were able to apply the nomenclature of Baxter et al. (2019) to five of these and assign individual peaks to previously identified compounds (Fig. S2). The abundance of brGMGTs relative to brGDGTs increases during the PETM. The proposed temperature indicators based on brGMGTs show mixed results, with some showing a clear response to the PETM (Fig. 7e) while others do not (Fig. 7d).

5 Discussion

5.1 IsoGDGT provenance

5.1.1 Contributions of soil-derived isoGDGTs

As noted by Sluijs et al. (2006), late Paleocene samples yield anomalously high abundances of GDGT-3, likely derived from a terrestrial source. We therefore consider the late Paleocene temperature estimates unreliable. To assess the temperature change during the PETM, Sluijs et al. (2006) developed a TEX₈₆ calibration without this moiety, termed TEX${}_{\mathrm{86}}^{\prime }$. However, TEX${}_{\mathrm{86}}^{\prime }$ has not been widely used outside the Paleogene Arctic because the anomalous abundances of GDGT-3 have not been recorded elsewhere. High contributions of GDGT-3 from terrestrial input would also be associated with an increase in the abundance of other isoGDGTs. Indeed, recent TEX₈₆-based global SST compilations and comparisons to climate simulations for the PETM excluded the Paleocene ACEX data because the TEX${}_{\mathrm{86}}^{\prime }$ calibration complicates the comparison to other regions where it has not been applied (Frieling et al., 2017; Hollis et al., 2019).

Input of soil organic matter is consistent with Willard et al. (2019), who established that the brGDGT assemblage is dominantly soil-derived as opposed to being produced in the coastal marine environment. This observation is based upon the weighted average number of rings in the tetramethylated brGDGTs (#rings_tetra), which generally does not exceed 0.4 to 0.7 in the global soil calibration dataset (Sinninghe Damsté, 2016). In the ACEX record, #rings_tetra is < 0.21 (Willard et al., 2019), consistent with a dominant soil source. This indicates that (1) brGDGT abundances, (2) brGDGT distributions and (3) the BIT index are reliable indicators of the relative supply of terrestrially derived isoGDGTs into the marine basin. The Paleocene section of the dataset also stands out regarding its relation between the BIT index and TEX₈₆ (Fig. 4), which confirms its anomalous nature.

Figure 4 Comparison between BIT index values and TEX₈₆ for various intervals spanning the upper Paleocene and lower Eocene of ACEX Hole 4A.

During the PETM, TEX₈₆ values are higher (due to warming) and BIT values lower. This was attributed to sea level rise during the hyperthermals, resulting in a more distal position relative to the terrestrial GDGT source (Sluijs et al., 2006; Sluijs et al., 2008a). The interval between 371.0 and 369.0 mcd (i.e., above the PETM and below ETM2) also stands out. This interval was previously recognized by Sluijs et al. (2009) to reflect an open marine environment, with a dominance of marine palynomorphs and algal biomarkers. They also found that high BIT values correspond to low TEX₈₆ values within that interval and therefore implemented a subjective threshold value of 0.3, above which TEX₈₆-derived SSTs were considered unreliable. Although the relation between BIT and TEX₈₆ exhibits considerable scatter, the new analyses support the notion that a higher influx of terrestrial isoGDGTs lowers TEX₈₆ values. The linear regression (Fig. 4; excluding the one outlier with high TEX₈₆ and BIT values in the top right of the plot because it has highly anomalous distributions; ΔRI = 0.61) yields an R² of 0.26 that explains a portion of the variation (Fig. 4). The nature of this influence is determined by the relative abundance of terrestrial isoGDGTs and their TEX₈₆ value. The TEX₈₆ value at the terrestrial endmember of BIT = 1, assuming various types of regressions, centers around 0.5. The remainder of the data does not show a clear relation between BIT and TEX₈₆, although some of the lowest TEX₈₆ values correspond to high BIT values, suggesting that the terrestrial endmember contributed isoGDGT assemblages with relatively low TEX₈₆ values in other intervals as well.

Figure 5 The abundance of various isoGDGTs relative to the total brGDGT abundance in modern peats (a, b) and Paleogene lignites (a; Eq. 9) used to assess potential isoGDGT contributions to the ACEX samples. The box represents standard 25 %–50 %–75 % quantiles, and whiskers represent the box limits plus or minus 1.5× the interquartile range (IQR). Any data outside that range are given as circles. The number of measurements per dataset are as follows. Modern peats: 473 (most isoGDGTs have been identified in ±430 of those); modern peats above 15 ^∘C: 141 (all except one have isoGDGT data); lignites: 58 (all of which have isoGDGT data but only 29 have available and quantifiable crenarchaeol isomer data).

The relatively low degree of cyclization in the early Eocene contrasts starkly with the high degree of cyclization during the late Paleocene (Fig. 6e). This implies that the distribution of terrestrial isoGDGTs varies strongly between the latest Paleocene and early Eocene within our studied section.

The impact of soil-derived isoGDGTs also emerges from the ring index approach of Zhang et al. (2016; see Sect. 2.2.6 and Fig. 6). The difference between the ring index and TEX₈₆ at the onset of the PETM is mainly controlled by crenarchaeol, which is comparatively low in abundance in the Paleocene but highly abundant in the PETM. This increase is likely associated with sea level rise during the PETM because crenarchaeol is predominantly produced in the marine realm. It is also consistent with a drop in BIT index values and the relative abundance of terrestrial palynomorphs (Sluijs et al., 2008a). The approach of Zhang et al. (2016) also confirms that many isoGDGT distributions exhibit an anomalous relation between TEX₈₆ and the ring index relative to the modern core top dataset, with ΔRI values > 0.3 (Fig. 6). Importantly, all samples with ΔRI values > 0.3 have BIT values above 0.35, indicating that contributions of soil-derived iso-GDGTs dominate non-temperature effects in the distributions. We therefore discard TEX₈₆-derived SSTs for samples with BIT values > 0.35.

Figure 6 Branched and isoprenoid GDGT records across the upper Paleocene and lower Eocene of ACEX Hole 4A. (a) Carbon isotope stratigraphy (total organic carbon record from Sluijs et al., 2006, 2009; marine organic matter record from Sluijs and Dickens, 2012), (b) BIT index (Eq. 2), (c) indices indicative of anaerobic archaeal methanotrophy (MI (Eq. 3) and GDGT-2 ∕ crenarchaeol) and methanogenesis (GDGT-0 ∕ crenarchaeol), (d) GDGT-2 ∕ GDGT-3 ratio, (e) ring index (Eq. 5) and D-ring index, and (f) TEX86 (Eq. 1) calibrated to sea surface temperature using a nonlinear calibration TEX${}_{\mathrm{86}}^{\mathrm{H}}$ calibration (Kim et al., 2010) and the BAYSPAR method, which is based on a linear calibration (Tierney and Tingley, 2014).

We also develop a crude model to further constrain the potential contribution of terrestrially derived isoGDGTs. First, we determine the abundance of isoGDGTs relative to brGDGTs in modern peat samples (Naafs et al., 2017) and early Paleogene lignites (fossil peat) (Naafs et al., 2018b; the isoGDGT data are published here). Although there is no reason to assume that peat was a major component of the hinterland (Willard et al., 2019), the aforementioned datasets can provide an estimate of the potential contribution from terrestrial isoGDGTs to our study site. The raw signal intensity of brGDGTs in the ACEX samples is used to estimate the potential contribution of terrestrially derived isoGDGTs to the samples. To this end, we use the fractional abundance of the various isoGDGTs in the modern peat and Paleogene lignite datasets (Fig. 5). Then, we estimate the abundance of these terrestrially derived isoGDGTs in our ACEX samples by scaling this fraction to the measured abundances of brGDGTs and isoGDGTs in our ACEX samples, following

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathrm{Terrestrial}& \mathrm{fraction}\mathrm{of}\mathrm{isoGDGT}x\\ & =\left(\mathrm{fraction}\mathrm{of}\mathrm{isoGDGT}x\mathrm{in}\mathrm{terrestrial}\mathrm{test}\mathrm{dataset}\right.\\ & \left.\times {\displaystyle \frac{\mathrm{sum}\left(\mathrm{brGDGTs}\right)}{\mathrm{abundance}\mathrm{of}\mathrm{isoGDGT}x}}\right),\end{array}\end{array}$$

where x represents the specific analyzed GDGT (see the Supplement for an example of these calculations).

This leads to estimates of the potential relative contributions of the individual isoGDGTs derived from land to the ACEX samples based on the entire modern peat dataset (Naafs et al., 2017), modern peats from regions with mean annual air temperature (MAT) exceeding 15 ^∘C (Naafs et al., 2017) and Paleogene lignites (Naafs et al., 2018b; this paper, Figs. 5 and S1). This approach implies that crenarchaeol and the crenarchaeol isomer are almost exclusively from the marine realm. However, GDGT-1, GDGT-2 and GDGT-3 in our study site may be derived from the terrestrial realm (Fig. 5), especially in specific stratigraphic intervals (Fig. S1). In the most extreme cases, the modeled contributions of terrestrial isoGDGTs are higher than the measured isoGDGT abundances (i.e., terrestrial fraction > 1). This is principally seen in isoGDGT-2 and GDGT-3, especially when we employ the Paleogene lignite database. This particular assumption clearly overestimates the abundance of terrestrially sourced isoGDGTs in our setting. However, the temporal trends obtained using modern peats, subtropical modern peats and Paleogene lignites are essentially identical and give some indication of which isoGDGTs are most likely to be impacted by terrestrial input and across which intervals. Interestingly, this approach also suggests that particularly GDGT-3 is shown to be strongly affected (Fig. 5), which qualitatively matches the distributions in the ACEX samples. This is principally because GDGT-3 is the least abundant marine isoGDGT included in our analyses, whereas it is often as abundant as GDGT-1 and GDGT-2 in terrestrial settings (Fig. 5).

5.1.2 Contributions of methanotrophic or methanogenic archaea?

The depositional environment at the study site included ample (export) production, sediment organic matter content and low-oxygen conditions at the sediment–water interface (Sluijs et al., 2006; Stein et al., 2006; Stein, 2007; Sluijs et al., 2008b; Sluijs et al., 2009; März et al., 2010). This may have been suitable for abundant methanogenic and methanotrophic archaea, potentially contributing to the sedimentary isoGDGT assemblage. However, our GDGT-2 ∕ crenarchaeol values (< 0.23; Fig. 6) are far below values that suggest significant isoGDGT contributions of methanotrophic Euryarchaeota as described by Weijers et al. (2011). MI values (maximum observed 0.31) are also generally below proposed cutoff values (0.3–0.5, Zhang et al., 2011) that suggest such contributions. Finally, GDGT-0 ∕ crenarchaeol ratios (< 1.4) remain below the cutoff value of 2 throughout the section (Fig. 6), also making a significant isoGDGT contribution from methanogens highly unlikely (Blaga et al., 2009). Collectively, relative contributions of isoGDGTs from methanogenic and methanotrophic archaea seem low despite the low-oxygen environment, suggesting a relatively high flux of pelagic isoGDGTs.

5.1.3 Contributions of deep-dwelling archaea?

Taylor et al. (2013) showed that GDGT-2 ∕ GDGT-3 ratios correspond to depth of production, with high values (> 5) in deep waters (> 1000 m). We record low values (1–4) between ∼ 390 and ∼ 371.2 mcd (Fig. 6), which supports dominant production in the surface ocean based on the modern calibration dataset (Taylor et al., 2013). However, the overlying interval (∼ 371 to ∼ 368.3 mcd) has much higher (average 7.4) and variable GDGT-2 ∕ GDGT-3 values, with peak values of 10–14. Such values suggest significant contributions of isoGDGTs produced at water depths of several kilometers according to the analyses by Taylor et al. (2013).

However, all paleoenvironmental information generated based on the sediments and tectonic reconstructions of Lomonosov Ridge – a strip of continental crust that disconnected from the Siberian margin in the Paleocene – has indicated a neritic setting of the drill site at least up to the middle Eocene (e.g., O'Regan et al., 2008; Sangiorgi et al., 2008; Sluijs et al., 2008a; Sluijs et al., 2009). At ∼ 371.2 mcd a drop in the BIT index and a change in the palynological assemblages correspond to an interval of greenish sediment suggestive of pronounced amounts of glauconite. These changes are consistent with local relative sea level rise, causing a somewhat more distal position relative to the shoreline. However, the sediment remains dominantly siliciclastic and organic terrestrial components, particularly pollen and spores, still remain abundant, indicating a shallow setting (Sluijs et al., 2008a, 2008b). Increased contributions of isoGDGTs produced at depth would be expected to have caused a systematic cold bias, but based on linear regression analysis the large variability in GDGT-2 ∕ GDGT-3 ratios is unrelated to the recorded variability in TEX₈₆ values. The high GDGT-2 ∕ GDGT-3 ratio values therefore cannot be explained by contributions of deep-dwelling archaea.

In a study of the last 160 kyr in the South China Sea, Dong et al. (2019) found that very high GDGT-2 ∕ GDGT-3 ratios (∼ 9 but up to 13) correspond to high values in nitrogen isotope ratios, interpreted to reflect low contributions in diazotroph N₂ fixation and enhanced upwelling. In our record, the high GDGT-2 ∕ GDGT-3 ratios are associated with normal marine conditions and the dinocyst assemblages are not indicative of upwelling conditions (Sluijs et al., 2009). Unfortunately, the available nitrogen isotope record (Knies et al., 2008) does not cover this interval in sufficient resolution to assess a relation with diazotroph activity. The increase in the GDGT-2 ∕ GDGT-3 ratio correlates with a strong drop in BIT index values and an increase in normal marine dinocyst species (Sluijs et al., 2009), but a shift to more open marine environment does not explain the high ratio values. As such, the cause of the high GDGT-2 ∕ GDGT-3 ratios in this interval remains unclear, but we consider it highly unlikely to relate to contributions of deep-dwelling Thaumarchaeota.

5.1.4 Oxygen concentrations and ammonium oxidation rates

A variety of nonthermal factors can impact TEX₈₆ values, including ammonium and oxygen concentrations and growth phase (Elling et al., 2014; Qin et al., 2014; Hurley et al., 2016). Across the studied interval of the ACEX core, several intervals of seafloor and water column anoxia have been identified based on organic and inorganic proxies, notably during the PETM and ETM2 (Sluijs et al., 2006; Stein et al., 2006; Sluijs et al., 2008b; Sluijs et al., 2009; März et al., 2010).

Particularly suspect is an interval of low TEX₈₆ values that marks the middle of the ETM2 interval, directly following a ∼ 4 ^∘C warming at its onset (Sluijs et al., 2009). This interval is also marked by the presence of sulfur-bound isorenieratane (Sluijs et al., 2009), a derivative of isorenieratene. This biomarker is produced by the brown strain of green sulfur bacteria that require light for photosynthesis and free sulfide, indicating euxinic conditions in the (lower) photic zone (Sinninghe Damsté et al., 1993). We also record a concomitant shift in several methane-related indicators, GDGT-2 ∕ GDGT-3 ratio values and the ΔRI. A mid-ETM2 cooling signal has not been recorded at other study sites, and this interval marks the occurrence of pollen of thermophilic plants such as palms and baobab (Sluijs et al., 2009; Willard et al., 2019). Therefore, the low TEX₈₆ values were suggested to reflect thaumarchaeotal depth migration to the deeper chemocline due to euxinic conditions (Sluijs et al., 2009), similar to the modern Black Sea (Coolen et al., 2007; Wakeham et al., 2007) and the Mediterranean Sea during sapropel formation (Menzel et al., 2006).

More recent work has indicated that the isolated marine thaumarchaeotal species Nitrosopumilus maritimus produces lower TEX₈₆ values with higher ammonia oxidation rates (Hurley et al., 2016) and O₂ concentrations (Qin et al., 2015). Although this observation is difficult to extrapolate to the total response of the thaumarchaeotal community in the marine environment on geological timescales, lower O₂ availability should lower oxidation rates, leading to higher TEX₈₆ values (Qin et al., 2015; Hurley et al., 2016). However, we record a drop in TEX₈₆ values with the development of anoxia during ETM2. The nature of the anomalously low cyclization in the ETM2 isoGDGT assemblage, which passes all quality tests regarding GDGT distribution (Fig. 6), therefore remains elusive.

5.2 Origin and environmental forcing of brGMGTs

The relative abundances of brGMGTs in our samples are surprisingly high. On average, they comprise 25 % of the total branched GDGT and GMGT assemblage. The limited literature on modern occurrences implies that both terrestrial and marine sources may have contributed to the brGMGT assemblage. Data from marine sediments (Liu et al., 2012) and the water column (Xie et al., 2014) clearly show production within the marine realm. Their occurrence in modern peats (Naafs et al., 2018a), lake sediments (Baxter et al., 2019) and Paleogene lignites (Inglis et al., 2019) might also imply transport from land to marine sediments. A soil-derived source is currently unsupported, as they were most often below the detection limit in recent studies of geothermally heated soils (De Jonge et al., 2019) and a soil transect from the Peruvian Andes (Kirkels et al., 2020). The relative brGMGT abundances we record are close to the maximum found in modern peats (Naafs et al., 2018a). However, significant input of peat-derived organic matter into our study site is inconsistent with the low input of peat-derived Sphagnum spores (Willard et al., 2019). Alternatively, the high abundance of brGMGTs could also be related to subsurface production in marine sediments. An analogous process was invoked by Naafs et al. (2018a) to explain the very high abundance of brGMGTs in an early Paleogene lignite. Collectively, however, we surmise that production in the marine realm may be an important contributor to the brGMGT pool in our setting.

Several factors may contribute to the rise in the abundance of brGMGTs relative to brGDGTs across the PETM. Higher relative abundances of brGMGTs in modern peats generally occur at higher mean annual air temperatures (Naafs et al., 2018a), so this signal could relate to warming during the PETM if their origin at the study site is terrestrial. However, since we consider it likely that a large part of the brGMGT assemblage is of marine origin, the rise in brGMGT abundance likely relates to the previously recorded (Sluijs et al., 2006; Sluijs et al., 2008b) sea level rise during the PETM at the study site. This is consistent with the increase in marine brGMGT production relative to terrestrial brGDGT supply to the study site (Fig. 7b). This is consistent with the inverse correlation between brGMGT abundance and the BIT index (Fig. 7b). Lastly, if the production of marine brGMGTs was focused in oxygen minimum zones (Xie et al., 2014), the development of low-oxygen conditions in the water column based on several indicators, such as the presence of isorenieratane (Sluijs et al., 2006), might have increased the production of brGMGTs in the water column. It is also possible that all of these factors contributed to the changes in the abundance of brGMGTs relative to brGDGTs across the PETM.

Figure 7 Branched GMGT records across the upper Paleocene and lower Eocene of ACEX Hole 4A. (a) Carbon isotope stratigraphy (total organic carbon record from Sluijs et al., 2006, 2009; marine organic matter record from Sluijs and Dickens, 2012). (b) Fraction of brGMGTs of the total branched GDGTs and GMGTs and BIT index (Eq. 2). (c) MBT'_5me record (Willard et al., 2019) using the calibrations of De Jonge et al. (2014) as well as BayMBT (Dearing Crampton-Flood et al., 2020; median value; see the Supplement for uncertainty) and TEX${}_{\mathrm{86}}^{\mathrm{H}}$. (d) MBT_acyclic (Eq. 6) and H-MBT based on all isomers detected with m∕z 1020 and m∕z 1034 (H-MBT all; Eq. 7) and based on H1020a and H1034b (H-MBT H1020a, H1034c). (e) The brGMGT-1 record (Eq. 8).

The brGMGT-I proxy does not produce temperature trends similar to those seen in TEX₈₆ or MBT'_5me (Fig. 7d). If the majority of the brGMGTs are of marine origin, this indicates that brGMGTs produced in the marine realm do not respond to temperature as was hypothesized based on the African lake dataset by Baxter et al. (2019).

Also, the application of the H-MBT_acyclic index (equation 7) appeared problematic because, similar to Baxter et al. (2019), we identified several more isomers than Naafs et al. (2018a, who developed this index) detected in their peat samples. It therefore remains unclear which of our peaks should be used to calculate the H-MBT_acyclic index values. We therefore show the two plausible options. For the first, we use all peaks with m∕z 1020, 1034 and 1048 (H-MBT-all in Fig. 7e) within the expected retention time window. However, based on our chromatography, we consider it more likely that the dominant peaks identified by Naafs et al. (2018a) at m∕z 1020 and 1034 represent H1020c and H1034b, respectively, and therefore use only those in addition to the single identifiable peak at m∕z 1048 as a second option (H-MBT; H1020c, H1034b) in Fig. 7e. Both options show a clear rise across the PETM, although the HMBT (H1020c, H1034a) shows a larger signal and somewhat better correspondence in absolute values to MBT_acyclic though with more scatter. A close correspondence between MBT_acyclic and HMBT has also been found in a lignite that has been assigned to the PETM (Inglis et al., 2019).

If the dominant source of the brGMGTs was marine throughout the record, the increase in methylation possibly relates to warming. This would not be unprecedented as marine-produced brGDGTs show an increase in methylation as a function of temperature (Dearing Crampton-Flood et al., 2018). Sollic et al. (2017) also suggest that archaeal-derived isoprenoid GMGTs produced in marine sediments incorporate additional methyl groups at higher sediment temperatures. Water column oxygen concentrations and pH also changed at our site during the PETM, which potentially affected distributions. Extensive evaluation of brGMGT distributions in modern samples is therefore required to assess the proxy potential.

5.3 Uncertainty in TEX₈₆-based SST estimates

5.3.1 Uncertainty based on calibration dataset

To calculate SSTs, we use (1) the BAYSPAR method (Tierney and Tingley, 2014), which assumes a linear relationship between TEX₈₆ and SST, and (2) ${\mathrm{TEX}}_{\mathrm{86}}^{\mathrm{H}}$ (Kim et al., 2010), which assumes a nonlinear relationship between TEX₈₆ and SST. Differences between these calibrations are smaller than the calibration errors (Fig. 6) because the TEX₈₆ values in the ACEX dataset all fall within the range of the modern core-top calibration. Taken together, both indices imply that mean annual SSTs varied between 18 and 28 ^∘C in the early Eocene, providing strong evidence for remarkable early Eocene warmth in the Arctic region.

The TEX${}_{\mathrm{86}}^{\mathrm{H}}$ calibration has a calibration error of 2.5 ^∘C (residual mean standard error; RMSE) (Kim et al., 2010). The BAYSPAR method yields possible values that range ± 6 ^∘C from the most probable value (Fig. 6), but these uncertainty estimates are more comparable than is immediately apparent, as this analysis takes a 90 % confidence interval compared to the 68 % probability of RMSE. All of the calibrations and methods to obtain values and uncertainties are based on a modern core-top dataset and thus implicitly include potential confounding factors such as seasonality and depth of production and export. However, there is no (quantitative) constraint on any of these parameters in the calibration dataset. This is particularly important for the studied region because it represents a polar endmember of the marine environment with highly seasonal production and export as well as potentially high seasonality in temperature. In the modern ocean, relations between SST and TEX₈₆ in the Arctic and ice-proximal Southern Ocean settings differ from the global ocean. This is attributed to a change in viscoelastic adaptation to temperature at the low end and/or a change in the thaumarchaeotal community (Kim et al., 2010; Ho et al., 2014; Tierney and Tingley, 2014). This may mask potential confounding factors that may be specifically relevant to polar environments. This is important here, where the polar regions were ice-free and the functioning of physical, chemical and biological ocean systems was fundamentally different from the present day. This uncertainty is not accounted for using traditional regression analyses or Bayesian techniques, and quantification of uncertainty in non-analog climates remains extremely difficult.

5.3.2 Constraints from independent proxy data

Independent proxy data may provide additional constraints. The appearance of the dinoflagellate cyst genus Apectodinium during the PETM and ETM2 in the Arctic basin (Sluijs et al., 2006; Sluijs et al., 2009; Harding et al., 2011) provides qualitative support for pronounced warming and apparent subtropical conditions. Recent efforts to quantify the paleoecological affinities of this now extinct genus have suggested a required minimum temperature of ∼ 20 ^∘C (Frieling et al., 2014; Frieling and Sluijs, 2018). Although this value is partly based on TEX₈₆ data from the ACEX cores, it is supported by data from an epicontinental site in Siberia (Frieling et al., 2014).

A second line of independent proxy evidence includes vegetation reconstructions. As indicated above, the TEX₈₆ results are qualitatively consistent with the ample evidence for thermophilic plants and animals in the Arctic (e.g., Heer, 1869; Schweitzer, 1980; Greenwood and Wing, 1995; Uhl et al., 2007; Suan et al., 2017). Particularly valuable are minimum winter temperature tolerances for specific plant species. Palynological analyses have indicated the presence of palm and baobab pollen within the PETM and ETM2 intervals in the ACEX cores (Sluijs et al., 2009; Willard et al., 2019). Modern palms are unable to tolerate sustained intervals of frost, and sexual reproduction is limited to regions where the coldest-month mean temperature (CMMT) is significantly above freezing (Van der Burgh, 1984; Greenwood and Wing, 1995). This threshold was recently quantified to be ≥ 5.2 ^∘C (Reichgelt et al., 2018). The presence of baobab within the PETM interval and ETM2 also indicates mean winter air temperatures of at least 6 ^∘C (Willard et al., 2019). Importantly, these plants were not encountered in the intervals outside the PETM and ETM2, suggesting that background coldest-month mean air temperatures were potentially too low (< 6 ^∘C) to support megathermal vegetation.

Pollen of palms and Avicennia mangroves was recently identified in time-equivalent sections in Arctic Siberia (Suan et al., 2017). Although the details of the stratigraphic framework for these records may be somewhat problematic, these findings indicate elevated CMMT estimates on land (> 5.5 ^∘C) and in the surface ocean (> 13 ^∘C) during the late Paleocene and early Eocene (Suan et al., 2017).

Apparently conflicting evidence comes from the occurrence of glendonites and erratics in specific stratigraphic levels in Paleocene and Eocene strata in Spitsbergen, interpreted to reflect “cold snaps” in climate (Spielhagen and Tripati, 2009). Some of these stratigraphic levels are very close to (or even potentially within) the PETM, considering the local stratigraphic level of the PETM (Cui et al., 2011; Harding et al., 2011). However, glendonites and erratics have not been found at the exact same stratigraphic levels as thermophilic biota (Spielhagen and Tripati, 2009). The formation and stability of ikaite (the precursor mineral of the diagenetic glendonites) in Spitsbergen were dependent on relatively low-temperature, arguably persistent near-freezing seawater temperatures in the sediment (Spielhagen and Tripati, 2009). However, glendonite occurrences in other settings (e.g., Mesozoic sediments in midlatitude regions; Teichert and Luppold, 2013) have also recently been linked to methane seeps (Morales et al., 2017). Therefore, the specific temperature constraints implied by glendonites under such conditions are the subject of debate. Future work should apply temperature reconstructions based on the geochemical composition of the glendonites and biomarkers or biota on corresponding strata to assess whether glendonite occurrence is related to colder climates.

The estimate on seasonal minima provides an important constraint on Arctic climatology during the PETM and ETM2. Most likely, the palms and baobabs grew close to the shore, where the relative heat of the ocean kept atmospheric temperatures relatively high during the winter. If minimum winter SSTs were in the range of the SST reconstructions based on the nearby Avicennia mangrove pollen (Suan et al., 2017), which for open ocean settings would perhaps amount to ∼ 10 ^∘C, then summer SST must have soared to at least 30 ^∘C in summer if TEX₈₆-based SST reconstructions of ∼ 20 ^∘C truly reflect the annual mean. It would imply an SST seasonality of ∼ 20 ^∘C, much higher than any modern open marine setting. In the present-day Arctic Ocean, heat is seasonally stored and released in sea ice melting and freezing, and sea ice cover insulates the ocean and reflects a high amount of sunlight, resulting in a seasonal cycle of not more than 1.5 ^∘C, even in ice-free regions (Chepurin and Carton, 2012). However, coupled model simulations have indicated that the future loss of sea ice will greatly enhance the seasonal SST range to up to 10 ^∘C in 2300 given unabated CO₂ emissions (Carton et al., 2015). With year-round snow- and ice-free conditions as well as even stronger summer stratification during the Eocene due to higher greenhouse gas concentrations and freshwater supply through an enhanced hydrological cycle (Pierrehumbert, 2002; Carmichael et al., 2017), a nearshore 20 ^∘C seasonal cycle in Arctic Ocean SST may not be unrealistic, although it remains inconsistent with current-generation fully coupled, relatively low-resolution model simulations (e.g., Frieling et al., 2017).

Constraints from the total pollen assemblages in the ACEX cores based on a nearest living relative approach suggest Arctic mean annual temperatures on land of 13–18 ^∘C and summer temperatures significantly exceeding 20 ^∘C during the PETM and ETM2 (Willard et al., 2019). Although these estimates come with much larger uncertainty than winter temperatures and may suffer from the non-analogous setting, they are generally lower than our TEX₈₆ values. The brGDGT-based paleothermometer MBT'_5me (De Jonge et al., 2014) also suggests lower mean annual air temperatures than reported from TEX₈₆ (Willard et al., 2019, Fig. 7). These data, derived from the same HPLC–MS analyses as the isoGDGT data presented here, indicate mean annual air temperatures averaging ∼ 18 ^∘C during the PETM, with a residual mean calibration error of 4.8 ^∘C. This value is ∼ 7 ^∘C lower than earlier estimates based on a slightly different method and analytical procedure as well as a smaller modern calibration dataset (Weijers et al., 2007a). However, recent inclusion of data from Indian soils in the MBT'_5me calibration dataset improved the proxy at the higher-temperature end (Dearing Crampton-Flood et al. 2020). The new calibration, BayMBT, results in air temperatures ∼ 3 ^∘C higher than that of De Jonge et al. (2014), which was used by Willard et al. (2019; Fig. 7). The high calibration uncertainty of MBT (see the Supplement) implies that the biomarker-derived air temperature and SST reconstructions are within error.

5.4 State of constraints on Paleocene–Eocene Arctic temperatures

To unlock the unique premise of Eocene climates for testing the skill of current-generation fully coupled climate models under high greenhouse gas forcing, proxy data and models are ideally approached separately. Among the most important implications of the Arctic temperature estimates are reconstructions of the meridional temperature gradients. Importantly, not a single simulation using an IPCC-class model of early Paleogene climate has produced Arctic annual mean sea surface temperatures close to the ACEX TEX₈₆-based reconstructions without unrealistically high tropical SSTs (Lunt et al., 2012). Recent simulations using the Community Earth System Model (CESM) versions 1 (Frieling et al., 2017; Cramwinckel et al., 2018) and 1.2 (Zhu et al., 2019) with Eocene boundary conditions produced climates that correspond to SST reconstructions in many ocean regions based on several proxies, but they still produced cooler mean annual SSTs for the Arctic Ocean than suggested by TEX₈₆ (Frieling et al., 2017; Cramwinckel et al., 2018; Zhu et al., 2019). TEX₈₆ also indicates SSTs higher than in these model simulations at several sites along the Antarctic margin (Bijl et al., 2009; Bijl et al., 2013). The question thus remains of whether the conversion of TEX₈₆ values to mean annual SST using any modern core-top calibration for high-latitude Paleogene locations is valid and whether the climate models still significantly underestimate polar temperatures. Certainly, if interpreted as mean annual SST, TEX₈₆-based estimates are high compared to the few available additional estimates, notably based on vegetation, but the latter also suffer from similar uncertainties (e.g., Hollis et al., 2019).

A few biases might lead to underestimates of meridional temperature gradients as indicated from TEX₈₆. First, the flat Eocene temperature gradient implied by TEX₈₆ was suggested to result from erroneously calibrating the proxy to SST rather than to the temperature of the subsurface (Ho and Laepple, 2016). The rationale is that the meridional temperature gradient is smaller in deeper waters than it is in the surface. However, the idea was contested for multiple reasons, including the fact that sediments at most Eocene study sites, such as the ACEX site, were deposited at a depth of less than 200m, making the application of a deep subsurface (> 1000 m) calibration inappropriate (Tierney et al., 2017). Moreover, recent analyses have indicated that the TEX₈₆ signal dominantly reflects the temperature of the top 200 m of the water column (Zhang and Liu, 2018).

Secondly, as suggested previously (Sluijs et al., 2006), if TEX₈₆ were biased towards any season in the non-analog Arctic Ocean, it would be the summer, the dominant season of organic matter export towards the seafloor through fecal pelleting or marine snow aggregates. Vegetation suggests very high winter continental coldest-month mean air temperatures of at least 6–8 ^∘C (Sluijs et al., 2009; Suan et al., 2017; Willard et al., 2019), coastal coldest-month mean SSTs of > 13 ^∘C (Suan et al., 2017), and terrestrial mean annual and warmest-month mean temperature on land of 13–21 ^∘C and > 20 ^∘C, respectively (Suan et al., 2017; Willard et al., 2019) (see Sect. 5.3.2). These estimates are closer to the most recent model simulations and lower than the existing TEX₈₆ (e.g., Frieling et al., 2017; Zhu et al., 2019). If TEX₈₆-implied SST of ∼ 25 ^∘C is skewed towards a summer estimate, this would decrease the model–data bias regarding the meridional temperature gradient estimates. Given the current uncertainties in the use of TEX₈₆ for the non-analog Arctic Ocean, we cannot independently constrain this.

6 Conclusions

We analyzed isoGDGT and brGMGT (H-shaped brGDGT) distributions in sediments recovered from the Paleocene–Eocene Thermal Maximum (PETM; ∼ 56 Ma) to Eocene Thermal Maximum 2 (ETM2; ∼ 54 Ma) interval on Lomonosov Ridge, Arctic Ocean, using state-of-the-art analytical procedures, compared them to the original dataset (Sluijs et al., 2006; Sluijs et al., 2009) and interpreted the results following the currently available TEX₈₆ proxy constraints.

Although contributions of isoGDGTs from land complicate TEX₈₆ paleothermometry in some stratigraphic intervals, temperature was the dominant variable controlling TEX₈₆ values. Background early Eocene SSTs exceed ∼ 20 ^∘C, and peak warmth occurred during the PETM and ETM2. However, uncertainty estimates of these SSTs based on the non-analog modern ocean remain complex. Temperature constraints from terrestrial vegetation support remarkable warmth in the study section and elsewhere in the Arctic basin, notably coldest-month mean temperatures around 10 ^∘C at least within the PETM and ETM2. If TEX₈₆-derived SSTs of ∼ 20 ^∘C truly represent mean annual SSTs, the seasonal range of Arctic SST might have been of the order of 20 ^∘C. If SST estimates are entirely skewed towards the summer season, seasonal ranges of the order of 10 ^∘C may be considered comparable to those simulated in future ice-free Arctic Ocean scenarios.

We find abundant brGMGTs, which appear to be predominantly produced in the marine realm at the study site. Their abundance increases during the PETM, likely due to sea level rise and perhaps due to warming and a drop in seawater oxygen concentrations. Although speculative, an increase in brGMGT methylation during the PETM may be a function of temperature, but a relation between brGMGT distribution and environmental parameters including temperature has yet to be confirmed.