Impact of deoxygenation and hydrological changes on the Black Sea nitrogen cycle during the Last Deglaciation and Holocene

The marine nitrogen (N) cycle profoundly impacts global ocean productivity. Amid rising deoxygenation in marine environments due to anthropogenic pressures, understanding the impact of this on the marine N cycle is vital. The Black Sea's evolution from an oxygenated lacustrine basin to an anoxic marine environment over the Last Deglaciation and Holocene offers insight into these dynamics. Here, we generated records of the organic biomarkers heterocyte glycolipids (HGs), crenarchaeol, and bacteriohopanetetrol, associated with various water column microbial N-cycle processes, which indicate a profound change in Black Sea N-cycle dynamics at ∼7.2 ka when waters became severely deoxygenated. This transition substantially reduced Thaumarchaeota-driven nitrification and enhanced loss of bioavailable nitrogen through anaerobic ammonium oxidation (anammox). In contrast, other climatic changes over the Last Deglaciation and Holocene, such as freshwater input, water-level variations, and temperature changes, did not impact these processes. Cyanobacterial nitrogen fixation in surface waters proved more responsive to changes in salinity, which affected species composition, and associated water column stratification, which reduced the vertical transport of nutrients. Our results indicate that future deoxygenation in certain marine environments may enhance bioavailable nitrogen loss by anammox and reduce nitrification by Thaumarchaeota, while enhanced stratification may increase cyanobacterial nitrogen fixation in the surface waters.

1 Introduction

The marine nitrogen (N) cycle is a significant control of biological productivity in the global oceans. It is directly connected to the fixation of atmospheric carbon dioxide (CO₂) and carbon export from the ocean's surface, influencing atmospheric CO₂ levels over geological timescales (Falkowski et al., 1998). As the marine N cycle is strongly regulated by biology, the (de)oxygenation of the ocean determines the microorganisms involved in these biogeochemical cycles and the aerobic/anaerobic pathways that occur. Under anoxic conditions, the loss of bioavailable nitrogen is substantial, attributed to anaerobic ammonium oxidation (anammox) and denitrification (Kuypers et al., 2003; Dalsgaard et al., 2012). With deoxygenation in marine environments increasing due to anthropogenic climate and environmental changes (i.e. Keeling et al., 2010; Bopp et al., 2013) and with research linking deoxygenation to changes in the marine N cycle (Kalvelage et al., 2013; Naafs et al., 2019), it is important to enhance our understanding of how the marine N cycle may respond to future deoxygenation and what the associated feedbacks on carbon fixation might be.

Marine basins that have experienced changes in oxygenation in the past can provide perspective on the current deoxygenation of modern global oceans and the associated feedbacks in the marine N cycle, particularly on timescales beyond the observational record. Today, the Black Sea is the world's largest permanently stratified anoxic basin, with limited connection to the global ocean through the Bosporus Strait, and its redox gradient is a hotspot of diverse microbial populations and metabolisms, experiencing many of the same crucial microbial biogeochemical cycle processes as the global ocean (Kusch et al., 2022). However, over the Last Deglaciation and Holocene (approximately the last 20 ka), the Black Sea has experienced large hydrological changes. The basin was an oxygenated freshwater lacustrine environment during the Last Glacial Maximum (LGM) (Schrader, 1979) and experienced many environmental changes during the subsequent deglaciation, including temperature changes (Bahr et al., 2005, 2008; Ion et al., 2022), water-level variations (Ivanova et al., 2007; Nicholas et al., 2011; Piper and Calvert, 2011), and changes in freshwater input into the basin, through the melting of both Eurasian ice sheets and alpine glaciers after the LGM and changes in regional precipitation (Bahr et al., 2005, 2006, 2008; Badertscher et al., 2011; Shumilovskikh et al., 2012). It became reconnected to the global ocean at ∼9.6 ka when post-glacial sea-level rise caused an initial marine inflow (IMI) over the Bosporus sill (Aksu et al., 2002; Major et al., 2006; Bahr et al., 2008; Ankindinova et al., 2019), leading to enhanced salinity of the water column (Marret et al., 2009; Verleye et al., 2009; Filipova-Marinova et al., 2013) and euxinic deep waters developing in the basin after ∼7.2 ka (Arthur and Dean, 1998; Eckert et al., 2013). Thus, sedimentary records of the Black Sea may provide a unique perspective of the impacts of deoxygenation, changing temperature, and salinity on the marine N cycle.

Diagnostic lipid biomarkers of microbes preserved in the geological record can offer a unique insight into past changes in the N cycle (Rush and Sinninghe Damsté, 2017, and references cited therein; Elling et al., 2021; van Kemenade et al., 2023), although the limitations of extrapolating modern findings to ancient climates must be acknowledged, as past ecosystems may have operated under different dynamics that are not fully captured by contemporary analogues. Nitrogen-fixing heterocytous cyanobacteria play a crucial role in transforming nitrogen gas (N₂) into bioavailable nitrogen (NH₃) and sustaining primary productivity in both marine and freshwater environments (Villareal, 1992; Ploug, 2008). Identification of their diagnostic biomarkers, heterocyte glycolipids (HGs), in the geological record are a widely used proxy for exploring past changes in nitrogen fixation by these microbes (Bauersachs et al., 2009, 2010; Sollai et al., 2017; Bale et al., 2019; Elling et al., 2021; Pérez Gallego et al., 2025). The structure of the sugar moiety of HGs is a useful indicator of palaeoenvironmental conditions, as HGs with a hexose (C₆) headgroup are typically found in free-living heterocystous cyanobacteria (Bauersachs et al., 2009; Wörmer et al., 2012), while HGs with a pentose (C₅) headgroup are associated with endosymbiotic heterocystous cyanobacteria in marine diatoms (diatom–diazotroph associations, DDAs) (Schouten et al., 2013; Bale et al., 2015). Consequently, both hexose and pentose HGs have been applied as specific palaeo-biomarkers for the presence of N₂-fixing cyanobacteria in marine and lacustrine geological records (Bauersachs et al., 2010; Sollai et al., 2017; Bale et al., 2019; Elling et al., 2021). Nitrification, the microbial two-step conversion of ammonia (NH₃) and/or ammonium (${{\mathrm{NH}}_{\mathrm{4}}}^{+}$) to nitrate (${{\mathrm{NO}}_{\mathrm{3}}}^{-}$), is a central part of the marine N cycle. Archaea of the phylum Thaumarchaeota (also known as Nitrososphaerota) are among the most abundant and widespread marine prokaryotes (Karner et al., 2001; Francis et al., 2005), playing a crucial role in nitrification in the Black Sea (Lam et al., 2007) by aerobically oxidizing ammonia to nitrite (Könneke et al., 2005; Wuchter et al., 2006). As Thaumarchaeota are the exclusive producers of the membrane-spanning lipid, crenarchaeol (Sinninghe Damste et al., 2002), this biomarker can be used to identify Thaumarchaeota in the geological record and explore the palaeo-marine N cycle. Another critical part of the N cycle is the loss of bioavailable nitrogen to N₂. Under anoxic conditions, bioavailable nitrogen (${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, NH₃, and ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$) can be lost through two processes in subsurface waters: anammox (van de Graaf et al., 1997; Kuypers et al., 2003) and denitrification (Kuenen and Robertson, 1988). It is possible to explore past changes in anammox activity in the sedimentary record using the unique ladderane fatty acids (Sinninghe Damste et al., 2002), but these are relatively poorly preserved in sediments (Jaeschke et al., 2007). Alternatively, the ratio of bacteriohopanetetrol (BHT)-34S (which is ubiquitously synthesized by aerobic bacteria) to the later-eluting stereoisomer BHT-x (which is predominately synthesized by marine anammox bacteria, i.e. Ca. Scalindua spp.) (Rush et al., 2014; Schwartz-Narbonne et al., 2020; van Kemenade et al., 2023) can be used to trace past anammox activity. Denitrification is performed by a large range of organisms (Knowles, 1982), but, at present, there are no associated diagnostic lipid biomarkers (Rush and Sinninghe Damsté, 2017).

In this study, we used lipid biomarkers of microbes widely involved in the N cycle across various marine and freshwater environments in combination with other geochemical records from a sediment core located in the western Black Sea spanning the Last Deglaciation and Holocene (∼20 ka–present) to better constrain and assess the sensitivity of the marine N cycle under changing hydrological and oxygenation conditions and explore its potential links to broader global climate dynamics.

2 Regional setting

The Black Sea is a large meromictic marginal basin connected to the Mediterranean Sea via the Turkish Straits (the Bosporus Strait, the Sea of Marmara, and the Dardanelles Strait) (Fig. 1). The Black Sea has a net outflow into the Aegean Sea via the Turkish Straits and is primarily supplied by three major rivers: the Danube, the Dnieper, and the Don. With freshwater flowing out of the basin and dense, highly saline waters flowing in, the water column is highly stratified with respect to salinity (density). An oxygenated colder surface layer (0–50 m) overlies warmer, anoxic, sulfidic, hypersaline deep waters (100–2300 m), separated by a suboxic layer (50–100 m) (Murray et al., 1989, 1995). The general circulation of Black Sea surface waters is a basin-scale cyclonic boundary current encompassing large eastern and western cyclonic gyres, with several smaller, anticyclonic coastal eddies (Fig. 1) (Özsoy and Ünlüata, 1997). The modern Black Sea water column is characterized by the presence of strong redox gradients. In the basin's western gyre, the oxic zone (0–75 m depth range) has an oxygen concentration of ∼121 µmol kg⁻¹ and salinity of 19.4 psu at 50 m depth, with no detectable sulfide concentrations (Sollai et al., 2018; Bale et al., 2021). The suboxic zone lies below (75–115 m depth range), where salinity increases with depth, and traces of sulfide are detected at the bottom of this layer at ∼110 m (Sollai et al., 2018; Bale et al., 2021). Beneath this lies the euxinic zone (115–2000 m depth range), where salinity and sulfide concentrations increase with depth, reaching 22.3 psu and ∼400 µmol L⁻¹ at a depth of 2000 m (Sollai et al., 2018; Bale et al., 2021). Ammonium concentrations are low in the oxic and upper suboxic zones (<0.1 µmol L⁻¹), increasing below 90 m, with highest concentrations (∼100 µmol L⁻¹) at 2000 m water depth (Sollai et al., 2018). Nitrite concentrations are highest in the oxic zone, peaking at 50 m (∼0.08 µmol L⁻¹), decreasing through the suboxic zone, with the exception of a peak at 85 m, while, in the euxinic zone, nitrite is below the limit of detection (Sollai et al., 2018). Nitrate concentrations are highest in the lower oxic and upper suboxic zones, peaking between 70–80 m (∼2.5 µmol L⁻¹), decreasing with depth, and reaching the limit of detection below 105 m (Sollai et al., 2018; Bale et al., 2021).

Figure 1 Map of the Black Sea basin, showing the major surface circulation, the depth contours (m b.s.l.), and the location of core 64PE418 (adapted from Giorgi Balakhadze, English Wikipedia, 2016).

3 Methods

During the cruise with the RV Pelagia in April 2017, piston core 64PE418 (235 cm length) was recovered from 1970 m below sea level (m b.s.l.) depth in the Black Sea (42°56 N, 30°02 E) (Fig. 1). A total of 44 sediment samples were taken at 5 cm intervals along the depth of the core.

3.1 Biomarker extraction and analysis

Lipids were extracted from these samples using a modified Bligh and Dyer extraction method as described previously (Bale et al., 2021). Using a mixture of methanol (MeOH), dichloromethane (DCM), and phosphate buffer ($\mathrm{2}:\mathrm{1}:\mathrm{0.8}$, v:v), the sediment was twice extracted ultrasonically (10 min). The combined supernatants were phase-separated by adding DCM and phosphate buffer to create a solvent ratio of $\mathrm{1}:\mathrm{1}:\mathrm{0.9}$ (v:v). The organic phase was collected, and the aqueous phase was re-extracted three times using DCM. All extraction steps were then repeated on the residue but with a mixture of MeOH, DCM, and aqueous trichloroacetic acid solution (TCA) pH 3 ($\mathrm{2}:\mathrm{1}:\mathrm{0.8}$, v:v). Finally, the organic extracts were combined and dried under an N₂ gas stream. A deuterated betaine lipid {1,2-dipalmitoyl-sn-glycero-3-O-4'-[N,N,N-trimethyl(d9)]-homoserine; Avanti Lipids} internal standard was added to each sample before filtering the extract through 0.45 µm cellulose syringe filters (4 mm diameter; BGB, USA). Extraction blanks were performed alongside the sediment extractions, using the same glassware, solvents, and extraction methodology, but without sediment. Analysis of the extracts was performed using the following UHPLC-HRMS reversed-phase method. An Agilent 1290 Infinity I UHPLC instrument was used, equipped with a thermostatted auto-injector and column oven, and coupled to a Q Exactive Orbitrap MS instrument with an Ion Max source using a heated electrospray ionization (HESI) probe (Thermo Fisher Scientific, Waltham, MA). Separation was achieved using an Acquity BEH C18 column (Waters, 2.1 mm×150 mm, 1.7 µm) maintained at 30 °C. The eluent composition was (A) MeOH/H₂O/formic acid/14.8 M NH₃aq [$\mathrm{85}:\mathrm{15}:\mathrm{0.12}:\mathrm{0.04}$ (v:v)] and (B) IPA/MeOH/formic acid/14.8 M NH₃aq [$\mathrm{50}:\mathrm{50}:\mathrm{0.12}:\mathrm{0.04}$ (v:v)]. The elution programme was 95 % A (for 3 min) followed by a linear gradient to 40 % A (at 12 min) and then to 0 % A (at 50 min), which was maintained until 80 min. The flow rate was 0.2 mL min⁻¹. Positive-ion HESI settings were as follows: capillary temperature, 300 °C; sheath gas (N₂) pressure, 40 arbitrary units (AU); auxiliary gas (N₂) pressure, 10 AU; spray voltage, 4.5 kV; probe heater temperature, 50 °C; S lens, 70 V. Lipids were analysed with a mass range of $m/z$ 350–2000 (resolving power 70 000 ppm at $m/z$ 200), followed by data-dependent tandem MS/MS (resolving power 17 500 ppm), in which the 10 most abundant masses in the mass spectrum were fragmented successively. Optimal fragmentation was achieved with a stepped normalized collision energy of 15, 22.5, and 30 (isolation width, 1.0 $m/z$) for IPL analysis (Bale et al., 2021) and 22.5 and 40 (isolation width 1.0 $m/z$) for BHP analysis (Hopmans et al., 2021). The Q Exactive was calibrated within a mass accuracy range of 1 ppm using the Thermo Scientific Pierce LTQ Velos ESI Positive Ion Calibration Solution. During analysis, dynamic exclusion was used to temporarily exclude masses (for 6 s) to allow the selection of less abundant ions for MS/MS.

Biomarkers were identified based on their retention time, exact mass, and fragmentation spectra. Integrations were performed on (summed) mass chromatograms of relevant molecular ions ($[\mathrm{M}+\mathrm{H}{]}^{+}$, $[\mathrm{M}+{\mathrm{NH}}_{\mathrm{4}}{]}^{+}$, and $[\mathrm{M}+\mathrm{Na}{]}^{+}$) and, in the case of crenarchaeol, also the second isotope peak for each of the three adducts. Due to the coelution of BHT-34S, BHT-x isomer, and an unknown nitrogen containing a compound with the same mass, identification and integration of BHT-34S and BHT-x were conducted using the $m/z$ 529.462 dehydrated in-source product ($[\mathrm{M}+\mathrm{H}{]}^{+}$-H₂O). Isoprenoidal glycerol dialkyl glycerol tetraether (isoGDGT) crenarchaeol, monohexose crenarchaeol, and a crenarchaeol isomer were all integrated and combined as “crenarchaeol”. The lipid biomarker records are all presented as peak area per gram of total organic carbon (TOC).

3.2 Total organic carbon, total nitrogen, and δ¹⁵N_bulk measurements

Freeze-dried sediments were analysed for TOC, total nitrogen (TN), and bulk δ¹⁵N (δ¹⁵N_bulk) using a Thermo Scientific Flash EA DELTA V Plus IRMS instrument. Flow was 100 mL min⁻¹, and the temperatures for oxidation, reduction, and the oven were 900, 680, and 40 °C, respectively. Nitrogen isotopic measurements were calibrated to atmospheric air (AIR), and values are expressed in per mille (‰) units. Inorganic carbon was removed from the sediment prior to TOC analysis using HCl (2 mol), cleaned with bi-distilled water, then freeze-dried.

3.3 Age model

Accelerator mass spectrometry (AMS) ¹⁴C ages of bulk organic matter were measured from core 64PE418 (n=7) to create a chronology on the 64PE418 depth scale. Samples were weighed and freeze-dried at NIOZ. The AMS ¹⁴C measurements (${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$) were determined using a Compact Carbon AMS system at the Poznań Radiocarbon Laboratory, Poland. The sediment samples were pre-treated with 0.25 M HCl (room temperature overnight, then 80 °C, 1+ h) and rinsed with deionized water until reaching a pH of 7. Samples were then combusted in closed (sealed under vacuum) quartz tubes, together with CuO and Ag wool (900 °C, 10 h). The CO₂ released was then dried in a vacuum line and reduced with H₂ using 2 mg of iron (Fe) powder as a catalyst. The obtained carbon and Fe mixture was then pressed into an aluminium holder (Czernik and Goslar, 2001). The measurement was performed by comparing intensities of ionic beams of ¹⁴C, ¹³C, and ¹²C measured for each sample and for standard samples (with “Oxalic Acid II” used as modern standard and “coal” used as background standard of ¹⁴C-free carbon). In each AMS run, 30–33 samples of unknown age were measured, alternated with measurements of 3–4 samples of modern standard and 1–2 samples of background standard. The measured ${}^{\mathrm{14}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$ ratios were corrected for isotopic fractionation and reported as conventional radiocarbon age according to Stuiver and Polach (1977).

Seven bulk organic matter ¹⁴C dates were used in the production of the age model for core 64PE418 (Table S1 and Fig. S1 in the Supplement). Six of these were from this core, with an additional bulk organic carbon ¹⁴C date from the widely acknowledged Unit I/II boundary of core KNR 134-08 BC17, which was used to further refine the age model for the upper part of the core (Jones and Gagnon, 1994). Core KNR 134-08 BC17 was sourced from the same location and water depth as 64PE418, and this boundary was identified in our core using the same significant colour and elemental changes described in previous studies (Figs. 2 and 3) (i.e. Arthur and Dean, 1998; Bahr et al., 2005). While seven ¹⁴C measurements were conducted on core 64PE418, one was excluded from the age model due to an age reversal (142.5 cm), likely due to the presence of reworked material, which would have made the age model out of line with other studies (Fig. S2 in the Supplement). Variable reservoir ages were added to our calibration (Table S1), using those calculated by Kwiecien et al. (2008) for intermediate water depths in the Black Sea over the Last Deglaciation and Holocene which were deemed the most suitable for our site location. The ¹⁴C dates were calibrated using the Marine20 calibration curve (Heaton et al., 2020) for the upper three samples (24.5, 39, 76.5 cm), which reflect the period after the infiltration of marine water; this is based on the colour and elemental changes in the core which indicate that these samples fall within Units I and II (Arthur and Dean, 1998; Bahr et al., 2005). The lower four samples (118.5, 158.5, 183.5, and 217.5 cm) were calibrated using the IntCal20 calibration curve (Reimer et al., 2020), as they reflect the period prior to the marine infiltration when the Black Sea was a lacustrine environment, as indicated by colour and elemental signatures in the core (Arthur and Dean, 1998; Bahr et al., 2005). Using the R code CLAM (Blaauw, 2010), the age–depth model was created based on the seven ¹⁴C dates. Our age model shows that the 64PE418 biomarker records span the last 19.5 ka, with an average resolution of ∼450 years. The following transitions are identified in our core by colour (Fig. 2) and elemental changes (Fig. 3) and are dated by our age model as follows: the onset of the IMI (138 cm) is at 9.6 ka ± 237 yr, the boundary of Unit II/III (96 cm) is dated at 7.2 ka ± 202 yr, and the Unit I/II boundary (39 cm) is dated at 2.6 ka ± 402 yr. Despite the complexity of dating Black Sea cores, due to lack of suitable material for dating and dispute over reservoir age corrections (Kwiecien et al., 2008; Soulet et al., 2011; Yanchilina et al., 2017), the dates of these boundaries align well with previously published calibrated ages for these transitions (i.e. Jones and Gagnon, 1994; Ankindinova et al., 2019; Huang et al., 2021), as shown in Fig. S2.

Figure 2 Scan of core 64PE418 showing colour changes and the depth of unit boundaries. Unit boundaries are defined according to Arthur and Dean (1998) and are identified by colour changes and XRF core scan changes in Ti and Ca (Fig. S1).

Figure 3 Changes over the last 19.6 ka in the elemental content of core 64PE418, as measured through calibrated XRF core scanning (ppm) using the methods described in Hennekam et al. (2020).

4 Results

4.1 TOC, TN, and colour changes

Sedimentary bulk TOC (%), bulk TN (%), and δ¹⁵N_bulk (‰) range between 0.3 %–22.8 % for TOC, 0.05 %–1.9 % for TN, and 5.2 ‰–0.0 ‰ for δ¹⁵N_bulk (Fig. 4). There are significant colour changes in the core, as shown in Fig. 2, which correspond to changes in TOC, TN, and the elemental composition (Fig. 3). In the lower part of the core (19.5–9.6 ka), values are relatively low for TOC and TN, at ∼0.84 % and ∼0.10 %, respectively. At 9.6 ka, there is an appreciable change in the elemental composition of the core, with increases in $\mathrm{Ti}/\mathrm{Ca}$, K, and V and a decrease in $\mathrm{Mn}/\mathrm{Al}$, which corresponds with a transition to darker sediments and an increase in TOC and TN to ∼2.41 % and ∼0.26 %, respectively. At 7.2 ka, there is another major change in the colour and bulk elemental composition of the core, with an increase in redox-sensitive elements U, V, and Mo and a decrease in $\mathrm{Ti}/\mathrm{Ca}$ and K (Fig. 3), which corresponds with darker sediments and increasing TOC values. TOC peaks between 6.6–4.6 ka (∼21 % for TOC and ∼1.7 % for TN), declining towards the top of the core. δ¹⁵N_bulk shows a general decline in values from the upper to the lower part of the core. This decline is small between 19.5–7.7 ka (4.9 ‰–3.3 ‰), before a more significant decrease to 1.2 ‰ at 6.6 ka (3.3 ‰–1.2 ‰). Values increase to 3.7 ‰ at 6.1 ka before declining to 0.0 ‰ at 3.9 ka, increasing slightly towards the top of the core to values of 1.3 ‰.

Figure 4 Geochemical records from Black Sea core 64PE418 of (a) TOC (%), (b) TN (%), (c) δ¹⁵N_bulk (‰), (d) hexose HGs (peak area per g TOC), (e) pentose HGs (peak area per g TOC), (f) BHT-x ratio, (g) crenarchaeol (peak area per g TOC), and (h) isorenieratene (peak area per g TOC).

4.2 Biomarkers

We examined a number of lipid biomarkers related to the N cycle in Black Sea core 64PE418 (Fig. 4). HGs were identified in all samples (with the exception of 215 cm (16.4 ka)). These include HGs with a hexose (C₆) headgroup, i.e. hexose C₂₆ diol, hexose C₂₈ diol, hexose C₂₈ triol, and hexose C₃₀ triol, which are specific to free-living cyanobacteria found in predominately freshwater and brackish environments (Bauersachs et al., 2009; Wörmer et al., 2012). In addition, we detected those with a pentose (C₅) headgroup, i.e. pentose C₃₀ diol, pentose C₃₀ triol, and pentose C₃₂ triol, which are specific to cyanobacteria symbiotic with marine diatoms (Schouten et al., 2013; Bale et al., 2015). Hexose HGs are present throughout the core, increasing substantially in abundance between 9.6–6.6 ka, reaching maximum values at 9.6 ka. Pentose HGs are detected from 4.3 ka onwards, increasing in abundance at the top of the record and coinciding with a low abundance of hexose HGs. Crenarchaeol, a marker for Thaumarchaeota, was identified throughout our record, showing high values in the early part of the record ($\sim \mathrm{1.1}\times {\mathrm{10}}^{\mathrm{14}}$ peak area per g TOC) until 6.9 ka, abruptly shifting to lower values of $\sim \mathrm{3.9}\times {\mathrm{10}}^{\mathrm{13}}$ peak area per g TOC thereafter. The BHT-x ratio, a biomarker for anammox bacteria, is low in the early part of our record (<0.3) due to a low abundance of BHT-x. The BHT-x ratio increases after 6.9 ka to values around 0.3 due to a higher abundance of BHT-x and a lower abundance of BHT with a 34S stereoconfiguration.

Finally, to reconstruct levels of oxygen in the subsurface waters of the Black Sea, isorenieratene was identified (as described in Bale et al., 2021). Isorenieratene is a marker of the brown strains of the photosynthetic green sulfur bacteria, Chlorobiaceae, which are anoxygenic photoautotrophs that require light and hydrogen sulfide (H₂S); their presence indicates photic zone euxinia, whereby anoxic sulfidic waters reach the photic zone (Sinninghe Damste et al., 1993; Koopmans et al., 1996). Isorenieratene was identified in many of our samples after 9.5 ka, peaking between 5.6–4.3 ka (reaching 3.39×10¹² per g TOC at 5.6 ka), but was not detected between 3.9–2.7 ka.

5 Discussion

Based on clear changes in TOC (Fig. 4), colour, and elemental signatures (Figs. 2 and 3), we divided core 64PE418 into three widely acknowledged units, in line with previous studies (Jones and Gagnon, 1994; Arthur and Dean, 1998; Bahr et al., 2005). Unit III spans ∼20–7.2 ka, covering the period where the Black Sea was a lacustrine environment, disconnected from the global ocean, and also the transition interval, where the basin moved towards a marine environment after the IMI over the Bosporus sill at ∼9.6 ka (Aksu et al., 2002; Major et al., 2006; Bahr et al., 2008; Ankindinova et al., 2019). Unit II (∼7.2–2.6 ka) and Unit I (∼2.6 ka–present) span the period where the Black Sea had become an anoxic brackish-to-marine environment. Diagenesis and preservation biases can affect lipid biomarker records, but, since the records do not all mirror the oxygenation of the water column and since organic carbon contents remained relatively high, it is unlikely that diagenesis and preservation significantly influenced these records.

5.1 Oxic lacustrine phase (19.5–9.6 ka)

Throughout the Last Deglaciation and early Holocene (19.5–9.6 ka), TOC and TN levels are low, likely due to poor preservation of organic material, caused by the well-ventilated, oxygenated freshwater environment that existed in the basin at this time (Schrader, 1979). Isorenieratene is not detected during this period, while elements that accumulate in sediment under anoxic conditions (i.e. Algeo and Li, 2020) also remained low (i.e. U, V, Mo; Fig. 2), and this points to a well-oxygenated environment. Freshwater/brackish conditions prevailed throughout this time, as shown by previous studies (Fig. S3 in the Supplement; Filipova-Marinova et al., 2013; Ion et al., 2022; Huang et al., 2021). Throughout this period, the abundance of Thaumarchaeota, indicated by crenarchaeol abundance, and anammox, indicated by the BHT-x ratio, remained relatively steady (Fig. 4). This stability is remarkable, since the region experienced significant climatic changes which led to large variations in the surface water temperatures of the Black Sea, varying from ∼10 °C during the Bølling–Allerød, ∼7 °C during the Younger Dryas and ∼14 °C by the early Holocene (Ménot and Bard, 2012), along with changes in the input of freshwater into the basin due to regional precipitation variability and the melting of Eurasian ice sheets and alpine glaciers (Bahr et al., 2005, 2006, 2008; Badertscher et al., 2011; Shumilovskikh et al., 2012; Filipova-Marinova et al., 2013; Ion et al., 2022). The subsurface stability is likely due to the stratification of the basin, where significant climatic shifts primarily impacted the surface waters, and the limited vertical mixing minimized the influence on subsurface waters, creating a stable environment for the subsurface nitrogen cycle. In contrast, changes in HG abundance and distribution suggest that surface-dwelling nitrogen-fixing cyanobacteria were sensitive to surface water hydrological changes in the Black Sea over this period (Fig. 5). The dominant HG structure varies between hexose C₂₆ diol, hexose C₂₈ diol, and hexose C₃₀ triol, and, after 11 ka, hexose C₂₈ triol becomes present, which has been shown to be the major HG in members of the Rivulariaceae family (i.e. Calothrix sp.) (Bauersachs et al., 2009). The warmer, wetter conditions of the early Holocene may have provided a trigger for this change in HG abundance and composition. Indeed, an increase in the abundance of the genus Rivularia was also noted in coastal regions of SW India during this period, coinciding with an increasingly warm and wet climate (Limaye et al., 2017). Another cause for this shift may have been related to changes in nutrient availability, with members of the Rivulariaceae family typically occurring in environments with highly variable phosphorus availability (Whitton and Mateo, 2012).

Figure 5 Changes over time in relative abundance of hexose and pentose HGs present in core 64PE418 in the Black Sea.

5.2 Transition phase (9.6–7.2 ka)

In line with existing research (Arthur and Dean, 1998; Bahr et al., 2006, 2008), the IMI occurred at ∼9.6 ka, leading to a significant change in colour (Fig. 2) and elemental composition of the sedimentary record (Fig. 3), along with a substantial increase in the abundance of HGs (Fig. 4). This increase does not coincide with higher TOC content, suggesting that enhanced preservation of HGs was not the cause. It is possible that these lipid biomarkers were transported fluvially to this site from lakes within the catchment basin of the Black Sea due to the warm/wet conditions at this time (Göktürk et al., 2011; Shumilovskikh et al., 2012; Filipova-Marinova et al., 2013). This, however, appears unlikely, as our site is located a substantial distance from the mouths of major rivers (>230 km) and the BIT index remains low during this period (∼0.08; Bingjie Yang, personal communications, 2024), indicating only a minor contribution of terrestrial organic matter at our site (Hopmans et al., 2004). Furthermore, as the subsequent period (7–5.6 ka) was also warm and wet (Göktürk et al., 2011; Shumilovskikh et al., 2012; Filipova-Marinova et al., 2013), we would expect the continuation of this peak if the HGs were being sourced from surrounding lacustrine environments. Instead, these high values decline abruptly after 6.6 ka.

It is therefore likely that the peak abundance in nitrogen-fixing cyanobacteria is related to warmer Black Sea surface temperatures during the early to mid-Holocene (Bahr et al., 2008) in combination with surface water stratification (Bahr et al., 2006). This stratification will have slowed the upward supply of fixed nitrogen, reducing nutrient availability in the surface waters, thereby promoting nitrogen fixation by cyanobacteria (Hutchins and Fu, 2017). This stratification may have been driven in part by enhanced freshwater influx due to wetter conditions but may also have been triggered by the IMI through the Bosporus Strait at ∼9.6 ka (Major et al., 2006; Bahr et al., 2008; Ankindinova et al., 2019). This IMI likely led to the gradual salinization of the water column over this transition interval and intermittent build-up of anoxia in the water column. This, in turn, led to periods of higher preservation of organic matter compared to the preceding period, as indicated by the slight increase in TOC after 9.6 ka. The presence of isorenieratene after 9.4 ka indicates that anoxia reached the photic zone at intermittent periods during this transition interval, thereby providing sufficient conditions for the presence of the anoxygenic photoautotrophs, Chlorobiaceae. While the peak in nitrogen-fixing cyanobacteria occurs ∼2 ka before anoxia intermittently entered the photic zone, the initial influx of dense saline water may have led to some reduction in vertical circulation, which reduced the amount of fixed nitrogen upwelled to the upper water column, leading to the presence of nitrogen-fixing cyanobacteria at 9.6 ka. This also coincides with a change in the distribution of HGs in our record between 9.7–6.9 ka, where hexose C₂₈ diol and hexose C₃₀ triol increase in abundance and hexose C₂₈ triol declines in relative abundance and is no longer present after 9.1 ka, coinciding with the presence of isorenieratene. These changes may reflect a shift in species composition, linked to the gradual salinization and periodic anoxification of the water column after the IMI. The IMI at ∼9.6 ka appears, however, to have had little impact on the abundances of anammox and Thaumarchaeota. This is possibly because basin-wide water column stratification and the permanent build-up of anoxia did not occur until later in the record, meaning that neither process instantaneously reacted to the IMI at ∼9.6 ka.

5.3 Shift to anoxic brackish-to-marine mode of operation: a critical N-cycle threshold (∼7.2 ka to present)

After 7.2 ka, there was a substantial increase in TOC and TN and an abrupt shift in parts of the subsurface N cycle. The latter is shown by an increase in the BHT-x ratio, indicating an intensification of anammox, which is coeval with a decrease in crenarchaeol, indicating that there was a decline in Thaumarchaeota-driven nitrification. Studies have shown that, by ∼7.2 ka, anoxia had built up in the water column, as indicated by changes in redox elements (Fig. S2; Eckert et al., 2013; Wegwerth et al., 2018), and water column salinity had significantly increased (Fig. S3; Hiscott et al., 2007; Marret et al., 2009; Soulet et al., 2011; Filipova-Marinova et al., 2013), following the IMI from the Sea of Marmara at ∼9.6 ka (Major et al., 2002, 2006; Bahr et al., 2005, 2008; Ankindinova et al., 2019). This is supported by the presence of isorenieratene in our record during this time, which indicates that anoxia penetrated the photic zone. This water column anoxia likely led to the enhanced preservation of TOC and TN and triggered a shift in the subsurface N cycle, which crossed a threshold from an oxygenated lacustrine mode of operation to an anoxic brackish-to-marine mode of operation. The anoxic water column enabled anammox bacteria to expand their habitat from the anoxic sediments, where they were likely confined when the basin was an oxygenated freshwater environment, up into the suboxic/anoxic water column. This may therefore have commenced part of the modern-day N cycle in the Black Sea, where anammox activity occurs in the lower suboxic zone (∼100 m b.s.l.), where O₂ is (nearly) depleted and H₂S is absent (Jensen et al., 2008), with anammox bacteria consuming ammonium diffusing from the deep sea and utilizing the nitrite produced by both Thaumarchaeota and ammonia-oxidizing bacteria (AOB) (Kuypers et al., 2003; Lam et al., 2007). Consequently, it may be that the abundance of anammox bacteria increased as a result of the coupling to nitrite production by other microbes in the suboxic zone, whilst benefitting from ammonium diffusing upwards from the deep sea. The increased anammox after 7.1 ka likely indicates that more bioavailable nitrogen was lost from the system after the switch to the anoxic brackish-to-marine mode of operation. At the same time, Thaumarchaeota abundance declined, which may be in part due to the build-up of anoxia in the water column that reduced the niche of these aerobic microbes and the nitrification performed by them. Once these processes crossed a threshold from an oxygenated lacustrine mode of operation to an anoxic brackish-to-marine mode of operation, they appear to have remained steady for the remainder of the Holocene despite changes in the salinity of the basin (van der Meer et al., 2008; Mertens et al., 2012; Coolen et al., 2013) and significant changes in regional temperature and precipitation (Göktürk et al., 2011; Shumilovskikh et al., 2012; Filipova-Marinova et al., 2013). This shows that deoxygenation was the main driver of the observed change in anammox and in archaeal nitrification and that they were not affected by hydrological changes mainly occurring at the surface. While our record shows centennial-scale changes in the N cycle, we acknowledge that there may have also been rapid or short-term variations in N-cycle dynamics over this period that may not have been captured by the resolution of this record.

At 6.1 ka, the abundance of the HGs substantially declined, coinciding with an increase in δ¹⁵N_bulk, indicating a reduction in nitrogen fixation. This corresponds with the high-resolution record of Fulton et al. (2012), which also shows a decline in δ¹⁵N_bulk values at this time. As this decline in HG abundance and increase in δ¹⁵N_bulk does not coincide with a reduction in TOC, it is unlikely that reduced preservation of HGs played a role here. As nitrogen-fixing cyanobacteria inhabit the upper surface layer, it is likely that this change is linked to the salinization of the surface waters, with many studies demonstrating the disappearance of many freshwater mollusc, ostracod, and dinoflagellate cyst species at this time, which were replaced by an increased abundance of euryhaline Mediterranean species (Hiscott et al., 2007; Marret et al., 2009; Filipova-Marinova et al., 2013; Ivanova et al., 2015). At 6.1 ka, hexose C₂₆ diol and hexose C₂₈ diol are the only HGs present in the record, which may reflect the dominance of genera in the Nostocaceae family (i.e. Anabaena sp., Aphanizomenon sp., Nodularia sp., Nostoc sp.), as these members demonstrate a dominance of the hexose C₂₆ diol and also contain varying amounts of hexose C₂₈ diol (Gambacorta et al., 1999; Bauersachs et al., 2009). This distribution is similar to that of the Baltic Sea after ∼7.2 ka, when a series of weak intrusions of saline water led to the basin becoming fully brackish (Sollai et al., 2017). It is therefore possible that the peak in HGs in our Black Sea record between 9.6–6.9 ka represents a transition from the dominance of freshwater-tolerant nitrogen-fixing cyanobacteria to more brackish species, with brackish species dominating the surface waters after 6.6 ka. After 6.1 ka, δ¹⁵N_bulk gradually decreases, indicating a rise in nitrogen fixation, with this pattern also shown in previous records (Blumenberg et al., 2009; Fulton et al., 2012). While previous studies have shown riverine nitrogen to be a major source of fixed nitrogen in the modern Black Sea (McCarthy et al., 2007), due to the remoteness of our core site from the coast, our records are unlikely to have been significantly influenced by riverine input. It should be noted that a previous study (Fulton et al., 2012) suggested, based on compound specific measurements of pyropheophytin, that sedimentary δ¹⁵N in the Black Sea is primarily derived from eukaryotic algae rather than cyanobacteria, which exhibits a different fractionation of nitrogen isotopes. This means the use of δ¹⁵N_bulk as a nitrogen-fixation signal must be used with caution. HGs, however, are only derived from N-fixing cyanobacteria and are therefore an unambiguous biomarker of nitrogen fixation. Interestingly, at 4.3 ka, pentose HGs are detected, coinciding with lowest δ¹⁵N_bulk, indicating the presence of marine nitrogen-fixing cyanobacteria found in symbiosis with marine diatoms. This indicates that the surface water salinity had reached a threshold which enabled these marine microbes to survive, with research indicating salinity reached ∼17 ‰ during the deposition of Unit I (Ion et al., 2022), and freshwater/brackish species had disappeared by this time (Fig. S3; Filipova-Marinova et al., 2013). Indeed, reported increases in the number of euryhaline species at this time also point to the increasing salinity of the surface waters (Marret et al., 2009; Bradley et al., 2012), which may be linked to warmer/drier conditions that reduced freshwater influx and/or enhanced evaporation (Göktürk et al., 2011). Between 3.9–2.7 ka, isorenieratene is not detected in the samples, reflecting the findings of previous studies (Sinninghe Damsté et al., 1993). It has been suggested that this resulted from the erosion of the chemocline (Sinninghe Damsté et al., 1993), while other research shows a short reoccurrence of freshwater/brackish species (Fig. S3; Filipova-Marinova et al., 2013), which may indicate that enhanced freshwater input was responsible for lowering the chemocline below the photic zone. The disappearance of hexose HGs after 0.6 ka indicates that surface water salinities may more recently have become too high for the proliferation of brackish nitrogen-fixing cyanobacteria.

6 Conclusions

This study shows a relatively stable subsurface N cycle in the Black Sea over the Last Deglaciation and Holocene, with the exception of a critical threshold observed at 7.2 ka when the basin shifted from an oxygenated lacustrine environment to an anoxic brackish-to-marine basin. At this time, the loss of bioavailable nitrogen through anammox activity was enhanced and Thaumarchaeota-driven nitrification was reduced. Prior to and after this transition, the subsurface N cycle was remarkably stable despite various climatic and hydrological changes that impacted the basin during the deglaciation and the Holocene. Both the amount of nitrogen fixation by cyanobacteria and the composition of these microbes in the surface waters, however, appear to be much more dynamic and sensitive to hydrological changes over this period, responding in particular to salinity and temperature changes and stratification of the water column. Consequently, these records provide important insight into how future deoxygenation in marine environments may affect the microorganisms involved in the N cycle, possibly leading to enhanced loss of bioavailable nitrogen by anammox and reduced nitrification by Thaumarchaeota. Furthermore, in areas where enhanced water column stratification limits the supply of fixed nitrogen in the surface waters, localized enhanced cyanobacterial nitrogen fixation in the surface waters may occur. These changes may have associated feedbacks on nutrient cycling and carbon fixation, with implications for the future global carbon budget.