Diatom responses and geochemical feedbacks to environmental changes at Lake Rauchuagytgyn (Far East Russian Arctic)

This study is based on multiproxy data gained from a ¹⁴C-dated 6.5 m long sediment core and a ²¹⁰Pb-dated 23 cm short core retrieved from Lake Rauchuagytgyn in Chukotka, Arctic Russia. Our main objectives are to reconstruct the environmental history and ecological development of the lake during the last 29 kyr and to investigate the main drivers behind bioproduction shifts. The methods comprise age-modeling, accumulation rate estimation, and light microscope diatom species analysis of 74 samples, as well as organic carbon, nitrogen, and mercury analysis. Diatoms have appeared in the lake since 21.8 ka cal BP and are dominated by planktonic Lindavia ocellata and L. cyclopuncta. Around the Pleistocene–Holocene boundary, other taxa including planktonic Aulacoseira, benthic fragilarioid (Staurosira), and achnanthoid species increase in their abundance. There is strong correlation between variations of diatom valve accumulation rates (DARs; mean 176.1×10⁹ valves m² a¹), organic carbon accumulation rates (OCARs; mean 4.6 g m⁻² a⁻¹), and mercury accumulation rates (HgARs; mean 63.4 µg m⁻² a⁻¹). We discuss the environmental forcings behind shifts in diatom species and find moderate responses of key taxa to the cold glacial period, postglacial warming, the Younger Dryas, and the Holocene Thermal Maximum. The short-core data likely suggest recent change of the diatom community at the beginning of the 20th century related to human-induced warming but only little evidence of atmospheric deposition of contaminants. Significant correlation between DAR and OCAR in the Holocene interglacial indicates within-lake bioproduction represents bulk organic carbon deposited in the lake sediment. During both glacial and interglacial episodes HgAR is mainly bound to organic matter in the lake associated with biochemical substrate conditions. There were only ambiguous signs of increased HgAR during the industrialization period. We conclude that if increased short-term emissions are neglected, pristine Arctic lake systems can potentially serve as long-term CO₂ and Hg sinks during warm climate episodes driven by insolation-enhanced within-lake primary productivity. Maintaining intact natural lake ecosystems should therefore be of interest to future environmental policy.

1 Introduction

Today, northern and mountain regions warm faster than elsewhere on Earth, putting cold freshwater systems at risk for loss of ecosystem services (IPCC, 2021). Paleoenvironmental research, however, still lacks sufficient geographical coverage in the eastern Russian Arctic (Kaufman et al., 2020; Mckay et al., 2018; Sundqvist et al., 2014). Arctic lakes are powerful archives of past climate information because they respond rapidly to external forcing on their catchments (Biskaborn et al., 2021b; Nazarova et al., 2021; Subetto et al., 2017). Effects of both past climate changes and modern human impacts during the industrial period, including mercury contamination of pristine ecosystems, have been demonstrated for remote Siberian lake ecosystems (Biskaborn et al., 2021a). In paleolimnological research many studies are based on concentrations of fossil remains and geochemical compounds in the sediments, yet within the eastern Arctic, only a few studies have managed to accomplish the reconstruction of accumulation rates of these sediment constituents, possibly owed to limited age controls (Vyse et al., 2021).

The Pleistocene–Holocene (P–H) transition from glacial to interglacial climates commonly reveals the major change within biotic and geochemical sediment components and is well detectable in sufficiently old lake sediments. Shorter and less powerful climate events are less distinctly represented in low-accumulation systems, and their impacts on lake ecosystems in sparsely covered areas are not yet sufficiently understood (Kaufman et al., 2004; Subetto et al., 2017; Biskaborn et al., 2016; Renssen et al., 2012). One of the most known groups of photosynthetic organisms in Arctic lakes are diatoms (Smol and Stoermer, 2010). They are siliceous microalgae (Bacillariophyceae) that form opaline valves which preserve excellently in lake mud and allow identification up to the highest species levels by light microscope analysis (Battarbee et al., 2001). Diatoms as a group are one of the major primary producers in aquatic environments, contributing to the global net primary production at about 25 % (Smol and Stoermer, 2010). Diatom communities respond to numerous environmental forcings including hydrochemical changes, seasonal climate shifts (duration of ice cover), and physical habitat changes (Hoff et al., 2015; Douglas and Smol, 2010; Pestryakova et al., 2018; Herzschuh et al., 2013; Biskaborn et al., 2012, 2013; Palagushkina et al., 2017). Diatom productivity has been estimated from Si $/$ Al ratios (Vyse et al., 2020), valve concentrations (Biskaborn et al., 2012), and biogenic opal concentrations (Meyer et al., 2022). However, there is yet only sparse information available in the literature addressing the contribution of aquatic bioproduction, i.e., diatom primary producers, to accumulation rates of organic matter over different climate stages (Biskaborn et al., 2021b).

There is an ongoing discussion about the role of Arctic lakes in the carbon cycle. Thermokarst basins are believed to have switched from a net source to a sink during the mid-Holocene ca. 5000 years ago related to permafrost dynamics (Anthony et al., 2014). Glacial lakes are often larger and well oxygenized and thus are considered to strongly contribute to the modern CO₂ emission in the Arctic landscape (Tan et al., 2017; Wik et al., 2016). Differences in drivers of bioproductivity, e.g., land use in Europe (Vihma et al., 2016), accumulation rate, and preservation of sedimentary carbon, e.g., during the ice melt (Spangenberg et al., 2021), still lead to a high sink–source variability across temporal scales. To help gain insights into the fate of carbon accumulated in northern lakes, we provide a high-resolution study of a sediment core from Lake Rauchuagytgyn.

Paleoenvironmental records that include the Last Glacial Maximum (LGM) are sparse in the vast region of Chukotka (Vyse et al., 2020). Lake archives exceeding the LGM were published for Lake El'gygytgyn (Melles et al., 2007), Lake Ilirney (Vyse et al., 2020; Andreev et al., 2021), and Lake Rauchuagytgyn (Vyse et al., 2021), of which the latter is the subject of our study. There is still insufficient knowledge about feedbacks of specific climate events, such as the Younger Dryas (YD) cooling (Andreev et al., 2021; Kokorowski et al., 2008) or the Holocene Thermal Maximum (HTM) (Renssen et al., 2012), to lake primary producers. However, short-term and fast events that could compare to the pace of recent climate change are of specific interest to understand today's climate–ecosystem relationships.

Over the last few years influences of past and recent climate changes to diatom assemblage shifts have been investigated in lake records in Yakutia (Kostrova et al., 2021; Courtin et al., 2021; Biskaborn et al., 2021b), accompanied by lake ecosystem feedbacks and long-distance heavy metal contamination (Biskaborn et al., 2021a). This study was accordingly set up to test whether similar paleolimnological responses to climate and anthropogenic impacts exist in very remote areas in Chukotka. In our paper we present new diatom records and biogeochemical data based on the published chronological sediment records and climate reconstructions of Lake Rauchuagytgyn (Andreev et al., 2021; Vyse et al., 2021). Our objectives on Lake Rauchuagytgyn are to (1) reconstruct accumulation of diatoms since the last glacial in comparison to bulk organic carbon accumulation, (2) investigate the main drivers behind assemblage and bioproduction shifts, and (3) compare past natural and recent mercury loads to test for potential heavy metal contamination of remote pristine ecosystems.

2 Study site

The catchment of Lake Rauchuagytgyn (67.82^∘ N, 168.7^∘ E; elevation 625 m a.s.l.; surface area 6.1 km²; maximum water depth 36 m; catchment area 214.5 km²) is located in the northwestern Anadyr Mountains of Chukotka in the northeastern Russian Arctic (Fig. 1). The lake's main inflows are situated at the southern margin, and a few outflows drain the lake to the north and the sides. Glacial activity in the catchment is preserved by moraine structures north of the lake and in surrounding glacial cirques (Glushkova, 2011; Vyse et al., 2021). The basement of the study site consists of silicic–intermediate lithology, represented by Cretaceous andesite (Zhuravlev et al., 1999). The area is characterized by strong arctic continental climate with mean annual air temperatures of −11.8 ^∘C, and mean July and January temperatures are 13 and −30 ^∘C, respectively, while annual precipitation is at ca. 200 mm (Menne et al., 2012). Open herb and graminoid tundra, with tree occurrence only in lower elevations and close to rivers, characterizes the surrounding landscape (Huang et al., 2020; Shevtsova et al., 2020).

Figure 1 Study site. (a) Bathymetrical map of Lake Rauchuagytgyn with catchment area (boundary as black line, inflows as blue lines) and coring locations (long core EN18218, short core 16-KP-04-L19B). (b) Geographical overview map. Map based on ESRI (ESRI and GeoEye, 2019).

Hydrochemical data from July 2018 (Supplement S2) showed that the lake water had diluted freshwater with low conductivity (85.5 µs cm⁻¹), medium transparency (Secchi depth 3.9 m), slightly alkaline conditions (pH 7.8), and low dissolved organic carbon (0.9 mg L⁻¹).

3 Materials and methods

3.1 Fieldwork

Fieldwork and coring activities at Lake Rauchuagytgyn (Fig. 1) were performed by helicopter expeditions in July 2016 and July 2018. We used a handheld echo sounder and a UWITEC gravity corer (60 mm) to retrieve a short core 16-KP-04-L19B with 23 cm length in summer 2016 at a 31.0 m deep part of the lake at (67.7888 N, 168.7380 E). After a few hours the core was subsampled in 0.5 to 1 cm slices before being transported in dark and cool conditions. A longer parallel core from the same site and expedition was analyzed for pollen and radiocarbon-based chronology (Andreev et al., 2021).

In summer 2018 we used an Innomar SES-2000 compact parametric sub-bottom profiler to locate the coring location in the southern sub-basin and retrieved a long core (EN18218, ca. 6.5 m) at 29.5 m water depth using a UWITEC Niederreiter 60 mm piston coring system operated on a platform at anchor (67.7894 N, 168.7335 E). Coring, processing, and sediment geochemistry of this sediment core were already described by Vyse et al. (2021).

3.2 Chronology

For the chronology of the long core EN18218, we used LANDO (linked age and depth modeling). In the current version (v1.3), LANDO combines the output of five age–depth modeling software programs (Bacon, Bchron, clam, hamstr, Undatable) in a single interactive computing platform described in Pfalz et al. (2022). The advantage of this approach over a single age–depth model is that the combined model takes multiple age–depth uncertainty ranges into consideration and reduces biases towards overinterpretation. We have updated the published sedimentation rate (SR) values for EN18218 accordingly, based on the same 23 radiocarbon dates in Vyse et al. (2021) shown in Table 1. According to Vyse et al. (2021) we added the same age offset to the data derived from the surface sample (785 ± 31 years BP), which corresponds to 853 ± 31 years when the 2018 common era (CE) expedition year is taken into account.

Table 1 Radiocarbon dates from the sediment core EN18218 from Vyse et al. (2021) used to generate age–depth relationships and sedimentation rates in LANDO.

To date the CE over the industrial period in the short core 16-KP-04-L19B, freeze-dried subsamples were analyzed for ²¹⁰Pb and ¹³⁷Cs activities by direct gamma assay in the Liverpool University Environmental Radioactivity Laboratory, using ORTEC GWL series high-purity germanium (HPGe) well-type coaxial low-background intrinsic germanium detectors (Appleby et al., 1986). ²¹⁰Pb and ¹³⁷Cs were measured from their gamma emissions at 46.5 and at 662 keV, respectively. Accuracies of used detectors were determined using calibrated standard sources of known activity. The effect of self-absorption of low-energy gamma rays was used for corrections (Appleby et al., 1992). To model the chronology of the short core, we used the R package rPlum (Blaauw et al., 2021) to apply a Bayesian framework to determine the chronology based on ²¹⁰Pb measurements (Hunter et al., 2022; Aquino-López et al., 2018). We used 18 measurements of supported and unsupported ²¹⁰Pb (Table 2) within Plum, while assuming a varying supply of unsupported ²¹⁰Pb for the model (Fig. 2a, b). We corrected the Plum model by constraining it with the ¹³⁷Cs peak between 3 and 3.5 cm (Fig. 2c), which is attributed to the high point in atomic weapon testing in 1963 (Appleby, 2001; Hunter et al., 2022).

Figure 2 Fallout radionuclides in the short core 16-KP-04-L19B showing (a) total and supported ²¹⁰Pb, (b) unsupported ²¹⁰Pb, and (c) ¹³⁷Cs concentrations versus depth.

3.3 Biogeochemistry and mercury analysis

To gain information about the productivity in the lake we analyzed total organic carbon (TOC), total carbon (TC), and total nitrogen (TN) from 20 short-core samples 16-KP-04-L19B and from 66 samples from the long-core samples EN18218. A total of 25 samples below 220 cm in EN18218, however, revealed TN values below the detection limit (0.1 wt %) and are therefore not displayed. TOC and total inorganic carbon (TIC) were detected using a Vario soli TOC cube elemental analyzer (Elementar Analysensysteme GmbH) following combustion at 400 ^∘C for organic carbon and 900 ^∘C for TIC. The sum of TOC and TIC was used to estimate TC. TN was measured using a rapid MAX N exceed (Elementar Analysensysteme GmbH). The data were used to calculate the TOC $/$ TN ratio, using a factor of 1.167, which is the ratio of the atomic weights of nitrogen (14.007 amu, atomic mass unit) and carbon (12.001 amu), to obtain the atomic ratio TOC $/$ TN_atomic following Meyers and Teranes (2002). TOC from EN18218 was used from Vyse et al. (2021) to estimate TOC $/$ TN_atomic ratios for the long core.

Total mercury (THg) was analyzed in 20 samples from core 16-KP-04-L19B and 32 samples from core EN18218. We determined the THg in solid material by thermal decomposition, amalgamation, and atomic absorption spectrophotometry using a direct mercury analyzer (DMA-80 evo; MLS-MWS GmbH). The solid samples were weighted into 1 mL metal boats, which are then combusted at about 750 ^∘C under a flow of oxygen, and the Hg in the off-gases is trapped as an amalgam on a gold sieve. In a subsequent step, Hg is released, and its amount is determined by atomic absorption spectroscopy. We used the certified reference material BCR^®–142R (67 µg kg⁻¹ Hg) as reference material after every 18th measurement and four standards every beginning of a measuring day. The detection limit of the most sensitive cuvette was < 0.003 ng Hg. For each sample, we measured THg at least two times and up to four times if the results showed larger variations.

3.4 Diatom analysis

We analyzed diatoms in a total of 54 samples from the sediment core EN18218 and 20 samples from the surface core 16-KP-04-L19B, taken from 0.5 cm slices. For light-microscopy-based species identification we prepared diatom slides following the procedure described in Battarbee et al. (2001). We treated 0.1 g of freeze-dried sample material with hydrogen peroxide (30 %) for up to 5 h, added hydrochloric acid (10 %) to stop the reaction, and washed the sample with purified water. Finally, microspheres between 5×10⁶–8×10⁶ were added, according to the density of valves on the test slides, to estimate the concentration of diatom valves (DVC). Homogenized sediment suspension was transferred to cover slips placed in Battarbee cups to avoid species fractionation and mounted to slides using Naphrax. To identify diatoms to the lowest possible taxonomic level we used a ZEISS Axioscope 5 light microscope with an Axiocam 208 color camera attached, equipped with a Plan-Apochromat $\mathrm{100}\times /\mathrm{1.4}$ Oil Ph3 objective at 1000× magnification. We counted more than 300 diatom valves in each sample (Wolfe, 1997) in both sediment cores (mean 351 valves in EN18218; mean 375 valves in 16-KP-04-L19B). Diatom species identification was based on various literature including Hofmann et al. (2011) and Krammer and Lange-Bertalot (1986–1991) as well as online databases (i.e., http://www.algaebase.org, last access: 1 August 2022). Correct identification of species was supported by images from a scanning electron microscope and for some species with input from the diatom community online platform DIATOM-L (Bahls, 2015). During diatom analysis, valves were distinguished in pristine and non-pristine valves, and chrysophyte cysts were counted but not further identified.

Table 2 Fallout radionuclide concentrations in the short core 16-KP-04-L19B.

3.5 Data processing and statistics

For statistical analysis of downcore proxy data we used the R environment (R Core Team, 2016). Both cores, EN18218 and 16-KP-04-L19B, were statistically analyzed following the same procedure.

To create diatom zones along the cores we used the package rioja for constrained incremental sums-of-squares clustering (CONISS) based on Euclidean dissimilarity after log transformation of species percentage data to downweigh abundant species (Grimm, 1987). Attribution of diatom zones along the core depth was guided by CONISS results, while the total number of zones in the cores is referring to meaningful chronologies in the region (Andreev et al., 2021; Anderson and Lozhkin, 2015; Andreev et al., 2012).

We used the decorana function in the package vegan (Oksanen et al., 2020) to apply a detrended correspondence analysis (DCA) on percentage data and calculated gradient length in standard deviation units (SD EN18218 = 2.36; SD 16-KP-04-L19B = 1.25). According to the threshold suggested by Birks (2010), we chose principal component analysis (PCA) to reveal major trends in the data. In the PCA we also chose Euclidean distance but square-root transformation of the data to downweigh abundant species less aggressively than log transformation performed in CONISS. Before PCA, we filtered the species data to exclude rare species; i.e., only species present with ≥ 3 % in ≥ 2 samples were included in PCA. The bottom sample at 540 cm was too different from the rest of the more established diatom assemblage and was therefore excluded from the analysis.

We estimated diatom species richness (alpha diversity) based on Hill's N0 and N2 diversity (Hill, 1973) and performed rarefaction to correct richness estimates for differences of valve counts (Birks et al., 2016) using the vegan R package (Oksanen et al., 2020). The minimum base sum of all samples was n=304 for EN18218 and n=333 for 16-KP-04-L19B.

To estimate diatom valve dissolution, we calculated the F index following Ryves et al. (2001), providing a range between 0 and 1 in which 0 is poor and 1 is perfect preservation.

Pearson correlation matrices were generated using the cor function in the R package stats. To highlight significant correlations by p-value criterion, a significance test based on the upper-tail probability from the Pearson correlation coefficients was performed using the function cor.mtest in the R package corrplot. Only correlations that yielded p values < 0.05 were considered significant.

The TOC $/$ N_atomic ratios were calculated following Meyers and Teranes (2002) using weight ratios of TOC and N multiplied by 1.167. TOC from EN18218 was used from Vyse et al. (2021) to estimate TOC $/$ N ratios for the long core.

To calculate accumulation rates, we first computed dry mass accumulation rates (MARs; in g cm⁻² a⁻¹) using Eq. (1):

$$\begin{array}{}\text{(1)}& \mathrm{MAR}=\mathrm{DBD}\times \mathrm{SR},\end{array}$$

where DBD is dry bulk density (in g cm⁻³), and SR is sedimentation rate (in cm a⁻¹). We derived SR from age–depth modeling in a standard procedure according to Eq. (2) (Pfalz et al., 2022).

$$\begin{array}{}\text{(2)}& \mathrm{SR}\left(x_i\right)={\displaystyle \frac{\mathrm{depth}\left(x_i\right)-\mathrm{depth}\left(x_{i-\mathrm{1}}\right)}{\mathrm{age}\left(x_i\right)-\mathrm{age}\left(x_{i-\mathrm{1}}\right)}}\end{array}$$

The value x_i represents the layer of interest within a sediment core for which the SR calculation is necessary, while x_i−1 is the previous layers. Since the DBD measurement of the 16-KP-04-L19B surface sample was missing, we extrapolated the value from samples below by constructing a piecewise polynomial in the Bernstein basis using the Python package scipy (Virtanen et al., 2020). To determine uncertainty ranges for the individual accumulation rates we propagated the 2σ uncertainty range of SR into the MAR calculations.

Organic carbon accumulation rates (OCARs) were estimated by dividing TOC (wt %) by a factor of 100, then multiplying it by MAR, and converting it to g m⁻² a⁻¹ units by multiplying it by a factor of 10 000. We multiplied Hg values by MAR to estimate mercury accumulation rates (HgARs) but then converted HgARs to µg m⁻² a⁻¹ units by multiplying the HgAR values by a factor of 10. Diatom accumulation rates (DARs in valves m⁻² a⁻¹) were estimated following Birks (2010) by using MAR multiplied by the diatom valve concentration (valves g⁻¹) and concerting it to 10⁹ valves m⁻² a⁻¹ units.

The mean summer insolation was calculated for Rauchuagytgyn at 90^∘ from the vernal point (67.8^∘ N, 1365 solar constant) using QAnalySeries 1.5.1 and Earth's orbital parameters from Laskar et al. (2004). Pollen-reconstructed July temperatures (TJuly_pollen) and annual precipitation (APP_pollen) from Andreev et al. (2021) were resampled onto the core depths based on their published ages and our age-model output from EN18218. Mean July temperatures from the closest weather station OSTROVNOE (first observation 1936 CE; 68.12^∘ N, 164.17^∘ E; ID: RSM00025138; 98 m a.s.l.; 195 km W of Lake Rauchuagytgyn) were estimated from daily values retrieved from NOAA (https://www.noaa.gov, last access: 1 July 2022) and resampled onto the sample depths of 16-KP-04-L19B based on the plum age model.

4 Results

4.1 Chronology

The LANDO age–depth model based on radiocarbon dates for EN18218 (Fig. 3) shows an age range of 28 190 to 29 907 cal yr BP (weighted mean age: 28 950 cal yr BP) at 651.75 cm, which agrees with the age–depth model developed by Vyse et al. (2021) of about 29 000 cal yr BP at the core base. The weighted mean sedimentation rates of all LANDO models range over the entire core from 0.01 to 0.1 cm a⁻¹, which means that the maximum value is higher in some regions than previously reported (Vyse et al., 2021), with a maximum of 0.054 cm a⁻¹. However, both models agree on decreases in mean sedimentation rates below 0.02 cm a⁻¹ around 558–560 and 346–358 cm and increases of around 510–519 and 371–374 cm, where LANDO models suggest mean values above 0.065 cm a⁻¹. Additional decline in the mean sedimentation rate is found below 0.02 cm a⁻¹, approximately between 224–243 cm. Taking into account the 2σ confidence intervals of all models, the uncertainty for the sediment core ranges from minimum values of 0.002 cm a⁻¹ at 569 cm to maximum values of 0.575 cm a⁻¹ between 52–54 cm.

Figure 3 Generated output from LANDO for the sediment core EN18218 based on ¹⁴C data from Vyse et al. (2021). Left plot consists of a comparison between five age–depth models from different modeling codes indicated in the legend. Colored solid lines indicate the median age, while shaded areas represent their respective 1σ and 2σ ranges in the same colors with decreasing opacities. Panel (b) shows the calculated sedimentation rate with matching colors. Black circles in (a) indicate the mean calibrated ages of ¹⁴C bulk sediment samples based on an IntCal20 calibration curve (Reimer et al., 2020) and their 1σ uncertainty error bars.

The Plum age–depth model based on lead and cesium dates for the short core 16-KP-04-L19B (Fig. 4) reaches a maximum mean age of 1864 CE (uncertainty range: 1813–1905 CE) at 11 cm. Total ²¹⁰Pb activity reached equilibrium with the supporting ²²⁶Ra at a depth of around 7 cm (Fig. 2), which explains the larger uncertainty between 7 and 11 cm (Fig. 4). Unsupported ²¹⁰Pb concentrations vary irregularly with depth with significant non-monotonic features between 1–2.5 and 3–4.5 cm. Mean sedimentation rates range between 0.03 cm a⁻¹ (2–3 cm) and 0.139 cm a⁻¹ (0–1 cm), where the higher values can be explained by the lack of ²¹⁰Pb data for the first centimeter. The uncertainty range (2σ confidence interval) for the sedimentation rate over the entire short core lies between 0.025 and 0.192 cm a⁻¹.

Figure 4 Plum age–depth model for the sediment core 16-KP-04-L19B. The five upper panels show the Bayesian input parameters and their posterior distributions for Plum. Panel (b) consists of the age–depth model with its mean age in red and its 2σ confidence interval in grey, the unsupported ²¹⁰Pb concentrations (in Bq kg⁻¹) in blue with its 1σ uncertainty, and the supported ²¹⁰Pb concentrations (in Bq kg⁻¹) in violet. Panel (c) displays the mean sedimentation rate over depth as a dashed line and the 2σ confidence interval in grey.

4.2 Diatom species assemblages

Diatoms occurred upward of 541 cm (21.8 ka cal BP) in the long core EN18218 (Fig. 5b) and were found in all samples between 0 and 10.5 cm in the short core 16-KP-04-L19B (Fig. 5a). In total 204 different species were identified. The dominant taxa in the observed samples are represented by planktonic cyclotelloid species, Aulacoseira, and small achnanthoid species. Chrysophyte cysts only occurred with a few counts in two samples (251 and 311 cm in EN18218) and were therefore neglected. The valve dissolution F index in both the long (min 0.78, max 0.95, mean 0.90) and short core (min 0.86, max 0.98, mean 0.94) was generally high referring to an overall good valve preservation.

Figure 5 Relative abundance of diatom species. (a) Species assemblages in the short core 16-KP-04-L19B. Species percentage values are shown next to the calibrated mean ages (common era years, CE) and the core depth below the sediment surface. (b) Relative abundance of diatom species in the long core EN18218. Species percentage values are shown next to the mean calibrated ages before present and the core depth below the sediment surface. Diatom zones are established by CONISS clustering. Taxa present with ≥ 3 % in ≥ 2 samples were included in the graphs.

Mean diatom valve concentrations (DVCs) in EN18218 were 46.4 (2.3–134.7) 10⁷ valves g⁻¹, corresponding to diatom valve accumulation rates (DARs) of 176.1 (2.0–651.6) 10⁹ valves m⁻² a⁻¹. Modern DVCs found in the short core had a relatively higher mean of 82.0 (43.8–200.4) 10⁷ g⁻¹, corresponding to DARs of 199.2 (42.2–514.9) 10⁹ m⁻² a⁻¹.

The rarefied species richness, Hill's N0, varied between 11.8 and 42.2 (mean 29.4) in the long core and was slightly higher between 31.1 and 47.4 (mean 38.8) in the short core. The effective richness, Hill's N2, ranged between 1.5 and 8.4 (mean 3.1) in the long core and between 2.9 and 6.1 (mean 4.3) in the short core. A remarkable shift toward high effective richness is found in the long core between ca. 241 and 346 cm.

The first three directions in the principal component analysis explain ca. half of the data variance in EN18218; i.e., PC1, PC2, and PC3 explained 23.2 %, 17.8 %, and 11.8 %, respectively. The first three PCA axes from the short core assemblage data explained ca. two-thirds of the data variance, i.e., 28.9 %, 22.8 %, and 14.4  % for PC1, PC2, and PC3, respectively. PC3, however, was not included in the biplots shown in Fig. 6.

According to our cluster analysis we divided the cores into six diatom zones in core EN18218 and two zones in the surface core 16-KP-04-L19B and described the species in chronological order.

Figure 6 Biplots of the first two dimensions (PC1, PC2) generated by principal component analysis of diatom species filtered to ≥ 3 % in ≥ 2 samples from the long core EN18218 and the short core 16-KP-04-L19B. Color circles represent eco-taxonomical clusters with comparable environmental preferences. Colored sample depths indicate chronologies. The explained variance of each PC is indicated at the axis label in percentage.

Diatom zone 1: 541–426 cm (EN18218, 21.8–15.3 ka cal BP)

The oldest diatoms found in the long core (Fig. 5b) were dominated by planktonic Lindavia ocellata with ca. 76 % in the bottom sample (mean ∼ 42.9 in the zone), Lindavia cyclopuncta (∼ 40.8 %), Lindavia bodanica (∼ 2.4 %), and Aulacoseira valida appearing only in two samples in the middle of the zone (6.9 % and 12.8 %). Pliocaenicus costatus occurred more frequently from 441 cm onward (7.9 %). Benthic species started to appear with Achnanthidium minutissimum, Encyonema minutum, and Nitzschia palea with low abundance, followed by Staurosira construens, Encyonopsis descriptiformis, Psammothidium chlidanos, and Hannaea arcus.

Diatom zone 2: 426–366 cm (EN18218, 15.3–12.8 ka cal BP)

Dominant planktonic species were represented by L. ocellata (∼ 44.1 %), L. cyclopuncta (∼ 32.4 %), and P. costatus (∼ 4.9 %). Staurosira pinnata and Staurosira brevistriata occurred in addition to the benthic species found in zone 1.

Diatom zone 3: 366–346 cm (EN18218, 12.8–11.4 ka cal BP)

Zone 3 appeared as a small section (two samples) in which only L. ocellata remained frequent (∼ 61.6 %), while other planktonic forms were restricted to first-time occurrence of Aulacoseira subarctica (∼ 4.5 %) and less frequently P. costatus. Brachysira neoexilis (∼ 4.1 %) started to occur, while S. pinnata started to become more frequent (∼ 3.0 %).

Diatom zone 4: 346–241 cm (EN18218, 11.4–8.0 ka cal BP)

In the Early Holocene part of the core dominant planktonic species were L. cyclopuncta (∼ 37.7 %) and L. ocellata (∼ 14.8 %), accompanied by generally more frequent benthic forms represented by S. pinnata (∼ 7.9 %), B. neoexilis (∼ 6.8 %), and A. minutissimum (∼ 3.8 %).

Diatom zone 5: 241–166 cm (EN18218, 8.0–5.3 ka cal BP)

Zone 5 started with the dominance of L. cyclopuncta (∼ 57.7 %), while other Lindavia species disappeared. Instead, A. subarctica (∼ 4.8 %), P. costatus (∼ 3.2 %), and A. valida (∼ 2.9 %) occurred in higher frequencies. S. pinnata reached the highest values (∼ 8.4 %), while P. chlidanos and H. arcus also showed peaking values.

Diatom zone 6: 166–11 cm (EN18218, 5.3–1.1 ka cal BP)

The upper part of the long core EN18218 was characterized by the re-occurrence of L. ocellata (∼ 15.0 %) and the strong representation of A. subarctica (∼ 9.2 %), while L. cyclopuncta (∼ 53.8 %) remained the dominant species. Benthic Staurosira formed and A. minutissimum decreased.

Diatom zones 7–8 in the short core: 10.5–3.0–0.5 cm (16-KP-04-L19B, 1870–1970–2012 CE)

The strongest shift within the species assemblage of the short core is found at 3 cm (Fig. 5a). Overall in the short core, L. cyclopuncta (∼ 37.3 %) and A. subarctica (∼ 24.4 %) remain the most dominant species. L. ocellata appears with the highest abundance in zone 7, i.e., between 7.5 and 4.5 cm (∼ 17.5 %). L. cyclopuncta decreases at the beginning of zone 8, while P. chlidanos reaches the highest values up to 4.8 %.

4.3 Biogeochemical variables

We complemented geochemical data from core EN18218 (Fig. 7b) provided by Vyse et al. (2021) for TN and THg measurements. TN varied from > 0.1 to 0.25 wt %, with the highest values in the upper 100 cm of the core. Resulting TOC $/$ TN_atomic ratios ranged between 6.0 and 19.2, with a strong increase at 341 cm. THg in the same core ranged between 93.2 and 362.8 µg kg⁻¹, with the highest value in the sample at 600.25 cm and overall higher mean values above 321 cm (mean 198.6 compared to 141.8 µg kg⁻¹ below). Mean Hg accumulation rates (HgARs) estimated from these concentrations were 63.4 (11.8–138.6) µg m⁻² a⁻¹. Mean organic carbon accumulation rates (OCARs) estimated from TOC values published by Vyse et al. (2021) were 4.6 (0.8–12.7) g m⁻² a⁻¹. Low peaks in OCAR, DAR, and HgAR in both cores correspond to low sedimentation rates at 230, 350, and 550 cm in EN18218 and between 2 and 3 cm in 16-KP-04-L19B (Figs. 3 and 4).

Figure 7 Biogeochemical variables and statistical diatom indices since the Late Pleistocene from Lake Rauchuagytgyn. (a) Surface sediment core 16-KP-04-L19B covering the last ca. 150 years. (b) Long sediment core EN18218 covering the last 28 000 years. OCAR, organic carbon accumulation rate; DAR, diatom accumulation rate; HgAR, mercury accumulation rate; F index, diatom valve preservation index; TOC $/$ TN_atomic, total organic carbon to total nitrogen ratio; diatom species richness based on Hill numbers; PC1–2, main axis sample scores from the principal component analysis; light cyclotelloid Lindavia and euplanktonic Aulacoseira as sum percentages; planktonic to benthic species ratios. July temperatures (mean, min, max) in panel (a) are calculated from weather observation at the OSTROVNOE station (https://www.noaa.gov, last access: 1 July 2022); dotted red line indicates T increase and onset of the Anthropocene; pollen-reconstructed July temperatures T_July and annual precipitation P_ann in panel (b) were adopted from Andreev et al. (2021); OCAR and TOC $/$ TN_atomic were based on organic carbon concentrations from Vyse et al. (2021); and insolation was calculated from orbital parameters at the lake's latitude 67.8 N following Laskar et al. (2004).

Data from the short core 16-KP-04-L19B are presented for the first time (Figs. 5a and 7a). TOC ranged between 2.6 wt % and 3.5 wt %, with the highest values in the upper 3 cm. N varied between 0.28 wt % and 0.41 wt %, resulting in TOC $/$ TN_atomic ratios with only little fluctuation around the mean 10.7, which fits into the upper part of EN18218. THg in the short core varied between 162.4 and 244.7 µg kg⁻¹, with the highest values in between 4.5 and 3 cm. Mean OCAR and HgAR estimated from TOC and THg concentrations were 6.7 (2.7–11.5) and 46.0 (14.3–69.8) µg m⁻² a⁻¹, respectively.

5 Discussion

5.1 Ecological responses of diatom species to Late Quaternary environmental changes

Diatoms in Lake Rauchuagytgyn started to appear at 21.8 ka cal BP (Fig. 5b) with strong dominance of L. ocellata (Pestryakova et al., 2018), a planktonic and ultraoligotrophic to mesotrophic taxon, common in cold lakes (Wunsam et al., 1995). The first occurrence of diatoms was accompanied with the first increase of organic carbon accumulation (Fig. 7b). According to previous sedimentological work on the sediment core (Vyse et al., 2021), at this time glaciers retreated from the catchment, and unfrozen episodes became more frequent, leading to paraglacial deposition progressing in the lake basin. In the course of continued deglaciation since ca. 20 ka cal BP in Chukotka (Vyse et al., 2020) and Alaska (Elias and Brigham-Grette, 2013), the diatom assemblage developed progressively in diatom zone 1, enabling oligotrophic L. cyclopuncta (Scussolini et al., 2011) and a few benthic species to occupy ecological niches in the young and still cold lake ecosystem. Strong fluctuations, e.g., of tychoplanktonic A. valida, indicated unstable habitat conditions during that period. Pliocaenicus costatus is known in larger quantities from cold and strongly oligotrophic mountain lakes restricted to eastern Siberia (Cremer and Van De Vijver, 2006). Low abundance of benthic diatoms may result from thick ice due to long ice cover periods and reduced light penetration, as well as in-wash of clay during deglaciation (Vyse et al., 2021), leading to mildly transparent and narrow littoral zones in an overall deep basin. In diatom zone 2 the species richness increased strongly, and benthic diatoms became abundant (Hill's N0, planktonic $/$ benthic ratio in Fig. 7b), supporting a gradual climate amelioration equivalent to the Bølling–Allerød interstadial, which started ca. 15.5–15.0 ka cal BP (Wohlfarth et al., 2007; Obase and Abe-Ouchi, 2019; Andreev et al., 2021), facilitating shallow water habitats and thus more complex diatom communities due to longer growing seasons (Cherapanova et al., 2007).

Over the deglaciation period, in parallel to the development of catchment vegetation, the lake ontogeny was likely driven by changes in the load of dissolved organic carbon (DOC). As shown in lake evolution studies (Engstrom et al., 2000), young lakes in freshly deglaciated terrain have low DOC and rather alkaline conditions, which is reflected by the benthic species assemblage in the record, such as fragilarioid species successively accompanied by Encyonopsis descriptiformis and Brachysira neoexilis. Modern DOC measured in July 2018 (0.9 mg L⁻¹) clearly below the global lake average of 3.9 mg L⁻¹ (Toming et al., 2020) together with other hydrochemical parameters (Supplement Table S2) indicates an overall diluted and alkaline lake system, suggesting more depleted conditions in the past.

The short but remarkable diatom zone 3 is characterized by the same cold-adapted planktonic and parts of benthic species from the early deglacial period in diatom zone 1. Thus, in accordance with other findings from Chukotka (Anderson and Lozhkin, 2015) and, e.g., Lake El'gygytgyn (Andreev et al., 2012), the Rauchuagytgyn diatom assemblage provides evidence of an aquatic ecosystem response to climate cooling and drying between ca. 12.8–11.4 ka cal BP. Corresponding to the Younger Dryas (YD) period, our diatom data show disappearance of L. bodanica and L. cyclopuncta but relative increase of heavy Aulacoseira valves (Figs. 5b and 7b), indicating turbulent water conditions. Complex diatom responses within the YD associated with an increase of Aulacoseira species have been found in Lake Baikal (Mackay et al., 2022). In many boreal lakes YD cooling weakened lake thermal stratification, leading to turbulent conditions, resulting in similar diatom responses as observed in Lake Rauchuagytgyn (Neil and Lacourse, 2019).

The Pleistocene–Holocene (P–H) boundary is detected from the diatom assemblage change at ca. 346 cm in core EN18218, fitting well into the uncertainty range of 10.8–12.2 ka cal BP (Figs. 5b and 3). At the glacial–interglacial transition, the diatom community responded with a strong decrease of planktonic and light Lindavia species (Biskaborn et al., 2021b) that was accompanied with a decrease in both diatom and carbon accumulation rates (Fig. 7b). Mountain ice sheets that persisted in the catchment over the deglacial period vanished at the P–H boundary, leading to decreased water supply and lower lake levels over the Early Holocene. At that time, the effective species richness (Hill's N2) increased because relatively more benthic species reached higher percentages, while the pure richness (N0) slightly decreased.

The P–H is also characterized by a distinct increase of the first axis sample scores of the PCA (Fig. 7b), pointing to the most prominent increase of benthic diatom taxa in the record. The PCA biplot depicts grouping of planktonic Lindavia versus benthic Staurosira and Psammothidium species along the primary axis, while Aulacoseira species are oriented along the secondary axis (Fig. 6). This general shift to benthic communities can be explained by temperature-driven changes in the duration of the ice cover period. Longer open-water seasons in the Early Holocene promote light penetration and the availability of littoral habitats, while input of DOC and nutrients enhances benthic production in the littoral zone (Hu et al., 2018; Engstrom et al., 2000).

Fragilarioid taxa such as S. pinnata, S. construens, and S. brevistriata are known as typically small benthic pioneering forms in boreal shallow lakes (Valiranta et al., 2011; Biskaborn et al., 2012) that are often alkalophilous (Paull et al., 2017). Together with the increase of achnanthoid taxa, this assemblage indicates increased availability of littoral habitats. Increased chemical weathering of bare rocks during warmer and wetter interglacial conditions and the development of roots (Andreev et al., 2021) in fresh soils both led to enhanced ion supply (Herzschuh et al., 2013) and eventually increased alkalinity of the lake water.

At the P–H transition, abrupt high TOC $/$ TN_atomic values of around 15 (Fig. 7b) point to a higher contribution of less-degraded organic carbon, indicating an increase in benthic water plants and terrestrial plant material (Meyers and Teranes, 2002; Baird and Middleton, 2004), and thus provide evidence for shallower shores and the development of catchment vegetation due to maximal summer insolation and warm interglacial conditions (Fig. 7b). The increased role of water plants is supported by epiphytic B. neoexilis, E. minutum, A. minutissimum, and E. descriptiformis (Barinova et al., 2011; Hofmann et al., 2011).

Swann et al. (2010) reconstructed most favorable climate conditions known as the Holocene Thermal Maximum (HTM) at Lake El'gygytgyn (140 km to E) between 11.4 and 7.6 ka cal BP. However, based on the pollen data, Andreev et al. (2021) reconstructed the start of warmest conditions at ca. 8.0 ka cal BP for the Rauchuagytgyn region. At 8.0 ka cal BP in diatom zone 5, L. ocellata disappeared together with the low levels of L. bodanica, while Aulacoseira species started to establish themselves. Aulacoseira builds heavy and rapidly sinking frustules commonly found in deep boreal lakes (Laing and Smol, 2003). Euplanktonic A. subarctica is a pelagic species that requires turbulence to remain in the photic zone (Rühland et al., 2008; Gibson et al., 2003), while light cyclotelloid taxa prefer stratified water conditions (Rühland et al., 2015). We assume that high July temperatures continued but in addition the open-water seasons prolonged around 8–7 ka cal BP, as winters in Siberia gradually became warmer over the mid-Holocene and Late Holocene (Meyer et al., 2015). Early ice-out and the influx of meltwater during spring and summer associated with increased and longer spring circulation supported Aulacoseira species (Horn et al., 2011) and led to a distinct change in the Rauchuagytgyn species assemblage. For comparison, in recent times, a surface ice layer in Chukotka lakes builds up in October, reaching up to 1.8 m over the winter, and breaks up in early July after snowmelt started in mid-May (Nolan et al., 2002).

At ca. 6.4 ka cal BP the DAR, OCAR, and TOC $/$ TN_atomic ratios increased (Fig. 7b), while L. bodanica reappeared, and A. valida increased strongly, but small benthic Staurosira species retreated (Fig. 5b). We associate this change with the maturation of soils (Biskaborn et al., 2012), retreating woody vegetation (Andreev et al., 2021), and a shift in bioproduction driven by an increased supply of nutrients through increased river activity and increasing water levels (Buczkó et al., 2013). During the cooling of the Late Holocene the diatom zone-6 assemblage continued as a semi-pelagic cold-water community at intermediate-to-high water levels, high DAR but slightly decreased richness, and also decreasing terrestrial influence indicated by decreasing TOC $/$ TN_atomic values.

The last few decades are represented in the short core 16-KP-04-L19B about 220 m E of the long-core position. The surface sediments in this area of the lake were slightly different but also dominated by the same Lindavia and Aulacoseira species as compared to diatom zone 6 (Fig. 5a). In 1907 CE benthic taxa Psammothidium chlidanos and Pinnularia nodosa increased, accompanied by a slight shift from Aulacoseira to Lindavia species and decreasing OCAR, DAR, and HgAR but slightly increased TOC $/$ TN_atomic values (Fig. 7a). The Aulacoseira–Lindavia shift is tentatively supported in the PCA biplot in PC1 (Fig. 6). Even though the overall response to recent environmental changes in Rauchuagytgyn seems to be of minor extent, the timing and response correspond to warming at high latitudes observed during industrialization (Biskaborn et al., 2021a). Abrupt shifts in lake ecosystems were most frequently observed around 1950 CE (Huang et al., 2022) when the beginnings of human energy consumption and geochemical impacts initiated the (proposed) Anthropocene epoch (Syvitski et al., 2020). A strong warming in 1950 CE was also documented in air temperature observations in the weather station 195 km W of the study area. In Rauchuagytgyn, DAR decreased strongly at that time, while HgAR and TOC $/$ TN_atomic ratios fluctuated, having a negative influence on species richness N0 and N2 (Fig. 7a). At the boundary between diatom zones 7 and 8 there is a peak in Tabellaria flocculosa (1960–1985 CE), a species that can occur with both planktonic and benthic lifestyles (Heudre et al., 2021), indicating slightly acidic and nutrient-enriched environmental conditions, and may respond to unstable habitat conditions (Palagushkina et al., 2012). As a pennate planktonic diatom, Tabellaria often responds with increased abundance to atmospheric nitrogen deposition (Rühland et al., 2015), corresponding to increased nitrogen levels between 1970 and 1980 CE (Fig. 7a). After 1970 CE in diatom zone 8 Lindavia decreased but P. chlidanos and A. subarctica increased, possibly related to changed nutrient and mixing conditions (Rühland et al., 2015). Since the beginning of the 21st century DAR, OCAR, and Hg have increased again, while P. nodosa has decreased, pointing to either minor atmospheric influences on the lake hydrochemistry or natural short-term variation (Gibson et al., 2003).

5.2 Correlation between carbon, diatoms, and mercury accumulation

Accumulation rates (ARs) in sediment basins are generally uncertain due to limitations in precise age-model interpolations (Sadler, 1981). In addition, diatom concentrations are expressed as numbers of frustules (Battarbee et al., 2001) regardless of the weight and volume of the shells. Accordingly, one cannot directly infer biomass from count-based valve accumulation, as valves vary considerably in size among and within species (Birks, 2010). We showed above that the Rauchuagytgyn sedimentary record shows a tendency toward successional lake development in response to long-term changes in the landscape and ecosystem adaptation (Brenner and Escobar, 2009). Therefore, unknown deviations in the linkage between the mass of carbon stored and the number of diatom valves observed are likely to appear.

At the millennial timescale in the long core, OCAR is strongly correlated with HgAR (r 0.94, p<0.05) and significantly with DAR (r 0.64, p<0.005), while there is no significant correlation between diversity indices to HgAR (Supplement Fig. S1). The mean values at around 11.2 (and 12.9 in the Holocene) of TOC $/$ TN_atomic measured in the long core represent a mixture of (in simple words) high-N planktonic algae, medium-N benthic water plants, and low-N terrestrial vegetation input (Baird and Middleton, 2004). Phytoplankton produces TOC $/$ TN_atomic ratios of 4–10, whereas vascular land plants produce > 20 (Meyers and Teranes, 2002). Given the sparse vegetation cover around the lake (Huang et al., 2020) and the overall low TOC $/$ TN_atomic ratios, terrestrial input may play a minor role. We therefore assume that there is a strong contribution of algae to the bulk organic matter accumulated in the lake, which tends to be somewhat proportionate to the number of diatom valves. Accordingly, the TOC $/$ TN_atomic ratios correlate negatively with planktonic $/$ benthic ratios in the long core (r −0.57, p<0.05) but only slightly (insignificantly) in the short core (r 0.21, p>0.05). This mismatch could indicate that there is an anthropogenic nitrogen contribution from the atmosphere (Biskaborn et al., 2021a) that is in addition masked by short-term fluctuations and constraints of measurement precision in high-resolution samples from lakes under extreme environmental conditions. A tentative relationship between DAR and planktonic species could be detected in the short core (r 0.41, p<0.05), which could indicate that the widespread increase of planktonic species in high-latitude lakes as a response to global warming (Smol et al., 2005) contributed to the increase in diatom primary productivity.

Over the last few centuries visible in the short core there has also been a significant correlation between OCAR and HgAR. Atmospheric mercury, however, is not simply deposited in Arctic lakes, but instead there is a strong influence of limnological processes such as primary production and ice cover dynamics on mercury biogeochemical cycling (Korosi et al., 2018). The correlation between OCAR and DAR (r 0.65, p<0.05) apparently shows that diatoms play a role in these processes. In contrast to long timescales there is a significant negative correlation in the short core between HgAR and diversity estimates such as Hill's N0 (r −0.51, p<0.05) and N2 (r −0.39, p<0.05). Thus, contaminants during the industrial period could be assumed to have a stronger effect on the lake ecosystem than natural Hg supply before increased anthropogenic activity (Huang et al., 2022). Studies on deep permafrost soils in Siberia showed that average Hg concentrations of 9.7 µg kg⁻¹ could be used as a baseline for natural Hg concentrations (Rutkowski et al., 2021). However, as Hg binds to lake organic carbon (Braaten et al., 2018), lake bioproductivity is likely increasing the mercury load within sediments, explaining the overall high concentrations in older sections. Furthermore, we found a very good correlation between HgAR and OCAR during the cold glacial period (r 0.98, p<0.05) but obvious decoupling from diatoms as shown by a missing correlation between HgAR and DAR (Fig. 8; Supplement Fig. S1). Mercury in tundra catchments is closely related to non-vascular plants (Olson et al., 2019), and an external supply of plant organic matter was reported to represent the main source of cold-climate carbon deposition (Hughes-Allen et al., 2021). In Rauchuagytgyn, however, the higher amount of nitrogen detected in the pre-Holocene core section suggests one or both of the following two reasons: (1) within-lake aquatic production by algae other than well-preserved diatoms flourished during the glacial (Hernández-Almeida et al., 2015), and/or (2) the preservation of nitrogen was higher during the prolonged ice cover period (Kincaid et al., 2022) than during the interglacial.

Figure 8 Schematic drawing of the long-term processes leading to accumulation of diatom valves (DAR), organic carbon (OCAR), and mercury (HgAR) in Lake Rauchuagytgyn. Pearson correlation coefficients are indicated by r in black when correlation was significant (p<0.05) and in red italic when p>0.05 (insignificant).

5.3 Long-term ecosystem feedbacks to climate changes

Well-preserved and old diatom records in Chukotka provide the opportunity to study direct responses of natural lake ecosystems to regional climate changes (Cherapanova et al., 2007; Swann et al., 2010). The Lake Rauchuagytgyn sediment record provides insight into compositional changes of diatom assemblages in response to lake and catchment changes. The main changes observed are best represented by shifts within planktonic species and their proportions relative to benthic forms populating emerging habitats. Significant negative correlations were found at millennial timescales between planktonic $/$ benthic ratios and diversity estimates such as Hill's N0 (r −0.7, p<0.05) and N2 (r −0.53, p<0.05), indicating that long-term diatom diversity in Lake Rauchuagytgyn was closely related to lake ontogeny. The general catchment maturation accompanied by decreases of glacial ice-sheet influence led to the initiation of new ecological niches and thus diversification of, e.g., epiphytic species (Wilson et al., 2012; Rouillard et al., 2012).

Diatom accumulation rates on the other hand, recorded independently from lifestyles, show a clear relationship to lake bioproduction. Relationships between diatoms and organic carbon in lake sediments were also related to species alpha diversity in Lake Bolshoe Toko (Biskaborn et al., 2021b) and explained as stabilizing effects of well-developed species richness supporting primary biomass production (Marzetz et al., 2017). The correlation found between diatom valve and organic carbon accumulation rates during the interglacial in the Rauchuagytgyn sediment record may support that, during warm episodes, diatoms in high-latitude lakes with relatively small catchments are coupled with the bulk production of biomass (Fig. 8). Our study revealed a positive feedback mechanism between long-term climate amelioration and diatom-driven sink of organic matter. Compared to ocean systems, where fertilization projects attempted to force carbon burial by artificial diatom blooms (Yoon et al., 2018), lakes may possess a higher potential to withdraw carbon from the atmosphere because of lower carbon remineralization rates (Mendonça et al., 2017; Sobek et al., 2009). However, this may nowadays be questioned because whole-lake experiments and models have suggested a possible lagged response of lakes' natural resistance to anthropogenic stressors that could cause fast ecosystem switches (Pahl-Wostl, 2003). In turn, these potential alterations could possibly prevent carbon sink feedbacks, as observed in the remote and still pristine Rauchuagytgyn system. In this context, the observed positive relationship between sedimentary carbon and mercury suggests potential mitigation feedbacks of contamination stress accompanied by recent climate change. However, impacts of human-driven atmospheric stressors only seem a little pronounced in the short-core data, which is limiting possibilities to assign natural long-term mechanisms to present-day conditions. This is amplified by the fact that boreal lakes have either already passed important ecosystem thresholds or are about to exceed ecological tipping points upon further warming (Wischnewski et al., 2011) and are believed to not represent pristine ecosystems anymore (Smol et al., 2005).

6 Conclusions

Radiocarbon- and ²¹⁰Pb-dated sediment cores from Lake Rauchuagytgyn in the Far East Russian Arctic provide valuable archives of millennial- to decadal-scale lake ecosystem responses to regional environmental forcing of the last 22 000 years before today. Our main findings based on diatom species, organic carbon and nitrogen, and mercury analyses can be highlighted as follows:

• The Pleistocene diatom species assemblage reflects a planktonic community in a deep and cold lake with short growing seasons. The assemblage becomes more complex during a gradual climate amelioration at ca. 15 ka cal BP, similarly to Bølling–Allerød, leading to the successive development of benthic habitats. Diatom species temporarily returned to glacial conditions between ca. 12.8–11.4 ka cal BP, corresponding to the Younger Dryas.

• The Early Holocene diatom community reflects a shallower lake with larger littoral zones and higher alkalinity that we relate to prolonged ice-free periods and vegetation development in the catchment, supported by high carbon to nitrogen ratios. Gradual increasing Aulacoseira taxa indicate that winters became warmer over the mid-Holocene and Late Holocene, leading to earlier ice-out and longer spring circulation.

• During Late Holocene cooling, small benthic Staurosira taxa retreated due to soil maturation and increased water levels, facilitating a higher abundance of planktonic Lindavia and Aulacoseira species.

• The last few decades represented in a ²¹⁰Pb-dated short core only show vague evidence of recent change in the diatom community in 1907 CE, indicated by a slight increase of light Lindavia and decrease of Aulacoseira species, accompanied by shifts in the benthic community. Biogeochemical variables and diatom indices fluctuated strongly around 1950 CE.

• Diatom accumulation rates (DARs) and organic carbon accumulation rates (OCARs) do not correlate during the cold episode but show significant correlation during the warm interglacial when insolation was higher. The Rauchuagytgyn data suggest that during the Holocene (1) deposition of organic carbon was largely driven by within-lake bioproduction, and (2) diatoms reflect the activity of the gross primary producers of the lake.

• Mercury accumulation rates (HgARs) in the investigated sediments are strongly correlated to OCARs in both cold and warm episodes. As Hg accumulation is bound to organic matter, increased carbon sedimentation during warm climates and suitable biochemical substrate conditions facilitate Hg deposition.

• From our study we infer that bulk carbon accumulation is represented by climate-enhanced within-lake primary productivity. Thus, pristine boreal lake systems potentially can serve as long-term CO₂ sinks if short-term fluctuations are disregarded. Lake basins also represent disposal sites for heavy metal contaminants. Consequently, maintaining intact natural lake ecosystems should be a high priority in future environmental policy.