Deoxygenation and organic carbon sequestration in the Tethyan realm associated with the middle Eocene climatic optimum

As a result of anthropogenic greenhouse gas emissions, Earth’s oceans are warming (IPCC, 2019) and acidifying (Doney et al., 2009), and oxygen-deficient areas are expanding (Gruber, 2011; Breitburg et al., 2018). Reconstructing and analyzing periods of greenhouse gas–driven global warmth in the geologic past can uniquely aid our understanding of these processes. In particular, the globally warm Eocene Epoch (56–34 Ma; Cramwinckel et al., 2018; Evans et al., 2018; Westerhold et al., 2020), with its superimposed periods of transient warming, holds clues to the functioning of greenhouse Earth climates (Burke et al., 2018; Hollis et al., 2019). The middle Eocene climatic optimum around 40 Ma (Bohaty et al., 2009) represents such a period of transient warming, but with many questions surrounding its causes and consequences, raising fundamental uncertainties in our understanding of global carbon cycling in the past and present climate system (Sluijs et al., 2013; Henehan et al., 2020).

The middle Eocene climatic optimum stands out as a ~400 k.y. period of transient warming in the middle Eocene (Bohaty and Zachos, 2003; Bohaty et al., 2009; Westerhold and Röhl, 2013), interrupting the long-term Eocene cooling trend (Inglis et al., 2015; Westerhold et al., 2020). Over the past decades, indications of middle Eocene climatic optimum warming have been recorded at an increasing number of localities globally. This includes deep waters in Atlantic, Pacific, and Indian deep ocean sites (Bohaty and Zachos, 2003; Bohaty et al., 2009; Dawber and Tripati, 2011; Edgar et al., 2020). Warming also occurred in surface waters of the Atlantic Ocean (Edgar et al., 2010; Boscolo Galazzo et al., 2014; Cramwinckel et al., 2018, 2020a; Arimoto et al., 2020; Henehan et al., 2020), Pacific Ocean (Bijl et al., 2010; Cramwinckel et al., 2020b; Henehan et al., 2020), and Tethys Ocean (Spofforth et al., 2010; Giorgioni et al., 2019). Reconstructed warming at these sites is on the order of 3–6 °C. Middle Eocene climatic optimum warming was likely driven by increasing atmospheric CO₂ concentrations, as indicated by reconstructions based on stable carbon isotope ratios of alkenones (Bijl et al., 2010) and boron isotope ratios (δ¹¹B) of foraminiferal calcite (Henehan et al., 2020). Global stable carbon isotope ratios (δ¹³C) of sedimentary carbonate lack a distinct coherent signature during middle Eocene climatic optimum warming, indicating that the source of excess carbon had a similar isotopic signature as that of surficial carbon reservoirs (Bohaty et al., 2009). This points to volcanic degassing (δ¹³C of ~–5‰) rather than a more ¹³C-depleted organic-derived carbon source, such as sedimentary organic carbon (C_org; δ¹³C of ~–20‰ to ~–30‰) or methane (δ¹³C ~–30‰ to -60‰; Kump and Arthur, 1999; Dickens, 2001). Although there is no evidence for large igneous province emplacement during the middle Eocene, enhanced volcanism from mid-ocean-ridge or continental arc sources has been proposed as the cause of the middle Eocene climatic optimum CO₂ rise (outlined in Bohaty et al., 2009; van der Ploeg et al., 2018). Simple carbon cycle model simulations constrained by available proxy data indicate a carbon cycle imbalance on the order of 2000–4000 Pg C (Sluijs et al., 2013; van der Ploeg et al., 2018; Henehan et al., 2020). Several tectonics-related carbon release events of roughly middle Eocene age have been identified, including increased Pacific Rim arc volcanism, metamorphic decarbonation associated with Himalayan uplift (Kerrick and Caldeira, 1999), and mid-ocean-ridge volcanism related to East African (Bailey, 1992) or North Atlantic rifting (Mjelde et al., 2008). However, none of these has a well-resolved temporal link to middle Eocene climatic optimum warming. Recently, a continental arc flare-up in Iran related to Arabia-Eurasia convergence has been proposed as a potential source (Kargaranbafghi and Neubauer, 2018; van der Boon et al., 2021), since compilation of radiometric ages of volcanic rocks in the Iran-Azerbaijan region yields a distinct peak in eruption rates around 40 Ma. There is large uncertainty (300–12,500 Pg C) on the amount of carbon released by this middle Eocene volcanism in the Iran region due to the uncertain extent of the volcanic deposits, and the potential for greatly enhanced C release due to interaction with carbonate-rich sediments (van der Boon et al., 2021). Nevertheless, the temporal correlation could indicate a causal connection between Neotethys arc volcanism and the middle Eocene climatic optimum.

The mechanism(s) responsible for subsequent climatic recovery, in the aftermath of the middle Eocene climatic optimum, presumably through atmospheric CO₂ drawdown (Bijl et al., 2010; Sluijs et al., 2013; Henehan et al., 2020), remain(s) similarly enigmatic. Remarkably, the termination of the middle Eocene climatic optimum occurred on a time scale similar to that seen at the Paleocene-Eocene thermal maximum, for which silicate weathering and organic carbon burial were likely important (Ravizza et al., 2001; Pogge von Strandmann, 2021; Papadomanolaki et al., 2022). Osmium isotope compositions of marine sediments (van der Ploeg et al., 2018) and carbon cycle simulations coupled to CO₂ reconstructions (Caves et al., 2016) indicate a weakened silicate weathering feedback during the middle Eocene compared to, e.g., the Neogene and Paleocene. Therefore, it seems difficult to conceive that silicate weathering suddenly increased in strength during the middle Eocene climatic optimum when long-term weatherability of continental rocks was low (Caves et al., 2016; van der Ploeg et al., 2018). Interestingly, the middle Eocene climatic optimum recovery interval at the Alano di Piave section in Italy (central-western Tethys; Luciani et al., 2010; Agnini et al., 2011) records elevated C_org content, indicating a potential role for C_org burial (Spofforth et al., 2010), although the lack of well-constrained C_org accumulation rates, particularly from continental margin sites, hinders quantification of this process on a global scale. Conspicuously, earlier studies also postulated a link between global cooling during the middle–late Eocene (Beniamovski et al., 2003) and Eocene-Oligocene transition (Allen and Armstrong, 2008) and deposition of organic-rich deposits in the epicontinental Peri-Tethys basin to the north of the Tethys Ocean (Sachsenhofer et al., 2018). Of these organic-rich deposits, the regionally widespread middle–late Eocene Kuma Formation can be traced laterally (Fig. 1) to the Aral Sea in the east (Beniamovski et al., 2003), Bulgaria in the west (Krasheninnikov, 1986; Sachsenhofer et al., 2018), and potentially even across the Black Sea to time-equivalent organic-rich sediments in the Ukrainian Carpathians (Hnylko and Hnylko, 2019).

Collectively, several outstanding questions related to the middle Eocene climatic optimum warrant further investigation within the Peri-Tethys region: (1) Was there a causal link between volcanism in the Iran region and middle Eocene climatic optimum greenhouse warming? (2) How were middle Eocene climatic optimum temperature change, ocean deoxygenation, and C_org burial temporally and mechanistically related? (3) Did regional enhanced C_org burial in the Peri-Tethys contribute to post–middle Eocene climatic optimum carbon drawdown?

To address these questions, we investigated the Peri-Tethyan response to the climatic changes of the middle Eocene climatic optimum in terms of regional temperature, depositional environment, and carbon drawdown. We targeted the organic-rich middle Eocene Kuma Formation in the Belaya River section (Fig. 1; van der Boon et al., 2019; Popov et al., 2019b), located on the edge of the Scythian Platform in the North Caucasus, Russia, and performed a suite of integrated organic and inorganic geochemical, nannofossil, and marine palynological analyses. To explore the potential relation between continental arc volcanism in the Iran-Azerbaijan region and the middle Eocene climatic optimum (van der Boon et al., 2021), we generated sedimentary mercury (Hg) content data—a relatively novel proxy for assessing volcanism in the geologic record that has hitherto been primarily applied to reconstruct large igneous province volcanism (Sanei et al., 2012; Jones et al., 2019), but not yet applied to relatively smaller-scale arc volcanism. To assess climate change and ocean anoxia, we employed lipid biomarker paleothermometry and redox-sensitive trace elements and iron speciation. Finally, to investigate potential carbon burial, we constructed a chronobiostratigraphic age model and integrated lithologic, carbonate, and C_org contents, and stable isotope data. We extrapolated our findings for approximate quantification of carbon drawdown throughout the Kuma Formation by using the extensive mapping of this formation in the Peri-Tethys area (Sachsenhofer et al., 2018).

Our study area was located along the Belaya River in the North Caucasus region of Russia (44.3665°N, 40.1970°E; paleolatitud ~42°N at 40 Ma according to Paleolatitude.org; van Hinsbergen et al., 2015), ~25 km south of the town of Maikop (Fig. 1). The studied section contains a succession of middle–late Eocene sediments divided into the Cherkessk, Keresta, Kuma, and Belaya Glina Formations (Zakrevskaya et al., 2011; Beniamovski, 2012), which were deposited in a continental shelf setting (van der Boon et al., 2019). Here, we defined the base of the section (0 m) as the sharp transition from the green clays of the Cherkessk Formation to the white marls of the Keresta Formation (Fig. 2). The overlying Kuma Formation (~6.8–52 m) is represented by brown organic-rich marls, while the Belaya Glina Formation (~52–112.5 m) is again characterized by white marls. Further upward, the succession continues into the Maikop Series. The Kuma Formation at the Belaya River has been dated as Bartonian in age based on biostratigraphy of planktonic foraminifera (regional Crimea-Caucasus zones PF13–PF14; Beniamovski, 2012) and calcareous nannofossils (zones NP16–NP17/CP14; Martini, 1971; Okada and Bukry, 1980; van der Boon et al., 2019; Popov et al., 2019b), in line with the age range inferred for the Kuma Formation within the broader region (ca. 44–36.5 Ma; Beniamovski et al., 2003). The Kuma Formation is generally inferred to have been deposited under low-oxygen conditions (Beniamovski et al., 2003), as the Peri-Tethys was progressively losing its connectivity to the global ocean (Palcu and Krijgsman, 2022), evolving toward the very restricted Oligocene Paratethys and associated deposition of the organic-rich Maikop Group (Sachsenhofer et al., 2017; Popov et al., 2019a). We targeted the Kuma Formation at the Belaya River section for analysis and the interval between 26 and 37 m in particular, in which the middle Eocene climatic optimum should be expressed under the assumption of continuous sediment accumulation. This interval also contained the darkest-colored and therefore presumably most C_org-rich sediments.

In total, 237 samples (sample codes BX14–BX250) were collected from the Kuma Formation (~6.8–52 m) using a gasoline-powered hand drill during field work in the summer of 2017. The average sampling resolution was chosen at ~10 cm for the interval between 26 and 37 m (suspected middle Eocene climatic optimum interval) and ~30 cm for the rest of the succession. The sample set spanning 20.85–49.55 m (156 samples; BX90–BX245) was selected for subsequent analyses, with the resolution per method chosen depending on the desired scope and available material. To ensure sufficient age-depth constraints, a wider interval (6.86–57.20 m; BX15–BX260) was selected for calcareous nannofossil analysis, encompassing the full Kuma Formation and the lowermost part of the Belaya Glina Formation.

For palynological analysis, 40 samples were processed. Between 3 and 10 g (depending on C_org content, see section “Organic Carbon, Carbonate, and Nitrogen Contents, and Organic Carbon Isotopes”) aliquots of freeze-dried, lightly crushed sediment were treated with 30% HCl and ~38%–40% HF to dissolve carbonates and silicates, respectively. The residue was sieved using 15 and 250 μm nylon mesh sieves, with ultrasonic bath steps to break up agglutinated organic matter. The resulting 15–250 μm palynomorph fraction was mounted on glass microscope slides. A general characterization of palynofacies, as well as a scan for stratigraphic marker species of dinoflagellate cysts (dinocysts), was performed with light microscopy, using the taxonomical classification of Williams et al. (2017), with the exception of wetzelielloid taxa (Bijl et al., 2017).

Calcareous nannofossil analysis was carried out on 43 samples. Samples were prepared as smear slides with standard methods (Bown and Young, 1998) and investigated with light microscopy at a magnification of 1250×. A preliminary qualitative estimation of the abundance and preservation state was performed for all the study samples. The assemblage counts were executed on a minimum of 300 specimens (modified after Thierstein et al., 1977), while for biostratigraphic aims, we counted the number of specimens of the considered taxon present in an area of 1 or 2 mm², in the latter case normalized to 1 mm² (modified after Backman and Shackleton, 1983). This extension of the counting allowed us to better define the abundance pattern of biostratigraphically important taxa with rare and/or sporadic occurrences that are often not found in standard assemblage counts. The biostratigraphic zonations were adopted from Martini (1971), Okada and Bukry (1980), Fornaciari et al. (2010), and Agnini et al. (2014). The types of biohorizons used were base (B), base common (Bc), top (T), and top common (Tc), as described in Agnini et al. (2014). The taxonomic concepts followed Aubry (1984, 1988, 1989, 1990) and Perch-Nielsen (1985), except for Dictyococcites and Cribrocentrum, for which we followed Agnini et al. (2014).

Glycerol dialkyl glycerol tetraether (GDGT) lipids were extracted from ~10 g aliquots of powdered, freeze-dried sediment (47 samples) with a Dionex accelerated solvent extractor (ASE 350), using dichloromethane (DCM):methanol (MeOH) (9:1 volume mixture) as solvent. The extracts were separated into apolar, ketone, and polar fractions over an Al₂O₃ column with elution using hexane:DCM (9:1), hexane:DCM (1:1), and DCM:MeOH (1:1), respectively. A known amount of synthetic C₄₆ glycerol trialkyl glycerol tetraether (GTGT) standard was added to the polar fraction, which was subsequently dried, redissolved in hexane:isopropanol (99:1), and filtered over a 0.45 μm polytetrafluoroethylene filter to a concentration of ~3 mg mL^–1. The GDGT-containing filtrate was analyzed using ultrahigh-performance liquid chromatography–mass spectrometry (UHPLC-MS) on an Agilent 1260 Infinity series high-performance liquid chromatography system coupled to an Agilent 6130 single-quadrupole mass spectrometer in selected ion monitoring mode at Utrecht University, Utrecht, Netherlands. Chromatographic separation of target compounds was achieved on two Waters BEH HILIC silica columns (2.1 × 150 mm, 1.7 μm) preceded by a guard column (2.1 × 5 mm, Waters) packed with the same material. Solvents, elution scheme, and instrument settings were conducted according to Hopmans et al. (2016). We scanned the resulting chromatograms for presence of isoprenoid GDGTs (isoGDGTs) and branched GDGTs (brGDGTs) (Schouten et al., 2013), but also the more recently recognized branched glycerol monoalkyl glycerol tetraethers (brGMGTs, or H-shaped brGDGTs; Liu et al., 2012; Naafs et al., 2018; Baxter et al., 2019). Tetraether index of tetraethers consisting of 86 carbons (TEX₈₆) values were calculated from isoGDGT abundances following Schouten et al. (2002) and converted to sea-surface temperature (SST) using both an exponential (Kim et al., 2010;) and a linear (O’Brien et al., 2017) calibration. Based on long-term observation of the in-house standard, the analytical precision for TEX₈₆ is ±0.3 °C. GDGT abundances (peak areas) are provided in the Supplemental Material¹, in order to facilitate recalculation of SSTs from our GDGT data using any different calibration.

For measurements of the bulk carbonate oxygen and stable carbon isotope ratios (δ¹⁸O_carb and δ¹³C_carb), powdered freeze-dried sediments were analyzed using a Thermo Finnigan GasBench-II carbonate preparation device coupled to a Thermo Finnigan Delta-V isotope ratio mass spectrometer (IRMS) at Utrecht University, aiming for 50–100 μg of carbonate per measurement. Analytical errors were ±0.1‰ for δ¹³C and δ¹⁸O based on long-term reproducibility of in-house carbonate standards. Values are reported relative to the Vienna Peedee belemnite (VPDB) standard.

In order to analyze the elemental composition and carbon isotopes of organic matter (δ¹³C_org), 0.2–0.3 g aliquots of powdered freeze-dried weighted sediment were decalcified using 1 M HCl. Residues were oven-dried at 60 °C and weighed again to obtain approximate CaCO₃ weight percentages. Total nitrogen (N) and C_org contents, as well as stable carbon isotopic ratios of C_org (δ¹³C_org), were determined using a Thermo Scientific Flask 2000 elemental analyzer coupled to a Thermo Scientific Delta V Advantage via a Conflo IV elemental analyzer–isotope ratio mass spectrometer (EA-IRMS). Analytical errors were <0.1 wt% for C_org, <0.02 wt% for N, and ~0.1‰ for δ¹³C_org based on long-term reproducibility of in-house standards. Values are reported relative to the VPDB standard.

To determine elemental composition, for each sample, ~125 mg aliquots of powdered freezedried sediment were weighed in Teflon vessels, after which 2.5 mL of 40% HF and 2.5 mL of mixed acid (HClO₄:HNO₃, 3:2) were added. Subsequently, the vessels were closed and left overnight on a hotplate at 90 °C. The next day, the lids were removed, and the temperature of the hotplate was increased to 140 °C to fume off the acids. After ~4 h, the remaining residues were dissolved in 25 mL of 4.5% HNO₃ and left overnight on the hotplate at 90 °C. The next day, the vessels were left to cool down to room temperature (~20 °C), after which they were weighed to determine the dilution. The extracts were analyzed by inductively coupled plasma–optical emission spectrometry (ICP-OES; SPECTRO ARCOS). Accuracy (recovery) was at least 95%, and average analytical uncertainty was at most 5% for the elements studied here, based on duplicates and laboratory reference material (ISE-921).

Subsamples for sedimentary Fe speciation analysis were taken during the Belaya River field campaign and stored under an oxygen-free atmosphere. In the laboratory at Utrecht University, these subsamples were freeze-dried and then powdered in an argon gas–filled glovebox. Between 50 and 80 mg aliquots of powdered sample material were weighed in centrifuge tubes and sequentially extracted (Table 1) following the method presented by Poulton and Canfield (2005). The entire extraction procedure (freeze-drying, powdering, extraction) was conducted under oxygenfree conditions. The Fe content in all extracts was determined colorimetrically using the 1,10-phenanthroline method (APHA, 2005). Average analytical uncertainty, based on duplicates, was <1 μmol/g for Fe_carb, Fe_ox1, and Fe_mag and 6 μmol/g for Fe_ox2 (Table 1). To determine the amount of Fe present in pyrite, the diffusion method given by Burton et al. (2008) was used. Between 250 and 400 mg aliquots of the anoxic, powdered and freeze-dried sediment were weighed in centrifuge tubes in an argon gas–filled glovebox. Subsequently, 10 mL of acidic chromium(II) chloride was added, while at the same time, a smaller centrifuge tube, containing 7 mL of alkaline zinc acetate solution, was placed in the centrifuge tube holding the sediment sample. The produced sulfide gas was then trapped in the alkaline zinc acetate solution, forming a precipitate of zinc sulfide. The amount of sulfur in the precipitates, i.e., chromium-reducible sulfur (CRS; generally assumed to represent pyrite, i.e., FeS₂), was then determined by iodometric titration (APHA, 2005). Average analytical uncertainty, based on duplicates, was ~20 μmol/g. The highly reactive Fe fraction (Fe_HR) was determined by summing the different fractions of the sequential Fe extraction (i.e., Fe_carb, Fe_ox1, Fe_ox2, and Fe_mag), and the Fe associated with CRS, i.e., Fe_py, was determined by dividing the amount of CRS by two.

For analysis of Hg contents, samples were prepared following the protocol described in Percival et al. (2017). Hg contents were measured using a Lumex RA-915 portable mercury analyzer paired with a PYRO-915 pyrolysis unit at the University of Oxford. Hg is reported both as concentrations and normalized to our measured C_org contents, as organic matter is generally assumed to be the most common Hg-binding sedimentary component (Sanei et al., 2012). Comparison with a high-mercury (290 ppb) paint-contaminated soil standard (National Institute of Standards and Technology NIST 2587) showed that typical analytical uncertainty was below 10%.

The palynofacies of all samples is dominated by poorly preserved, dark-colored, amorphous organic matter and algal- and plant-derived debris. Sparse palynomorphs include prasinophytes (e.g., Tasmanites and Cymatiosphaera), organic-walled dinoflagellate cysts (or dinocysts), and sporomorphs. Unfortunately, poor preservation of the material in combination with huge dilution by amorphous organic material hampered determination to the species level, and hence limited use for biostratigraphy. Because of this, we did not perform full quantitative palynological counts, but rather we performed a qualitative characterization of the assemblage. Dominant dinocyst taxa comprising a large part of the assemblages were Enneadocysta spp., Spiniferites spp., Cordosphaeridium spp., Cleistosphaeridium spp., Batiacasphaera spp., Hystrichokolpoma spp., and some Protoperidinioids and wetzelielloids. Peridinioid dinocysts were overall low in abundance relative to gonyaulacoid. The observed combination of taxa is indicative of a middle-shelf, middle Eocene setting (Pross and Brinkhuis, 2005; Sluijs et al., 2005; Frieling and Sluijs, 2018).

Since calcareous nannofossils are generally abundant throughout the studied section, showing high diversity and mostly good preservation, they formed the primary basis for our age model (Fig. 3; Table 2). Abundance patterns of index species are provided in Supplemental Figure S1 (see ^{footnote 1}). We used the biohorizons from the standard calcareous nannofossil biozonations of Martini (1971) and Okada and Bukry (1980) where possible. However, due to the low reproducibility and reliability of some of these biozonations when applied to the Belaya River section, we integrated the standard framework with the biozonation schemes of Fornaciari et al. (2010; Mediterranean Nannoplankton Paleogene [MNP] zones) and Agnini et al. (2014; Calcareous Nannofossil Eocene [CNE] zones), as listed below from old to young.

According to Backman (1987), the top of Nannotetrina spp. is marked by the top of Nannotetrina cristata and lies in the lowermost part of zone NP16 or subzone CP14a, close to the Bc (base common) of Reticulofenestra umbilicus (Agnini et al., 2014). In the Belaya River section, two specimens of Nannotetrina alata group (gr.) were observed in the lower part of the succession, but no forms ascribable to N. cristata were observed. The presence of very sporadic N. alata gr. and the absence of N. cristata suggest that the specimens observed have been reworked. This idea is also supported by the presence of other reworked taxa, and therefore we did not use this biohorizon in our age-depth model. The Bc of Cribrocentrum reticulatum is used to define the base of zone CNE14 and usually lies within chron C19r (Agnini et al., 2014). In the studied succession, C. reticulatum appears (base) from sample BX 50 (stratigraphic position 12.75 m) and increases upward in abundance (Fig. S1). The base of Sphenolithus predistentus is reported in the upper part of zone CNE14, and the abundance pattern in the study succession from sample BX 110 (stratigraphic position 26.15 m) precedes the base of D. bisectus, as previously observed in other publications (Fornaciari et al., 2010; Toffanin et al., 2011; Agnini et al., 2014). Unfortunately, the precise positioning of this biohorizon is difficult to define because of the scattered abundance of the taxon close to its appearance and the presence of transitional forms. At low to midlatitudes, the top of Sphenolithus furcatolithoides is used to define the base of subzone MNP16Ba and is consistently found in the upper part of zone CNE14 (Agnini et al., 2014) and the middle-upper part of chron C18r (Fornaciari et al., 2010). In the Belaya River section, S. furcatolithoides shows a continuous abundance pattern up to sample BX 135 (stratigraphic position 29.22 m), where it disappears (Fig. S1). The Bc of Dictyoccocites bisectus has been used as base of zone CNE15 and of subzone MNP16Bb. This taxon has been characterized by an ambiguous taxonomy (Wei and Wise, 1989); however, recent clarification (long axis >10 μm) has pointed out quite a consistent and reliable position (upper part of chron C18r) for its first occurrence (Fornaciari et al., 2010; Agnini et al., 2014). In the Belaya River section, the base of D. bisectus has been observed at sample BX 145 (stratigraphic position 30.35 m; Fig. S1). The Tc (top common) of Sphenolithus spiniger marks the base of subzone MNP16Bc (Fornaciari et al., 2010). In the study section, S. spiniger displays a continuous pattern and relatively high abundance up to sample BX 225 (stratigraphic position 41.35 m), where it decreases significantly (Tc), but the final extinction does not occur before the last investigated sample (BX 260, stratigraphic position 57.20 m; Fig. S1). The base of S. obtusus serves to define the base of subzone MNP17A, although this event is not always found to be synchronous (Fornaciari et al., 2010), possibly due to sporadic occurrences close to its first appearance (e.g., Toffanin et al., 2013). In the Belaya River section, S. obtusus is sporadic from BX 150 up to sample BX 225 (stratigraphic position 41.35 m), where it increases in abundance and is continuously present (Bc; Fig. S1). The recognition of the Bc of S. obtusus and the Tc of S. spiniger invalidates the subzone MNP17A. In order to use the biozonation of Fornaciari et al. (2010), we decided to merge subzones MNP17A and MNP16Bc.

An earlier study indicated that the Kuma Formation sediments at the Belaya River section have a very weak magnetic signal, prohibiting reliable polarity determination and thus magnetostratigraphy (van der Boon et al., 2019). As dinocyst identification was generally not possible to the species level, our age model was based primarily on calcareous nannofossil biostratigraphy and biochronology. Based on the above discussion of nannofossil datums, the tie points that constrained our age model were: the Bc of C. reticulatum, the top of S. furcatolithoides, the base of D. bisectus, and the top and Tc of S. spiniger (Fig. 3; Table 2). The base of S. obtusus was not used because of the previously discussed low reliability of this biohorizon. Biochronological estimates were recalibrated against the Geologic Time Scale 2012 (GTS12; Gradstein et al., 2012). The age of the base (42.7 Ma) and the top (39.7 Ma) of the section were extrapolated assuming linear sedimentation rates. This implies that the section encompasses ~3 m.y., partially covering the Lutetian and Bartonian and including the middle Eocene climatic optimum, which is in line with earlier work based on multiple groups of microfossils (Popov et al., 2019b). Our chronology indicates that the middle Eocene climatic optimum at the Belaya River section is complete, within the limits of the biozones. Based on average sedimentation rates in between age-depth tie points, this section is characterized by lower sedimentation rates in the older interval (6.86–30.12 m), with values of 9.3–10.3 m/m.y. Sedimentation rates are much higher in the upper part of the section, increasing almost fivefold to 45.7 m/m.y. from Bc D. bisectus (40.35 Ma) upward. In the uppermost part of the section (40.56–57.2 m), the average sedimentation rate remains high but slightly decreases to 38.4 m/m.y. This results in an average sampling resolution of ~6 k.y. for the suspected middle Eocene climatic optimum interval (sample spacing ~10 cm) and ~15 k.y. for the surrounding interval (sample spacing ~30 cm). Dinocysts present in our palynological slides support the nannofossil-based Lutetian–Bartonian age constraints, but they do not provide any additional precise age-depth information.

GDGT fractions of the Kuma Formation consist almost entirely of isoGDGTs, with brGDGTs being fully absent and brGMGTs being present in low amounts. Several indicator ratios were calculated to appraise isoGDGT sourcing. The methane index (MI; Zhang et al., 2011), GDGT-2/crenarchaeol (Weijers et al., 2011), and GDGT-0/crenarchaeol (Blaga et al., 2009), used to investigate contributions by methanotrophic and methanogenic isoGDGT producers, were all within the ordinary range for production by marine Thaumarchaeota (Supplemental Material), supporting the use of TEX₈₆ as SST proxy. Next to this, the relative abundance of crenarchaeol regio-isomer (f_{Cren′:Cren′+Cren}; O’Brien et al., 2017) and the ring index (RI; Zhang et al., 2016) were assessed to identify anomalous GDGT distributions. Two samples were discarded, as ΔRI was higher than 0.3 (Supplemental Material). Because brGDGTs were entirely below detection limit, resulting branched versus isoprenoid tetraether (BIT) index values, a measure for terrestrially derived GDGTs relative to marine-produced GDGTs (Hopmans et al., 2004; Zell et al., 2013), were 0. Notably, ratios between GDGT-2/GDGT-3 ([2]/[3]; Taylor et al., 2013) were high throughout the record, fluctuating above (12 samples) and below the threshold value of 5.0 that was suggested to indicate isoGDGT contributions by deeper-dwelling marine Archaea. TEX₈₆ values ranged between 0.76 and 0.86, recording gradual warming followed by more rapid cooling (Fig. 4). This resulted in a range of reconstructed SSTs between 30.6 °C and 34.2 °C using an exponential calibration (Kim et al., 2010; calibration standard error ±2.5 °C), and a much higher range of 33.8–39.5 °C using a linear calibration (O’Brien et al., 2017; calibration standard error ±2.0 °C).

Despite the absence of brGDGTs, structurally related brGMGTs were detected in the Belaya River section, at around 1%–2% of total GDGTs. BrGMGTs are characterized by an additional covalent bond between the alkyl chains that is thought to be an adaptation to maintain membrane stability under heat stress (Morii et al., 1998). Indeed, the abundance of brGMGTs relative to that of brGDGTs increases with temperature in peats and lake sediments (Naafs et al., 2018; Baxter et al., 2019, 2021). In the Kuma Formation, two isomers, i.e., H1020a and H1020b (nomenclature following Baxter et al., 2019), of the brGMGT labeled as H-1020 in Liu et al. (2012) and Naafs et al. (2018) were reliably detected in all samples, without a clear stratigraphic trend in the ratio between the samples.

In order to facilitate a chronological description and discussion of the results, we subdivided our record into paleoclimate regimes, primarily based on the TEX₈₆ temperature record and overall lithological changes (Figs. 4 and 5). The TEX₈₆ record was characterized by stable background temperatures in regime I transitioning to gradual warming in regime II, which we associate with the onset of the middle Eocene climatic optimum. The interval of subsequent cooling was subdivided into regimes III and IV, in which regime IV was furthermore characterized by a shift in lithology, namely, a sharp drop in carbonate content (Fig. 4). Regime V represents a return to “background” conditions, characterized by more variable temperatures. We note that there is an offset between middle Eocene climatic optimum warming in the Belaya River record as inferred from TEX₈₆ (ca. 40.7–40.3 Ma) and the global deep-ocean warming trend (ca. 40.5–40.1 Ma). However, we argue that an ~100–200 k.y. offset is perhaps expected from the current uncertainty in the biostratigraphic reference calibration (see “Calcareous Nannofossil Biozonation” subsection of the “Results” section), combined with our sampling resolution.

Lithology is relatively stable across the studied interval, consisting of bioturbated marls with high carbonate contents (~50%–80%) and relatively high C_org contents (~2%–4%; Figs. 4B–4C). The exception to this is regime IV, corresponding approximately to the darkest sedimentary facies, in which CaCO₃ drops to ~30% and C_org rises tô6%.

Recorded δ¹⁸O_carb values (Fig. 4D; Table 3) were overall rather low (–3‰ to –6‰; Fig. 4C; Table 3) compared to contemporaneous open-ocean records (Bohaty et al., 2009). Such values may suggest contributions of diagenetic calcite. However, the preservation of the calcareous nannofossils is remarkably good. These depleted values might also suggest a consistent diagenetic overprint by meteoric water (e.g., Gat, 1996), but we consider trends and variability to reflect a primary signal, modulated by changes in lithology. Values decrease from –3.3‰ to –4.0‰ in regime III, followed by a transient shift to very low values of ~–5.5‰ at ca. 40.25 Ma (regime IV). The δ¹⁸O_carb values subsequently return to –4.0‰ and gradually increase to –3.2‰ toward the top of the studied interval (regime V). By comparison, δ¹³C_carb values vary from ~1.0‰ to ~1.5‰ in regimes I–III (Fig. 4E; Table 3), followed by a transient shift to ~2.0‰–2.5‰ in regime IV. The youngest interval (regime V) is characterized by a gradual increase from ~1.5‰ tô2.0‰. The most prominent shifts, negative in δ¹⁸O_carb and positive in δ¹³C_carb, occur in regime IV, synchronous with the CaCO₃ decrease, pointing to a link with this lithological change.

In regimes I and II, molar C_org/N_tot ratios are roughly ~15 to ~20 (Fig. 4G; Table 3). Subsequently, they gradually increase to values ranging betweeñ20 and ~25 in regime V. The largest increase in values occurs throughout regimes III and IV, with a sudden step in-between the regimes. Measured δ¹³C_org varies around a value of –27.0‰ in regimes I–III (Fig. 4F; Table 3). Within regime IV, we observed a minor δ¹³C_org decrease to –27.5‰, followed by a rapid increase to values of –26.5‰. In the youngest regime V, δ¹³C_org values gradually return to –27.5‰. C_org/P_tot ratios are generally well above the Redfield ratio of 106 (Fig. 4H; Table 3). In regimes I–III, C_org/P_tot ratios are relatively stable, averaging ~160–200. At the start of regime IV, we recorded a rapid step in C_org/P_tot ratios from ~160 to ~270. Regime V retains higher values (~270) of C_org/P_tot compared to the older part of the record.

Concentrations of Mo are well above the crustal average (1–2 ppm; Taylor and McLennan, 1985) for the entire record (Fig. 5D; Table 4). Concentrations average around 40 ppm at the base of the record, after which they decrease to more stable average values of ~30 ppm throughout regimes I–II. A distinct maximum in Mo with values up to 90 ppm is recorded in regime IV. In regime V, values stabilize to around ~40 ppm. Normalizing Mo to Al (Mo/Al) dampens the peak in regime IV but still reveals a similar trend of high values at the base of the record and in regime IV. Patterns of Mo/Al are generally mirrored by other redox-sensitive trace elements such as Ni and V (Fig. 5E; Table 4; Supplemental Material). Concentrations of Mn are stable and low (generally <100 ppm) for the entire record and always far below the crustal (600 ppm; McLennan, 2001) and shale (850 ppm; Wedepohl, 1971, 1991) averages. Al-normalized Mn values are higher at the base of the record and toward the top of the record (Fig. 5C; Table 4). Between 30% and 45% of the total Fe (Fe_TOT) in the studied record is highly reactive (Fe_HR), while between 50% and 90% of this Fe_HR is Fe in pyrite (Fig. 5F; Table 3). Broadly, Fe_py/Fe_HR ratios are high (>0.8) in organic-rich regime IV. Fe_HR/Fe_TOT ratios are slightly lower (~0.32) in regime III than in regime IV (~0.35) (Fig. 5G; Table 3).

Mercury (Hg) concentrations are between 10 and 50 ppb throughout the record (Fig. 5H; Table 3) and are strongly and significantly correlated to the C_org measured in the same samples. The Hg-C_org relationship is roughly linear, as often observed in modern sediments, all else being equal (Fig. S2; Percival et al., 2021, and references therein). The linear nature of the Hg-C_org relationship allows normalization of Hg to C_org to detect Hg-cycle perturbations (Sanei et al., 2012) (Fig. 5I; Table 3). Hg/C_org ratios are low on average (6–9 ppb/%), with a total range between 2 and 14 ppb/%. While the record shows variability, Hg/C_org values lack distinctive positive anomalies or a systematic change within the phase of gradual warming (regime II).

Sediments of the Kuma Formation at the Belaya River are characterized by relatively high C_org contents (background ~3%; this study; Sachsenhofer et al., 2018), similar to other studied localities on the Russian Platform and Scythian Platform (e.g., Gavrilov et al., 2000; Beniamovski et al., 2003), suggesting high productivity and/or preservation throughout the entire interval of Kuma Formation deposition. In the palynological residues, the organic matter consists of a mix of marine palynomorphs and large amounts of algal and (higher) plant-derived debris. While fairly high ratios of C_org/N_tot (~20) might suggest pronounced contributions of terrestrial organic matter (Meyers, 1994), sedimentary C_org/N_tot values might also be elevated due to preferential loss of labile, nitrogen-rich components during early sedimentary diagenesis, as has been proposed for Cretaceous black shales (Junium and Arthur, 2007). Notably, there are no indications for high productivity in the dinocyst assemblage, as presumed heterotrophic peridinioid dinocysts are low in abundance. Taken together, high preservation is the most likely explanation for the high C_org observed here. BrGDGTs were below the detection limit in all samples, indicating a lack of erodible soils or soil mobilization in the hinterland, as well as low riverine and marine sedimentary brGDGT production and/or an enigmatic preferential transport process selecting against brGDGTs. The presence of brGMGTs despite the absence of brGDGTs suggests that these components were derived from a different source. Although a modern inventory of these components is required to fully assess their sourcing and proxy potential, H1020 was first identified in marine sediments (Liu et al., 2012) and was later shown to be abundant especially within the oxygen minimum zone of the marine water column (Xie et al., 2014). The latter study hypothesized anaerobic planktonic microbes as their main source, which might also be a possibility in our setting (see “Basin Restriction and Deoxygenation Associated with the Middle Eocene Climatic Optimum” in the “Discussion” section), analogous to the dominant marine source for brGMGT occurrences in Upper Paleocene and Lower Eocene sediments in the Arctic (Sluijs et al., 2020) and southwest Pacific Oceans (Bijl et al., 2021). In situ production in the sediments might be another possible source, as brGMGTs occur preferentially in the anoxic parts of lake sediments (Baxter et al., 2021).

We reconstructed tropically high SSTs in the Peri-Tethys throughout the middle Eocene, with background and peak SSTs of around 32 °C and 34 °C, respectively, when using an exponential (TEX₈₆^H; Kim et al., 2010) calibration, and 35 °C and 39 °C, respectively, when using a linear (Fig. 4; O’Brien et al., 2017) calibration. These temperatures are higher than coeval TEX₈₆-based temperatures from the equatorial and subtropical South Atlantic Ocean (Boscolo Galazzo et al., 2014; Cramwinckel et al., 2018). The first-order explanation for this difference is the context of the northern Peri-Tethys as an epicontinental or epeiric sea, a depositional environment that is typically a few degrees warmer than the mean ocean SST at the same latitude (e.g., Judd et al., 2020). It should be noted that, while our GDGT distributions point to marine Thaumarchaeotal production of the GDGTs comprising TEX₈₆, values of GDGT-2/ GDGT-3 are high throughout and cross the suggested threshold of 5 (Taylor et al., 2013) in the lowermost and uppermost parts of the section. Above this threshold, TEX₈₆ might additionally incorporate a signal that partially derives from deeper (>200 m) waters, possibly causing a bias in absolute reconstructed SST values. We, however, do not consider it likely that this affected our reconstructed temperature trends, as there are no large shifts in [2]/[3] ([2]/ [3] values are <5 for most of the record), and there is no correlation between TEX₈₆ and [2]/ [3] (Supplemental Material). Reconstructed high temperatures might be supported by high fractional abundance of brGMGTs relative to brGDGTs (%brGMGT), which can be used as a temperature indicator in modern peats and lakes (Naafs et al., 2018; Baxter et al., 2019). This relationship has, however, not been constrained for marine environments, and, in our record, it is very much driven by total absence of brGDGTs. Within the brGMGT distributions, the fact that H1020 isomers are present, while H1034 and H1048 are absent, would likewise point to high temperatures, if the same relationship holds as in peats (Naafs et al., 2018).

Unexpectedly, the gradual warming that pinpoints the middle Eocene climatic optimum in the Belaya River section is not mirrored by gradual changes in the other analyzed sedimentary properties, such as lithology and redoxsensitive trace elements. This indicates that regionally in the Peri-Tethys, middle Eocene climatic optimum warming was not linked to strong changes in ocean chemistry or ecosystems. Because very strong deoxygenation and ecosystem change characterize Paleocene-Eocene thermal maximum warming in this region (Dickson et al., 2014; Shcherbinina et al., 2016), we surmise that the lower rate and/or amplitude of warming during the middle Eocene climatic optimum was insufficient to result in similar regional effects. Based on global compilations of records, middle Eocene climatic optimum warming seems to have been slightly smaller in amplitude than that during the Paleocene-Eocene thermal maximum (Dunkley Jones et al., 2013; Frieling et al., 2017). More notably, the much longer duration of the middle Eocene climatic optimum warming trend (several hundreds of thousands of years) implies that the average rate of warming was one or two orders of magnitude lower (e.g., Zeebe et al., 2016).

Our Hg/C_org record from the Kuma Formation does not contain major spikes that anomalously stand out from the bulk of the values (Fig. 5I; Fig. S3). Combined Hg and C_org values of the Belaya River section overlap with those of other sections that are characterized by relatively stable and high (>2%) C_org, but that are not marked by large igneous province volcanism (e.g., Grasby et al., 2019), such as the Jurassic Kimmeridge Clay (Percival et al., 2015) and sediments from the past ~160 k.y. in the Peruvian margin upwelling zone (Fig. S2; Shen et al., 2020). Despite a suggestive temporal correlation between middle Eocene climatic optimum warming and a continental arc flare-up related to Arabia-Eurasia collision (van der Boon et al., 2021), our Hg record does not show evidence for the increased local or regional volcanism that has been hypothesized to have contributed to a carbon cycle imbalance on the order of 2000–4000 Pg C (Kargaranbafghi and Neubauer, 2018; van der Boon et al., 2021).

The Belaya River is relatively close to the region of proposed volcanism, which would promote the use of the Hg proxy in this way. On the other hand, it is not clear whether the scale of the arc flare-up perturbation to the Hg cycle would have been sufficient to register in the Hg sedimentary record. Hg contents have, over recent years, been regularly applied to track large igneous province activity in the geologic record (e.g., Sanei et al., 2012; Percival et al., 2017; Jones et al., 2019). While Hg emissions from arc volcanism are certainly significant, it should be noted that a continental arc flare-up might result in different local and global Hg signatures compared to large igneous provinces. Sedimentary Hg records spanning periods with fluctuating or elevated volcanism, but not related to large igneous provinces, are very limited in number (e.g., Roos-Barraclough et al., 2002). A recent data compilation of Hg emissions encompassing active volcanoes suggests that arc volcanoes emit more Hg per year than non-arc volcanoes (Edwards et al., 2021). However, as Hg emitted per volume of volcanic gas is indistinguishable, this mainly results from the higher overall fluxes from studied arc volcanoes (Edwards et al., 2021). Long-term mean lava discharge rates over the approximately million-year lifetimes of some studied large igneous provinces are comparable to the longterm global mean arc magma discharge rate (Mather and Schmidt, 2021), suggesting that large igneous provinces likely represent a more significant perturbation to the Hg cycle than a single regional arc flare-up. Inferring eruption activity in volcanic arcs is further complicated by the observations from ice cores, which show far-field Hg spikes are not necessarily scaled to the volume of individual eruptions, likely as a result of heterogeneous dispersal and initial Hg loading (Schuster et al., 2002).

As poorly constrained as the mechanisms governing Hg emissions from volcanoes are, if we assume the 2000–4000 Pg C perturbation was entirely released from volcanism, with a Hg/CO₂ mass ratio of 10^-8 to 10^-6 (Witt et al., 2008), this would suggest associated Hg emissions in the range of 10¹¹ to 10¹³ g. Even at the top of this range, this is significantly lower than estimates of Hg emissions from the Karoo-Ferrar large igneous province, on of order 10¹⁴ g (Percival et al., 2015), and would represent a reduced impact on the Hg cycle, and smaller likelihood of recording the perturbation, particularly given the shorter duration of the middle Eocene climatic optimum compared to a large igneous province event. Further work is required to understand whether Hg might be useful for constraining rates of non–large igneous province volcanism in the geologic record. Regarding the middle Eocene climatic optimum (and potentially similar warming events of similar duration and amplitude), we conclude that while an imbalance between volcanism and silicate weathering remains the most likely cause (e.g., van der Ploeg et al., 2018), pinpointing any exact regions of excess volcanic emissions remains challenging.

The largest changes in sediment composition (CaCO₃ decrease and C_org increase), the concomitant shifts in δ¹⁸O_carb and δ¹³C_carb, and the most prominent changes in redox-sensitive proxies all occur during the second half of middle Eocene climatic optimum cooling, in regime IV (Tables 3 and 4; Figs. 4 and 5). The increased sedimentary C_org content in this interval may relate to high marine primary productivity and/or high preservation, the latter of which may have been affected by factors such as bottom-water deoxygenation and sedimentation rate (e.g., Berner and Raiswell, 1983; Middelburg and Levin, 2009). According to our age model, sediment accumulation rates increased from the middle Eocene climatic optimum onward from ~10 to ~45 m/m.y. (Fig. 3), likely facilitating enhanced C_org preservation (Canfield, 1994; Algeo et al., 2013; Schoepfer et al., 2015). Moreover, the concentrations of redox-sensitive elements and Fe-speciation data indicate variable levels of oxygenation throughout the studied section (Fig. 5; Table 4). As brGMGTs might be produced by anaerobic planktonic microbes (Xie et al., 2014), their high abundance relative to brGDGTs might also be related to low-oxygen conditions. C_org/P_tot ratios average 221, i.e., significantly above the Redfield ratio of 106 (Fig. 4), indicating preferential regeneration of P from the sediments under anoxic bottom waters (Algeo and Ingall, 2007). Besides C_org and C_org/P_tot ratios, some sedimentary (trace) metal concentrations are sensitive to bottom-water redox conditions (e.g., Calvert and Pedersen, 1993; Tribovillard et al., 2006; van Helmond et al., 2018). In particular, sedimentary Mo has been shown to be an accurate, semiquantitative redox proxy, as sedimentary enrichment of Mo increases when sulfide concentrations in pore waters and eventually bottom waters increase (e.g., Helz et al., 1996; Algeo and Tribovillard, 2009). The high Mo concentrations recorded in the Belaya River section (Fig. 5; Table 4) are indicative of intermittently euxinic bottom waters (i.e., free sulfide in the water column) and perhaps even (more) permanently euxinic bottom waters during phases of maximum enrichment, based on comparison with the distribution of sedimentary Mo in modern sulfidic marine systems (Fig. 6; Scott and Lyons, 2012). In contrast to Mo, Mn is mobilized under reducing conditions; i.e., sediments deposited under anoxic conditions are depleted in Mn (e.g., Thomson et al., 1995; van Helmond et al., 2014). This is also the case for the entire Belaya River record (Fig. 5; Table 4), supportive of deoxygenated bottom waters. This is furthermore in agreement with the observation that benthic foraminifera are absent throughout and in the middle-upper part of the Kuma Formation at the Belaya River (Popov et al., 2019b) and Kheu River sections (Beniamovski et al., 2003), respectively, suggesting anoxic conditions.

We used the speciation of Fe between different mineral phases to further investigate the degree of deoxygenation throughout the record. Under reducing conditions, a larger fraction of the total sedimentary Fe pool is highly reactive, thus leading to higher Fe_HR/Fe_TOT ratios. Based on values in both modern and ancient sediments, a Fe_HR/Fe_TOT ratio of 0.38 was found to separate oxic-suboxic conditions from anoxic-euxinic conditions (Poulton and Canfield, 2011), although verification with other proxies is recommended (Raiswell et al., 2018, and references therein). In the studied interval, Fe_HR/Fe_TOT are close to 0.38 for the entire record (Fig. 5F), consistent with the previously discussed proxies indicating alternating to persistent reducing conditions. Intensification of reducing conditions during middle Eocene climatic optimum cooling (regime IV) is high-lighted by increases in redox-sensitive trace elements such as Mo and Ni, both in their concentrations in parts per million and when normalized to Al (Fig. 5; Table 4). In addition to Fe_HR/Fe_TOT ratios, the degree of Fe pyritization, as expressed in the Fe_py/Fe_HR ratio, can be used to reconstruct bottom-water redox conditions. In regime IV, this ratio exceeds 0.8 (Fig. 5G), indicating euxinic conditions were likely prevalent in the bottom waters (Anderson and Raiswell, 2004). Increased preservation of organic matter under low-oxygen conditions thus likely facilitated the sharp increase in C_org that is recorded in the section, in addition to the high sedimentation rates previously mentioned.

As this interval of intensified deoxygenation occurred during cooling following regional peak middle Eocene climatic optimum warmth, euxinia in the northern Peri-Tethys was not directly caused by ocean warming or associated sea-level rise. Interestingly, the slope of the Mo-C_org relationship decreases during regime IV (Fig. 6). While most of the Belaya River samples fall on the regression line of the modern Saanich Inlet, the data from the C_org-rich interval overlap with the (more restricted) Cariaco Basin. This change could relate to reduced availability of aqueous Mo due to increased hydrographic restriction of the deep waters (Algeo and Lyons, 2006), pointing to a transient increase in basin restriction as a possible cause for the transient changes in oxygenation and sedimentary chemistry during middle Eocene climatic optimum cooling. Part of the sharp decline in δ¹⁸O_carb in regime IV, coincident with Mo/C_org-based restriction, might therefore be explained by reduced surface-water δ¹⁸O, i.e., reduced salinity, caused by enhanced relative influence of freshwater input over marine connections. The interval of euxinia might therefore not have been directly tied to middle Eocene climatic optimum climate change; instead, it was likely a consequence of intensified ongoing tectonic restriction of the Peri-Tethys (Palcu and Krijgs man, 2022). Specifically, the intensified arc volcanism in the Iran-Anatolian region around middle Eocene climatic optimum times (Verdel et al., 2011; van der Boon et al., 2021) might well represent a step change in restriction of the southern gateways connecting the Peri-Tethys to the Indian Ocean (Palcu and Krijgsman, 2022), as regional mountain belts were uplifted (Keskin et al., 2008). The Peri-Tethys–Siberian (Turgay strait) and Siberian–Arctic (Kara strait) connections likewise show increased restriction around middle Eocene climatic optimum times (Palcu and Krijgsman, 2022). Interestingly, a shift to more arid conditions is recorded in dust records from central China at 40 Ma (Bosboom et al., 2014; Meijer et al., 2021), i.e., coincident with the middle Eocene climatic optimum, which Meijer et al. (2021) related to a major Peri-Tethys regression from the Tarim and Tajik basins around 41–40 Ma (Kaya et al., 2019). We surmise that increased restriction and euxinia in the Kuma basin can be mechanistically related to the observed middle Eocene climatic optimum–associated cooling if boreal surface inflow, e.g., from the North Sea, became dominant over the southern gateways, thereby suppressing middle Eocene climatic optimum warmth (Palcu and Krijgsman, 2022). Model simulations with sophisticated oceanography and bathymetry would be necessary to assess scenarios restricting different tropical and boreal gateways in the middle Eocene and their effects on Peri-Tethys SSTs and oxygenation. Nevertheless, this study highlights the importance of the interplay between global climate and regional gateway evolution in determining local climate and oceanographic change.

According to our age constraints, sedimentation rates at the Belaya River section increased almost fivefold during the middle Eocene climatic optimum, while C_org contents during middle Eocene climatic optimum cooling (regime IV) record an increase from ~3 wt% to ~5.5 wt% over an interval of ~55 k.y. (Fig. 4). Using an estimated dry bulk sediment density of 2.5 g cm^–3 together with the above values results in C_org mass accumulation rates of ~0.6 g C cm^–2 k.y.^–1 for this interval. The Eocene Peri-Tethys was prone to developing (intensified) low-oxygen conditions, as evidenced both by the high organic content of the Kuma Formation in general (Beniamovski et al., 2003) and also by a similar deoxygenation response to Paleocene-Eocene thermal maximum warming ca. 56 Ma (Dickson et al., 2014) and potentially during later hyperthermal events (Radionova et al., 2003), and during the Eocene-Oligocene transition (Sachsenhofer et al., 2017). In order to provide a ballpark figure of total C_org burial in the northern Peri-Tethys during middle Eocene climatic optimum cooling, we estimated the areal extent of the Kuma Formation to b ~1.5 × 10⁶ km², based on the area maps of Sachsenhofer et al. (2018). This is a conservative estimate, as organic-rich lateral equivalents of the Kuma Formation might extend further, such as past the Black Sea to the west (Hnylko and Hnylko, 2019). Using our values at the Belaya River section, we extrapolated the total C_org burial during middle Eocene climatic optimum cooling in the northern Peri-Tethys to be ~9.2 Tg C yr^–1, ~3% of the modern global C_org burial of ~300 Tg C yr^–1 (Supplemental Material). This results iñ500 Pg buried C total over the whole 55 k.y. interval. Expressing this as excess burial relative to pre–middle Eocene climatic optimum values, and extrapolating from the Belaya River C_org and sedimentation rates, yields excess C_org burial of ~8.1 Tg C yr^–1 and a total of ~450 Pg excess C_org buried. Such an amount might have played a significant role in middle Eocene climatic optimum CO₂ drawdown, but it was likely not of sufficient magnitude to explain the full recovery to a pre-event background global climate state. Yet, there are indications that middle Eocene climatic optimum–associated enhanced C_org burial might not have been limited to the Peri-Tethys region. A similar pattern of middle Eocene climatic optimum warming followed by enhanced sedimentary C_org content has been observed in sediments from the Alano di Piave section in Italy, central western Tethys (Fig. 1; Spofforth et al., 2010). At Alano, two δ¹³C-enriched, C_org-rich (~2 wt%) layers with a thickness of several meters occur in the post–middle Eocene climatic optimum interval. These layers are concurrent with planktic and benthic foraminiferal assemblages suggesting eutrophic sea-surface conditions, high export productivity, and oxygen-depleted bottom waters right after the zenith of middle Eocene climatic optimum warmth (Luciani et al., 2010; Boscolo Galazzo et al., 2013). This suggests that eutrophic and anoxic conditions facilitated production, preservation, and ultimately burial of C_org. The Monte Cagnero section in Italy records a similar δ¹³C-enriched post–middle Eocene climatic optimum layer with increased magnetic susceptibility, possibly likewise indicating C_org enrichment (Savian et al., 2014; Kochhann et al., 2021). The middle Eocene climatic optimum has also been recognized in other Tethyan sites, such as the Contessa Highway section, Italy (Jovane et al., 2007), the Elazig Basin, Turkey (Rodelli et al., 2017; Rego et al., 2018; Giorgioni et al., 2019), and the oil shales of Jordan (Hussein et al., 2014), but well-constrained C_org accumulation rates from multiple sites within and outside of the Tethys Ocean are lacking. For the present study, we therefore refrained from extrapolating C_org burial to outside of the Peri-Tethyan Kuma Formation. Still, we note that a subtle positive δ¹³C trend and maximum after the peak δ¹⁸O minimum in the global benthic stable isotope compilation (Bohaty et al., 2009; Westerhold et al., 2020) are consistent with a role for sequestration of relatively ¹²C-rich C_org from the exogenic carbon pool, resulting in atmospheric CO₂ drawdown and climatic recovery following the middle Eocene climatic optimum. Estimating the amount of global C_org burial during middle Eocene climatic optimum recovery and modeling its consequences for climate and carbon cycling represent promising targets for future studies in order to improve our understanding of climate-carbon feedbacks during events on similar timescales to the middle Eocene climatic optimum.

The middle Eocene climatic optimum (ca. 40 Ma) was a period of transient global warming lasting ~400 k.y. that interrupted long-term Eocene cooling. The controls on its initiation and termination are poorly understood but might offer key insights into Earth’s Eocene climate system. In this study, we targeted the organic-rich middle Eocene Kuma Formation in the Belaya River section, North Caucasus, in order to investigate the Peri-Tethyan response to the climatic changes of the middle Eocene climatic optimum in terms of regional temperature, depositional environment, and carbon drawdown. We found that the middle Eocene Peri-Tethys was characterized by high TEX₈₆-based SSTs, on which middle Eocene climatic optimum warming stands out as an additional 2–3 °C rise (Fig. 4). A nearby contemporaneous arc-volcanism flare-up, evidenced by deposits left in modern-day Iran and hypothesized to play a role in driving environmental changes across the middle Eocene climatic optimum, left no clear trace in the sedimentary Hg record in the Belaya River section. However, it is currently uncertain whether volcanism on this scale is expected to perturb the global Hg cycle sufficiently to be recorded in the geologic record. Intriguingly, we found that the most striking changes in sediment lithology and geochemistry did not occur during middle Eocene climatic optimum warming or peak warmth, but were concomitant with subsequent cooling. Based on sedimentary trace-element compositions and iron speciation (Fig. 5), middle Eocene bottom waters in the study area were consistently reducing, with an alternation of suboxic, anoxic, and euxinic conditions. The most reducing conditions, with presumably euxinic bottom waters, occurred during the middle Eocene climatic optimum cooling phase, which was furthermore characterized by enhanced C_org burial. Extrapolating from our record, excess C_org burial might have been ~8.1 Tg C yr^–1, resulting in a total ~450 Pg of excess C_org buried over middle Eocene climatic optimum cooling in the northern Peri-Tethys. As indications for enhanced C_org burial also exist in the western Tethys, this phenomenon might have played a quantitatively significant role in middle Eocene climatic optimum recovery. This illustrates that the closing Tethys Ocean might have affected global Paleogene climate by functioning as a large organic carbon sink.

This work was carried out under the program of the Netherlands Earth System Science Centre (NESSC), financially supported by the Dutch Ministry of Education, Culture and Science. We thank Sergey Popov and Larisa Golovina (Russian Academy of Sciences), Michael Morton and Sarah Davies (University of Leicester), Arjen Grothe and Kevin Vis (Utrecht University), and Stephen Vincent (Cambridge Arctic Shelf Programme, CASP, University of Cambridge) for field support and discussions. We thank Kirsten de Haan, Arnold van Dijk, Coen Mulder, and Natasja Welters (Utrecht University) for analytical support. T.A. Mather and J. Frieling acknowledge funding from a European Research Council (ERC) consolidator grant (ERC-2018-COG-818717-V-ECHO), A. Sluijs acknowledges ERC consolidator grant 771497 (SPANC), and C.P. Slomp and N.A.G.M. van Helmond acknowledge funding from Dutch Research Council (NWO) Vici grant 865.13.005. M.J. Cramwinckel and A. Sluijs thank the Ammodo Foundation for funding unfettered research by laureate A. Sluijs.