Diatom shifts and limnological changes in a Siberian boreal lake: a multiproxy perspective on climate warming and anthropogenic air pollution

Lake ecosystems are affected globally by climate warming and anthropogenic influences. However, impacts on boreal lake ecosystems in eastern Siberia remain underexplored. Our aim is to determine if shifts in diatom assemblages in a remote lake in eastern Siberia are related to climate warming, similar to observations in temperate regions, while also exploring how the ecosystem might be influenced by hydroclimate and human-induced air pollution. We analysed continuous sediment samples from a ²¹⁰Pb–¹³⁷Cs-dated short core from Lake Khamra (59.99° N, 112.98° E), covering ∼220 years (ca. 1790–2015 CE), following a multiproxy approach on the same sample material to provide a comprehensive record of environmental changes. Biogeochemical proxies include total organic carbon (TOC) and total nitrogen (TN) concentrations and corresponding stable isotopes of bulk sediment samples (δ¹³C, δ¹⁵N), as well as diatom silicon isotopes (δ³⁰Si_diatom), alongside light microscope diatom species analysis. The diatom assemblage at Lake Khamra is dominated by few planktonic species, primarily Aulacoseira subarctica and Aulacoseira ambigua. At ca. 1970 CE, we observe a major shift in diatom assemblages, characterised by a marked increase in the planktonic species Discostella stelligera and a decrease in both Aulacoseira taxa. We attribute these changes to recent global warming, which is likely associated with earlier ice-out and enhanced summer thermal stratification, consistent with similar observations in temperate lake ecosystems. A rapid increase in chrysophyte scales (Mallomonas) from the 1990s onward further supports an increasing thermal stratification of the lake driven by rising temperatures. Biogeochemical proxies indicate substantial limnological changes around 1950 CE, preceding the major shift in diatom communities, likely driven by hydroclimatic variability. Increased precipitation and weathering are further discussed in order to explain changing silica sources leading to decreasing δ³⁰Si_diatom after ∼1970 CE. Nevertheless, the interpretation of δ³⁰Si_diatom in lacustrine systems is complex, likely influenced by both in-lake biogeochemical processes and catchment dynamics. Indications of anthropogenic influences on Lake Khamra include a δ¹³C depletion, likely linked to fossil fuel combustion and emissions, coinciding with industrial growth in Asia and Russia. Nonetheless, we find no evidence for atmospheric nitrogen deposition. We conclude that the Lake Khamra ecosystem is severely affected by climate warming and shows indications of human influence. This emphasises the urgent need for comprehensive research to mitigate these impacts on remote lake ecosystems in order to secure natural water resources.

1 Introduction

Human activities, particularly the burning of fossil fuels, have contributed significantly to the “great acceleration” of Earth system changes especially since the 1950s (Steffen et al., 2015). Since then, the incidence of major ecological shifts in lake ecosystems globally has increased, driven by both climate change and anthropogenic impacts (Huang et al., 2022). These lake observations align with the proposed onset of the “Anthropocene” (Crutzen and Stoermer, 2000), though it has not been formally recognised as a distinct geological epoch (ICS, 2024).

Globally lake surface temperatures have risen in recent decades (O'Reilly et al., 2015; Hampton et al., 2018), having numerous consequences like changing thermal stratification properties, longer ice-free periods and changing lake ecosystems like shifting biological communities (Smol et al., 2005; Saros et al., 2012; Hampton et al., 2017; Woolway et al., 2020). Extensive research on lakes in North America and Europe (Smol et al., 2005; Smol and Douglas, 2007; Rühland et al., 2008; Kahlert et al., 2020) has shown that climate warming and human activities lead to significant ecological changes in remote high-latitude ecosystems. In Siberia, for example, Lake Baikal, a thoroughly studied system, has experienced documented ecosystem changes driven by recent warming, including increased water temperatures and shorter ice-cover periods (Todd and Mackay, 2003; Mackay et al., 2006; Hampton et al., 2008, 2014; Moore et al., 2009; Izmest'eva et al., 2016). However, a notable gap remains in understanding how these global changes impacted boreal lakes in remote and less-studied regions of eastern Siberia, such as the Republic of Sakha (Yakutia).

Yakutia, located in eastern Siberia (Fig. 1a), has experienced rapid climate warming, with annual air temperature trends showing an increase of 0.3–0.6 °C per decade since 1966 (Gorokhov and Fedorov, 2018). Most meteorological stations report the highest temperature increases in winter, along with a rise in overall precipitation rates over the last 50 years (Gorokhov and Fedorov, 2018). Furthermore, climate warming is associated with an increase in wildfire activity (AMAP, 2021), also observed in Yakutia, which is partly due to anthropogenic alterations to the ecosystems including traditional agro-industrial burning practices and the expansion of industrial areas (Kirillina et al., 2020). Evidence suggests that recent warming and human-induced pollution are impacting limnological conditions in Siberia. For instance, a remote lake in eastern Siberia, without direct human influence in the catchment, displays signs of industrialisation's legacy, with mercury contamination in its sediments and a trend of lake acidification (Biskaborn et al., 2021b).

Figure 1 (a) Location of Lake Khamra in eastern Siberia, with Yakutia highlighted in darker grey. (b) Catchment of Lake Khamra. (c) Bathymetry of Lake Khamra, main inflow in the south-west, main outflow in the north-east and coring location of EN18232-1 in the central part of the lake. Map layer credits: Esri, Maxar, Earthstar Geographics and the GIS User Community.

Diatoms, microscopic unicellular algae prevalent in nearly all aquatic environments, are effective indicators of various environmental disturbances such as climate change, acidification and nutrient enrichment (Smol and Stoermer, 2010). Their remains, called frustules, are composed of biogenic silica (SiO₂⋅n H₂O), which are well-preserved in lake sediments. Diatom valves are suitable for light microscopic identification up to the highest species level (Battarbee et al., 2001), as well as for oxygen and silicon isotope measurements, to reconstruct past climate and environmental conditions (Leng and Barker, 2006; Sutton et al., 2018; Frings et al., 2024).

Circum-Arctic studies reveal climate warming, resulting in a shortened duration of lake ice covers and an extended growing season, along with the thermal stratification of lake water, is the primary cause of the observed shift in diatom assemblages (Sorvari et al., 2002; Rühland et al., 2003, 2008, 2015; Rühland and Smol, 2005; Smol et al., 2005). Paleolimnological studies in eastern Siberia, including Yakutia, have utilised diatom assemblages, alongside diatom isotopes, biogeochemical proxies such as organic carbon and nitrogen, and mercury analyses, to investigate past environmental conditions. However, many of these studies predominantly cover periods within the Quaternary, primarily focusing on the Holocene, and investigate long-term trends (e.g. Cherapanova et al., 2006; Biskaborn et al., 2012, 2016, 2021c; Pestryakova et al., 2012; Firsova et al., 2021; Kostrova et al., 2021; Mackay et al., 2022; Biskaborn et al., 2023). While there is a growing body of work on recent environmental change in eastern Siberia, high-resolution studies covering shorter time periods in Yakutia remain rare. Besides a study in the remote boreal areas in Yakutia (Biskaborn et al., 2021b), notable study sites are located further south at Lake Baikal (Roberts et al., 2018), on Kamchatka (Jones et al., 2015), further west towards Europe (Palagushkina et al., 2020) and in the northern Urals (Solovieva et al., 2008). However, these sites are influenced by different climatic and ecological conditions that can differ from those in Yakutia.

In a previous study at the boreal Lake Khamra in south-west Yakutia, a diatom oxygen isotope record (δ¹⁸O_diatom) was established to reconstruct hydroclimatic anomalies and suggested increasing (winter) precipitation along with rising air temperatures from the 1970s onward (Stieg et al., 2024b). In addition to oxygen isotopes, silicon isotopes of diatoms have been used as indicators of changing productivity and nutrient utilisation. Predominantly, δ³⁰Si_diatom serves as a proxy for primary production in marine environments (De la Rocha et al., 1998; Sutton et al., 2018). In lacustrine systems, however, δ³⁰Si_diatom interpretation is more complex due to the combined influence of catchment processes (e.g. weathering, vegetation, soil development) and in-lake biogeochemical dynamics (e.g. silica cycling and sedimentation) (Sutton et al., 2018; van Hardenbroek et al., 2018; Frings et al., 2024). Research at Lake Baikal has explored the preservation and reliability of δ³⁰Si_diatom for paleo-reconstructions, as well as its link to diatom productivity and sedimentary processes (Panizzo et al., 2016, 2017, 2018). Despite these studies, spatial and temporal coverage of δ³⁰Si_diatom research in eastern Siberia remains limited. Furthermore, studies that analyse both oxygen and silicon isotopes of diatoms from the same samples are rare, with one example from northern Siberia focusing on the time period since the Last Glacial Maximum (Swann et al., 2010).

In this study, we analyse subfossil diatom assemblages in Lake Khamra to determine whether taxonomic changes are consistent with recent climate warming, as documented in many temperate lakes in North America and Europe. Additionally, we investigate potential limnological responses to climate and anthropogenic impacts. Our research addresses this by examining diatom assemblages and biogeochemical proxies (total organic carbon (TOC), total nitrogen (TN), δ¹³C, δ¹⁵N, δ³⁰Si_diatom) at decadal resolution from ca. 1790 CE onward. Using continuous sediment samples from a short core of Lake Khamra, eastern Siberia, we provide a detailed account of environmental changes. Our specific objectives are to (I) identify historical lake ecosystem changes within a continuous diatom assemblage record spanning the last ∼220 years, (II) evaluate the origin of organic material and distinguish between allochthonous and autochthonous sources in response to hydroclimatic variations using biogeochemical proxies, and (III) assess the impact of human-induced air pollution on the lake ecosystem.

2 Materials and methods

2.1 Study site

Lake Khamra (59.99° N, 112.98° E; 340 m a.s.l.) is located in the south-west of the Republic of Sakha (Yakutia), eastern Siberia, Russia, with the closest urban settlements, Peleduy and Vitim, situated 40 and 60 km to the south-west, respectively. The landscape is located in the transition zone of evergreen to deciduous needle leaf boreal forest within the ecoregion of the middle taiga (in Geng et al., 2022, after Stone and Schlesinger, 2003). According to field observations (Kruse et al., 2019; Miesner et al., 2022), the research area is covered by a mixed coniferous forest. The annual mean ground temperature lies between 0.0–1.0 °C at a depth of 10–20 m (Shestakova et al., 2021), indicating that the study site is within a discontinuous to sporadic permafrost zone (Obu et al., 2019). Geologically, the region is underlain by Cambrian bedrock composed of alternating dolomite and limestone, interspersed with silty Ordovician sandstone and patches of clayey Silurian limestone (Chelnokova et al., 1988).

The lake morphology is characterised by shallow shore areas and a nearly central deep part with a maximum water depth of 22.3 m (Fig. 1c). The lake has a surface area of 4.6 km² and a catchment size of 107.3 km². Lake Khamra is a hydrologically open system with a main inflow in the south-west and an outflow in the north-east and an estimated average long-term discharge of 1.1 m³ s⁻¹ (Messager et al., 2016). In winter, the lake is covered by ice and snow, observed with an average ice thickness of 0.5 m and an about 1 m thick snow cover in March 2020 (Biskaborn et al., 2021a). Temperatures measured in the water column during that time ranged from 2–2.7 °C. A surface water sample taken during a summer field campaign in 2018 revealed a pH value of 6.07 and a conductivity of 40 µS cm⁻¹ (Stieg et al., 2024b).

The nearest weather station is located in the town of Vitim (59.45° N, 112.58° E; 186 m a.s.l.; ECA (European Climate Assessment) station code 3235). Details and a climate diagram for the study site can be found in Stieg et al. (2024b). In relation to this, the regional climate is continental, expressed in extremely cold and dry winters (January: −28.8 °C mean temp., 1929–2018 CE) and warm and humid summers (July: +18.1 °C mean temp., 1929–2018 CE). The annual precipitation average at Vitim reaches 423 mm (1929–2018 CE), while the mean annual air temperature is at −5.0 °C.

2.2 Fieldwork, subsampling and core dating

Coring activities at Lake Khamra were conducted during a summer expedition in 2018 (Kruse et al., 2019). The sediment short core EN18232-1 was obtained with an UWITEC gravity corer (60 mm) from the central and deepest part of Lake Khamra (59.99091° N, 112.98373° E; 22.3 m), based on water-depth measurements determined with a surveying rope and a portable HONDEX PS-7 LCD digital sounder. The sediment core with a total length of 42 cm was sealed with water-absorbent floral foam and transported in a PVC tube to the Alfred Wegener Institute (AWI) in Potsdam, where it was stored dark and cool at 4 °C until further analysis. The sediment core was subsampled downcore in horizontal 1 cm continuous increments (n=39) in October 2021. To avoid potential contamination between depths due to mixing, rim material (<0.5 cm) was carefully removed from each sample layer. All subsamples were freeze-dried for at least 48 h before further processing.

We use the age–depth model for the short core EN18232-1 published by Stieg et al. (2024b), providing a comprehensive ∼220-year record (ca. 1790–2015 CE, mean ages) with a sub-decadal resolution according to mean ages (5.7±1.7 years). The chronology relies on analysis for ²¹⁰Pb, ²²⁶Ra, ¹³⁷Cs and ²⁴¹Am using direct gamma assay techniques with Ortec HPGe GWL series well-type coaxial low-background intrinsic germanium detectors, conducted at the Liverpool University Environmental Radioactivity Laboratory as described by Appleby et al. (1986). The age–depth model of the short core was calculated with a Bayesian accumulation model within the R package “rbacon” v2.5.8 (Blaauw and Christen, 2011; R version 4.1.1), based on the ²¹⁰Pb chronology shown in Fig. A1 in the appendix. In addition, data on the water content and the average dry bulk density (DBD) of the short core were determined from selected subsamples (n=24) using a 1 cm³ tool taken at regular intervals, with weights measured before and after freeze-drying. The mean DBD (g cm⁻³) was then used together with the sedimentation rates (SRs; cm a⁻¹) to calculate mass accumulation rates (MARs; $\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$) downcore (full details given in Stieg et al., 2024a).

2.3 Diatom analysis

Diatom slides were prepared for each of the 39 samples of the short core EN18232-1 to analyse the species assemblages. Preparation of slides for siliceous microfossils followed Battarbee et al. (2001). Aliquots of 0.10 g freeze-dried sample material were treated with H₂O₂ (30 %) at 90 °C for up to 8 h to remove organic matter. The reaction was stopped and carbonates were eliminated by adding HCl (10 %). After washing with purified water, 5 mL of microsphere suspension ($\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{microspheres}{\mathrm{mL}}^{-\mathrm{1}}$) was added to estimate the diatom valve concentration (DVC; Battarbee and Kneen, 1982). In addition, one drop of ammonia (NH₃) solution was added to each sample to prevent clumping of diatom valves. The homogenised sample solution was transferred onto cover slips using Battarbee cups and mounted to slides using Naphrax. A minimum of 350 diatom valves in each sample were counted along transects (mean: 391 valves) using a ZEISS AXIO Scope.A1 light microscope with a Plan-Apochromat $\mathrm{100}\times /\mathrm{1.4}$ Oil Ph3 objective at 1000× magnification. Diatom species were identified to lowest possible taxonomic level primarily using classical identification literature (Krammer and Lange-Bertalot, 1991, 1997b, a; Krammer et al., 1991; Hofmann et al., 2011), supported by regional (Russian) floras (Komarenko and Vasilyeva, 1975; Moiseeva and Nikolaev, 1992; Makarova, 2002). Diatom nomenclature and synonyms were applied using online databases such as Diatoms.org (https://www.diatoms.org; last access: 1 August 2024; Spaulding et al., 2021) and AlgaeBase (https://www.algaebase.org; last access: 1 August 2024; Guiry and Guiry, 2024) and for some species further supported with input by members of the diatom community online platform DIATOM-L (Bahls, 2015).

In addition, silicified chrysophyte Mallomonas scales were counted without further specification to calculate the Mallomonas index, which measures Mallomonas in relation to diatom cells (M $/$ D), to evaluate the degree of thermal stratification and the trophic status (Smol, 1985; Ginn et al., 2010). The ratio of planktonic-to-benthic diatoms (P $/$ B) of the most abundant species was calculated, excluding tychoplanktonic species, as an additional indicator of paleoenvironmental change. Habitat assignment followed mainly the online database Diatoms.org (https://www.diatoms.org; last access: 1 August 2024; Spaulding et al., 2021) and habitat information from Barinova et al. (2011). Furthermore, the abundances of the dominant Aulacoseira species (A. subarctica and A. ambigua) were summarised to better identify the assemblage changes in relation to the planktonic taxa Discostella stelligera.

2.4 Diatom purification process and silicon isotopes

In a previous study (Stieg et al., 2024b) which focused on the oxygen isotopes of diatoms of the sediment short core EN18232-1, a comprehensive description of the diatom isotope purification procedure of the samples is provided. The silicon isotope composition was determined from the same aliquots (Leng and Sloane, 2008). In general, the purification assessment of diatom samples includes wet chemical and physical preparation steps (Morley et al., 2004; Leng and Barker, 2006). Essentially, the cleaning procedure included the removal of organic matter with H₂O₂ (30 %), eliminating carbonates using HCl (10 %) and multiple heavy liquid separations (HLSs) with sodium polytungstate (SPT) solutions with decreasing densities (EN18232-1: 2.50–2.12 g cm⁻³) to separate the diatom valves from the heavy minerogenic fraction. The high purity of the diatom samples was supported by a contamination assessment using a JEOL M-IT500HR analytical scanning electron microscope (SEM) with an integrated energy-dispersive X-ray spectroscopy (EDS) system.

The laser fluorination method was applied using bromine pentafluoride (BrF₅) as a reagent to completely extract oxygen from the diatom SiO₂ (Clayton and Mayeda, 1963), after employing the inert gas flow dehydration method (Chapligin et al., 2010). The released silicon component combines with BrF₅ to form silicon tetrafluoride (SiF₄). The separation of SiF₄ from the bromine compounds was achieved using a cold trap, composed of a slush of ethanol cooled with liquid nitrogen (N₂) at −115 °C. SiF₄ was then transferred and trapped in Pyrex glass tubes for the subsequent silicon isotope analysis (Chapligin et al., 2010; Maier et al., 2013). For all samples (n=39) two SiF₄ extractions were performed (except two samples with only one extraction). The samples were measured 10 times on a MAT 252 IRMS at AWI Bremerhaven against a reference gas with known isotopic composition, calibrated to the reference quartz standard NBS-28. The data were processed for outliers by performing a Dixon test with a confidence level of 80 % (Dixon, 1953). The results are expressed as delta notation δ³⁰Si (${}^{\mathrm{30}}\mathrm{Si}/{}^{\mathrm{28}}\mathrm{Si}$) in per mil (‰). The calibration was done with the internal biogenic silica standard BFC with a known δ³⁰Si of +0.13 ‰ relative to NBS-28 (Leng and Sloane, 2008). Analytical precision determined by repeated analyses of the internal standard (n=12) was better than 0.07 ‰ (1σ of δ³⁰Si). A weighted mean value for each sample measurement was calculated with a standard deviation resulting from the internal or the external consistency, depending on which value is higher.

2.5 Biogeochemical proxies

The total organic carbon (TOC) and total nitrogen (TN) of 39 samples were determined to analyse historical changes in the organic matter content in the lake sediments. Freeze-dried subsamples were ground to obtain a homogeneous material. TOC was analysed using an Elementar soli TOC cube. Total carbon (TC) was computed as the sum of TOC and the total inorganic carbon (TIC), published in a previous study (Stieg et al., 2024b). Organic carbon accumulation rates (OCARs; $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$) were derived by dividing TOC by 100, multiplying it by the mass accumulation rate (MAR), and then multiplying it by 10 000 to convert to units of grams per square metre per annum ($\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$). TN was quantified by an Elementar rapid MAX N exceed. The measurement accuracy was 0.1 % for both devices, carried out at the Permafrost Carbon and Nitrogen Lab at AWI Potsdam. The TOC $/$ TN atomic ratio (C $/$ N), as an indicator of the organic matter source, was calculated by multiplying its mass ratio by 1.167, which is the ratio of the atomic weights of nitrogen and carbon (Meyers and Teranes, 2001).

The stable carbon isotope ratio (${}^{\mathrm{13}}\mathrm{C}{/}^{\mathrm{12}}\mathrm{C}$) of all freeze-dried and milled samples was examined, to gain additional information on past productivity and organic matter sources. Prior to isotope measurement, inorganic carbon was removed by a hydrochloric acid treatment (HCl, 1.3 mol L⁻¹). Carbon isotopes were measured by a ThermoFisher Scientific Delta V Advantage gas mass spectrometer equipped with a FLASH elemental analyser EA 2000, a CONFLO IV gas-mixing system and a MA200R autosampler at the ISOLAB of AWI Potsdam. Isotope ratios are given relative to laboratory standards of known isotopic compositions and relative to the Vienna Pee Dee Belemnite (VPDB) reference standard. Nitrogen isotope ratios (${}^{\mathrm{15}}\mathrm{N}{/}^{\mathrm{14}}\mathrm{N}$) of the 39 samples were analysed using the same gas mass spectrometer and expressed relative to the standard atmospheric nitrogen ratio (air). Both isotope ratios are given as delta notation (δ¹³C, δ¹⁵N) in per mil (‰), with 1σ standard errors better than ±0.15 ‰.

2.6 Data processing and statistical analyses

For data processing and statistical analyses, we used the R environment (R Core Team, 2024). Diatom species richness (alpha diversity) and evenness (effective richness) were assessed by calculating Hill's N0 and N2 (Hill, 1973), with a rarefaction of the dataset based on the count's minimum (n=352) to prevent biases of varying sample sizes, utilising the R package “vegan” (Oksanen et al., 2022). By categorising the diatom counts into dissolved and pristine valves, the F index was calculated to assess the diatom valve dissolution (Ryves et al., 2001; after Flower and Likhoshway, 1993). The F index has a total range between 0 and 1, whereby 1 equals a perfect preservation.

For further processing, the total counts of diatom species were converted to percentage data. With the R package “rioja” (Juggins, 2022) constrained incremental sums-of-squares clustering (CONISS) was performed to define diatom zones in the record, after square-root transformation and Euclidean dissimilarity calculation of the percentage data of all involved species using the “vegan” package (Oksanen et al., 2022), to ensure all species contribute meaningfully to the analysis, regardless of their abundance (Grimm, 1987). The diatom relative abundances and CONISS clustering are presented in a stratigraphic diagram using the function “strat.plot” of the “rioja” package (Juggins, 2022). The stratigraphy plot displays the predominant species only, which occur at least twice (≥2 samples) with a minimum abundance of 2 %.

The length of gradient variation in the percentage data was determined by a detrended correspondence analysis (DCA), calculated in standard deviation units (SD=1.03) using the “decorana” function of the R package “vegan” (Oksanen et al., 2022). Following the thresholds discussed in Birks (2010), a principal component analysis (PCA) was conducted on the most abundant species (≥2 % in ≥2 samples) after square-root transformation. Additionally, data of diatom isotopes and biogeochemical proxies for the same sediment samples (δ¹⁸O_diatom, δ³⁰Si_diatom, TOC, TN, C $/$ N, OCAR, δ¹³C, δ¹⁵N, mercury, diatom valve concentration, diatom accumulation rates, M $/$ D) were projected onto the PCA ordination plot of diatom data using the “envfit” function (Oksanen et al., 2022), enabling the visualisation within the taxonomical and ecological context.

Diatom accumulation rates (DARs) were calculated by multiplying the diatom valve concentration (DVC in 10⁷ valves g⁻¹) by MARs ($\mathrm{g}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$) according to Birks (2010) and converting it to ${\mathrm{10}}^{\mathrm{9}}\mathrm{valves}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$.

3 Results

3.1 Diatom assemblages and corresponding indices

Diatoms occur in all 39 samples of the short core EN18232-1. In total, 102 species were identified, with 12 species occurring with more than 2 % in ≥2 samples (Fig. 2). The valve dissolution F index, not displayed here, reveals overall a good valve preservation within the short core (mean: 0.84; min: 0.71; max: 0.92). The record is dominated by the planktonic-centric genus Aulacoseira, mainly represented by A. subarctica and A. ambigua, followed by Tabellaria flocculosa with a mean abundance of 8.9 %. According to the cluster analysis, we divided the diatom assemblage of EN18232-1 into five diatom zones (Fig. 2). The following percentage values refer to the mean abundance of the individual diatom species within the specific diatom zone. Only the most abundant species (n=12) are considered.

Figure 2 Relative abundances of the most dominant diatom species (≥2 % in ≥2 samples) of the short core EN18232-1 vs. sample depth and corresponding mean ages. Diatom zones are based on CONISS clustering shown on the right. Habitat preferences are shown at the top.

Diatom zone 1 (39–35 cm; ca. 1790–1815 CE)

The first zone of the diatom record is dominated by the planktonic species A. subarctica (43.8 %), followed by the facultatively planktonic T. flocculosa (13.2 %) and the planktonic A. ambigua (12.7 %). Besides the two Aulacoseira species, the planktonic species Asterionella formosa (7.0 %), D. stelligera (1.2 %), Fragilaria radians (0.4 %) and Stephanodiscus alpinus (2.5 %) are also represented, in which the latter appears almost only in this diatom zone. Facultative planktonic diatom species are further represented by Staurosirella pinnata (3.3 %) and Staurosira venter, whereby the latter has only low abundance with a maximum of 1 %. Overall benthic diatoms are the least represented, especially by Sellaphora seminulum (3.2 %), followed by Psammothidium subatomoides and Nitzschia cf. perminuta.

Diatom zone 2 (34–28 cm; ca. 1815–1850 CE)

Diatom zone 2 is characterised by the two most abundant species, A. subarctica (54.8 %) and A. ambigua (11.3 %), followed by T. flocculosa (10.1 %). A. formosa (5.4 %) show a slight decrease. Comparable to diatom zone 1, S. pinnata (3.2 %) represent facultative planktonic diatoms, and benthic diatoms are mainly represented by S. seminulum (3.8 %).

Diatom zone 3 (27–17 cm; ca. 1850–1910 CE)

In addition to the almost continuous dominance of A. subarctica (52.7 %) and A. ambigua (16.4 %), T. flocculosa shows an increase (11.4 %) and a peak around 1900 CE at the transition to zone 4 (max: 32 %), while both Aulacoseira taxa decline at this point. S. pinnata (3.5 %) and S. seminulum (2.2 %) also show a slight increase towards zone 4. A. formosa (2.6 %) continues to decline, while the facultative planktonic species S. venter (1.2 %) appears more frequently in this zone.

Diatom zone 4 (16–7 cm; ca. 1910–1970 CE)

Dominant planktonic species are represented by A. subarctica (55.7 %) and A. ambigua (18.1 %), both reaching a maximum in this zone, while A. formosa (2.8 %) stays at a reduced level. Facultative planktonic diatoms are represented rather equally by T. flocculosa (5.9 %) and S. pinnata (5.1 %), the latter reaching a maximum towards zone 5 (max: 12.0 %). Benthic diatoms S. seminulum and P. subatomoides show only minor occurrence.

Diatom zone 5 (6–1 cm; ca. 1970–2015 CE)

The uppermost diatom zone is characterised by a decrease in the two dominant Aulacoseira species (A. subarctica: 47.9 %; A. ambigua: 11.5 %) and a clear increase in the planktonic D. stelligera (6.5 %) as well as in A. formosa (5.5 %) and F. radians (2.6 %). The occurrence of S. pinnata (6.5 %) remains elevated in this zone. The abundance of T. flocculosa (4.1 %) declines towards the top to its lowest values in the record. Benthic diatoms are still underrepresented mainly based on S. seminulum (3.4 %).

Indices, derived from the diatom dataset, are presented in Fig. 3. The diatom valve concentration (DVC) increases towards the top of the core and ranges from $\mathrm{82.7}-\mathrm{533.8}\times {\mathrm{10}}^{\mathrm{7}}\mathrm{valves}{\mathrm{g}}^{-\mathrm{1}}$, with a mean of $\mathrm{212.2}\times {\mathrm{10}}^{\mathrm{7}}\mathrm{valves}{\mathrm{g}}^{-\mathrm{1}}$. This results in a mean diatom valve accumulation rate (DAR) of $\mathrm{464.1}\times {\mathrm{10}}^{\mathrm{9}}\mathrm{valves}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$ ($\mathrm{200.98}-\mathrm{1529.6}\times {\mathrm{10}}^{\mathrm{9}}\mathrm{valves}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$). The rarefied species richness, Hill's N0, has a mean of 21.6 and varies between 13.9 and 27.7. The evenness (effective richness), Hill's N2, results in a mean of 3.3 and spans 2.1–6.0. The combined abundance of A. ambigua and A. subarctica has a mean of 66.7 % and has a similar curve to the ratio of planktonic-to-benthic diatoms (P $/$ B). Downcore the P $/$ B ratio fluctuates along the mean value of 24.6 before it declines to values below the mean after ca. 1950 CE. The Mallomonas index (M $/$ D) has a mean of 3.5 and ranges from 0 to a clear maximum of 21.9 in the topmost sample (Fig. 3).

Figure 3 Statistical diatom indices of the short core EN18232-1. Diatom valve concentration (DVC), diatom accumulation rate (DAR), diatom species richness (Hill's N0) and evenness (Hill's N2), abundances of Aulacoseira species (A. ambigua and A. subarctica), planktonic-to-benthic species ratio (P $/$ B) of most abundant species excluding tychoplanktonic species, main axis sample scores of principal component analysis (PC1, PC2), and the Mallomonas index (M $/$ D). Vertical lines indicate the mean values. Horizontal lines indicate diatom zones (DZs), and the dashed line marks the time before and after 1950 CE.

The biplot of sample scores on the first and second axes of the principal component analysis (PCA) of diatom species data explains nearly 45 % of the total variance, whereby PC1 explains 25.7 % and PC2 explains 19.2 % (Fig. 4). TOC, TN and the Mallomonas index (M $/$ D) shows an association with the diatom species S. pinnata along the negative axis of PC1 and PC2 and with diatom samples after 1950 CE. Mercury, measured at the same sample material (Stieg et al., 2024b), is affiliated with the diatom species D. stelligera and F. radians and also points towards samples after 1950 CE. In addition, DVC and DAR point in the same direction but are both somewhat shorter. The C $/$ N ratio, δ³⁰Si_diatom and δ¹³C display opposing clustering to the previous sediment proxies along positive PC1 and PC2 axis. δ¹⁵N is associated with S. venter and A. subarctica along the negative PC2 axis. OCAR is positioned along the positive PC2 axis, associated with T. flocculosa and S. alpinus. δ¹⁸O_diatom does not exhibit a strong directional association with any particular diatom species.

Figure 4 Biplot showing the results of the principal component analysis (PCA) with its first two dimensions (PC1, PC2) of the short core EN18232-1 derived from data of the most abundant diatom species (≥2 % in ≥2 samples) and corresponding environmental proxies. The percentage of explained variance for each principal component (PC) is displayed on the axis label. Coloured sample depths indicate mean ages before and after 1950 CE.

3.2 Biogeochemical proxies and δ³⁰Si_diatom

The stratigraphic profiles of TOC and TN show a very similar pattern throughout the core (Fig. 5). TOC has a mean value of 9.2 %, and TN has a mean of 0.9 %. Around 1800 CE, at the lower part of the core, both proxies are above their mean. Between ca. 1830–1940 CE, values decline and stay below their mean, including their minima at about 1875 CE (TOC=7.5 %; TN=0.7 %). At the onset of the 1950s, TOC and TN values rise above the mean again and increase nearly continuously from the 1970s onward, both reaching their total maxima in the uppermost sample (TOC=12.5 %; TN=1.2 %). The C $/$ N ratio has a mean of 12.6 and varies by ±1.4. Until the 1950s, the C $/$ N ratio fluctuates slightly between 12.5 and 13, including two maxima at around 1800 CE and 1895 CE. Around 1950 CE, the C $/$ N ratio shifts to lower values of about 12. In the two uppermost samples, the ratio slightly increases again towards 12.5 (Fig. 5). The values of OCAR fluctuate slightly around the mean value of 20.5 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$ and only fall below it from the 1950s onward until a minimum of 11.5 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$ at around 1970 CE, before rising to the maximum value of 50.6 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$ in the youngest sample (Fig. 5).

Figure 5 Biogeochemical records of the short core EN18232-1: total organic carbon (TOC), total nitrogen (TN), TOC $/$ TN atomic ratio (C $/$ N), organic carbon accumulation rate (OCAR), diatom silicon isotopes (δ³⁰Si_diatom), and stable carbon (δ¹³C) and nitrogen (δ¹⁵N) isotope values of bulk sediment samples. Vertical lines indicate the mean values. Horizontal lines indicate diatom zones (DZs), and the dashed line marks the time before and after 1950 CE.

δ³⁰Si_diatom ranges between −0.29 to +0.52 ‰, with an overall mean of +0.24 ‰. Until ∼1970 CE, δ³⁰Si_diatom values stay mainly above the mean, including the absolute maximum of +0.52 ‰ around 1900 CE. From the 1970s onward, δ³⁰Si_diatom values steadily decline and reach the absolute minimum of −0.29 ‰ in the topmost sample (Fig. 5).

δ¹³C values are mainly above the mean of −30.6 ‰ before 1950 CE, ranging from −31.1 ‰ at about 1800 CE to the absolute maximum of −29.5 ‰ at ca. 1900 CE (Fig. 5). From the 1950s onward, δ¹³C continuously decreases towards more depleted values, including the absolute minimum of −32.3 ‰ at about 2000 CE. The δ¹⁵N record shows a generally low variability and follows mainly the mean value of +2.9 ‰ and varies in total by ±0.5 ‰, with slightly enriched values between approx. 1930 and 1970 CE and a subsequent minor decline (Fig. 5).

4 Discussion

4.1 Limno-ecological interpretation of Lake Khamra diatom assemblages

Owing to the deep central part of Lake Khamra (maximum depth: 22.3 m; Fig. 1c), planktonic diatom species have been prominent over the last ∼220 years, as reflected by their high abundance within the short core EN18232-1 (Fig. 2). The planktonic diatom genera Aulacoseira, mainly represented by A. subarctica and A. ambigua, is the most abundant in the sediment record and common in deep, cold, predominantly unstratified lakes in the Arctic tundra (Round et al., 1990; Laing and Smol, 2003). Aulacoseira taxa rely on adequate turbulence to stay suspended within the water column (Kilham et al., 1996; Gibson et al., 2003). Thermal stratification or a perennial lake ice cover can weaken the turbulence and lead to a settlement of heavier diatom cells (Gibson et al., 2003). Most Aulacoseira species develop dense resting cells as part of their life cycle, which must be resuspended from the sediment into the water column to sustain planktonic populations (Round et al., 1990). Temperature measurements taken under the ice in March 2020 (Biskaborn et al., 2021a) indicate very cold water temperatures at Lake Khamra (mean: 2.5 °C) and a well-mixed water column (Stieg et al., 2024b). Given that A. subarctica predominates throughout the record (Fig. 2), it suggests that the growth conditions at Lake Khamra generally favour this species and that a seasonal ice cover does not necessarily reduce its abundance. Furthermore, the snow-rich region and the associated inflow of cold snowmelt water (Stieg et al., 2024b) presumably lead to sufficient turbulence and nutrient upwelling in spring, promoting sustained circulation beneficial for A. subarctica (Horn et al., 2011). However, other factors, including atmospheric–surface temperature gradients and wind-induced mixing, are also likely relevant for sufficient lake overturn (Smol and Douglas, 2007; Winder and Sommer, 2012).

Facultative planktonic species, mainly epiphytic (attached to plants), are less abundant and primarily characterised by T. flocculosa (reaching up to 32 %) and S. pinnata (up to 12 %; Fig. 2). Besides the deep central part, Lake Khamra has a shallow littoral area (Fig. 1c), where macrophytes can provide a habitat for epiphytic diatom taxa. Purely benthic diatoms (N. cf. perminuta, P. subatomoides, S. seminulum) are of low abundance at Lake Khamra and reach a maximum of 3 % to 9 %. Light limitations due to the lake depth and annual ice cover, cold water temperatures, and the availability of nutrients likely regulate the abundance of benthic taxa.

The diatom assemblage of the short core EN18232-1 characterises a deep, open freshwater lake in the boreal forest biome with oligotrophic conditions over the last 220 years. A DAR consistently above the minimum of $\mathrm{201}\times {\mathrm{10}}^{\mathrm{9}}\mathrm{valves}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$ indicates continuous diatom accumulation, and good preservation environments are reflected by an overall high F index (mean: 0.84). These factors, along with the sub-decadal resolution of the record, make Lake Khamra well-suited for paleoecological studies.

4.2 Linking diatom assemblage shifts to hydroclimatic variability

Based on cluster analysis, the diatom assemblage of the short core was subdivided into five diatom zones (Fig. 2). The following discussion focuses on major changes within these zones. In diatom zones 1–4 (pre-1970s), assemblage shifts occur mainly in diatom zone 1 ($\sim \mathrm{1790}-\mathrm{1815}\mathrm{CE}$), followed by a rather stable phase until the 1970s.

From the onset of the record until ∼1815 CE, the dominance of the two planktonic diatom species A. subarctica and A. ambigua was lower than in the subsequent three diatom zones and remains similarly low as observed in the most recent zone after ∼1970 CE (Fig. 3). A hydroclimatic reconstruction from Lake Khamra, based on a δ¹⁸O_diatom record of the same sediment samples (Stieg et al., 2024b), indicates a rather dry and cool period at the onset of the record ($\sim \mathrm{1790}-\mathrm{1830}\mathrm{CE}$). This is reflected in elevated δ¹⁸O_diatom values, suggesting reduced winter precipitation (Fig. 6). It is assumed that less precipitation in winter, accumulated as snow, leads to reduced snowmelt water inflow into the lake. This could lead to reduced water turbulence in spring and hence less favourable conditions for Aulacoseira to bloom. Both species richness indices (Hill's N0 and N2) are notably higher during Aulacoseira reduction, showing a clear distinction between the lowest diatom zone and the subsequent ones. The decline in the dominant Aulacoseira species likely reduces competitive pressure, allowing other diatom species to become established and leading to a more diverse diatom assemblage. Noteworthy for the planktonic species is S. alpinus, which is primarily present only in the lowermost zone of the diatom assemblage. S. alpinus appears in oligotrophic lakes with very low water temperatures (Häkansson and Kling, 1989; Krammer et al., 1991), and this aligns with the assumed cold conditions at the onset of the record. In terms of nutrients, Stephanodiscus has a high demand for phosphorus, although it can still grow under low silica concentrations and limited light (Kilham et al., 1986). A reduced meltwater inflow due to overall relatively dry winter conditions could have led to both decreased overturn and a lower input of nutrients, including silica (Si). Aulacoseira requires high Si availability for their highly silicified frustules (Laing and Smol, 2003), so a reduction in Si input, combined with weaker turbulence, could have favoured other planktonic species like S. alpinus as stronger competitors. Overall, however, DARs in the lowest diatom zone are rather low, gradually increasing over time, likely due to the initially cool conditions.

After the dry and cold period at the onset of the record, both Aulacoseira species show a clear increase in abundance until the 1940s, reaching their highest abundance of the entire record (Fig. 3). The dominance is accompanied by a decline in effective species richness, as indicated by decreasing Hill's N2 (Fig. 3). The δ¹⁸O_diatom hydroclimatic reconstruction suggests a shift from dry to wet conditions around 1830 CE, followed by a prolonged rather stable period in which the lake water was primarily influenced by snowmelt until ∼1930 CE (Stieg et al., 2024b). Constant hydroclimatic conditions and sufficient snowmelt input presumably provided lake turbulence in spring and nutrient upwelling, beneficial for A. subarctica (Horn et al., 2011).

Around 1935 CE, Aulacoseira reached a maximum in abundance (∼85 %) but subsequently declined to a minimum of ∼50 % at ca. 1970 CE (Fig. 3), reflecting a period of high variability. This corresponds to a phase of hydroclimatic instability, as indicated by meteorological data from the nearest weather station in Vitim, which suggest a prolonged annual precipitation deficit (Fig. 6), particularly in early winter months around the 1950s. This further aligns with a reconstructed dry anomaly in the 1950s based on the δ¹⁸O_diatom values (details given in Stieg et al., 2024b) and coincides with increased wildfire activity in the vicinity of the lake, as inferred from a charcoal record from Lake Khamra (Glückler et al., 2021). A reduction in spring overturn, caused by reduced snowmelt water supply, seems a reasonable explanation for the decline in Aulacoseira. Comparable to the transition from diatom zone 1 to 2, we observe an increase in diatom species richness and especially in evenness between ca. 1950 and 1970 CE (Hill numbers, Hill's N2 peak in the 1970s; Fig. 3). This pattern coincides with the decline in the dominant Aulacoseira taxa, likely allowing other diatom species, like facultative planktonic T. flocculosa or S. pinnata, to compete more successfully.

4.2.1 Diatom shift at ca. 1970 CE: a climate warming signal

The most prominent shift in diatom assemblage occurs around 1970 CE, corresponding to the transition between diatom zones 4 and 5 (Fig. 2). Diatom zone 5 ($\sim \mathrm{1970}-\mathrm{2015}\mathrm{CE}$) is characterised by a substantial increase in the planktonic D. stelligera, rising from low relative abundances up to 8.5 %, while both Aulacoseira species show a minimum (Fig. 6). D. stelligera (basionym Cyclotella stelligera; Guiry and Guiry, 2024) belongs to diatoms of the Cyclotella sensu lato taxa (Saros and Anderson, 2015). The shift in diatom assemblage characterised by an increase in small cyclotelloid planktonic diatoms like the Cyclotella stelligera complex and a simultaneous decrease in Aulacoseira taxa and/or small benthic fragilariod taxa is revealed in other paleolimnological studies from the Northern Hemisphere (Rühland et al., 2003, 2008; Smol et al., 2005). Diatom zone 5 corresponds with a hydroclimatic phase of the δ¹⁸O_diatom record, indicating increased winter precipitation and more snowmelt water supply to Lake Khamra, supported by meteorological data (Fig. 6) (Stieg et al., 2024b). Despite increased water turbulence from inflow, which usually benefits Aulacoseira, its observed decrease is unexpected. Therefore, we infer that an additional factor is responsible for the decline. From the 1970s onward regional annual temperatures start to increase (anomalies of annual mean temperature, Vitim; Fig. 6). Rising temperatures in the context of global warming likely reduce the ice-cover duration and increase the possibility of thermal stratification of the lake water in the ice-free period, which allow D. stelligera to thrive (Smol et al., 2005; Rühland et al., 2015). Increased summer stratification is less favourable for Aulacoseira (Gibson et al., 2003; Rühland et al., 2008), although it remains dominant at Lake Khamra (>50 %; Fig. 3).

Additionally, we observe a rapid increase in the silicified chrysophyte Mallomonas scales since the 1990s, inferred from the Mallomonas index (M $/$ D) (Fig. 3). Chrysophytes are common in oligotrophic environments (Smol, 1985). Their motility, enabled by flagella, allows them to thrive in stratified lakes by maintaining their position in the photic zone. This gives them an advantage over non-motile, colonial diatoms such as Aulacoseira (Ginn et al., 2010; Mushet et al., 2017). The observed increase in chrysophytes at Lake Khamra suggests changes in the lake's mixing regime. Further it provides evidence for a likely longer ice-free period and enhanced thermal stratification during summer months in recent decades. Similar increases in scaled chrysophytes have been reported in other lake systems, associated with climate warming and increased thermal stability (Paterson et al., 2004; Ginn et al., 2010; Favot et al., 2024).

An increase in annual temperature and a shorter lake ice period extend the diatom growing season and change light and nutrient availability, favouring small planktonic diatoms (Rühland et al., 2015). Observed elevated diatom species richness values (Hill's N0; Fig. 3) are congruent with reported increased diatom community complexity and diversity during prolonged ice-free periods in Arctic lakes (Douglas and Smol, 2010; Rühland et al., 2015). Different studies have reported a rise in pennate (tychoplanktonic) planktonic Fragilaria populations as a response to climate warming (Rühland et al., 2013; Michelutti et al., 2015; Sochuliakova et al., 2018). S. pinnata and F. radians both group together with D. stelligera in the PCA biplot (Fig. 4), and F. radians also occur primarily in the uppermost diatom zone (Fig. 2), underlining a possible common reaction to warming at Lake Khamra.

In many Arctic lakes, comparable shifts in the diatom community occurred as early as 1850 CE (Smol et al., 2005), while the change at Lake Khamra occurred about 120 years later. The timing discrepancy regarding the onset of increase in small Cyclotella sensu lato species aligns with studied temperate regions in North America and Europe (Rühland et al., 2008), which also show a change from 1970 CE onward. Yakutian lakes are rare in circumpolar studies (Smol et al., 2005; Rühland et al., 2008). A similar shift in diatom composition has been observed at Lake Bolshoe Toko, a comparable remote boreal lake in eastern Siberia, after 1850 CE, interpreted as a response to climate warming (Biskaborn et al., 2021b). However, the much lower stratigraphic resolution at Lake Bolshoe Toko introduces age uncertainties, making direct comparison with the Lake Khamra record difficult. The observed age offset in diatom shifts could refer to different climate-related drivers. Further diatom assemblage shifts are observed at the south basin of Lake Baikal at ca. 1970 CE (Roberts et al., 2018), as well as at pristine Arctic and sub-Arctic lakes in the northern Urals with most distinct shifts after 1970 CE (Solovieva et al., 2008). Our data align with the later onset of diatom assemblage changes at lower latitudes compared to the earlier changes observed in Arctic regions, which is attributed to an earlier warming and a quicker response of Arctic ecosystems compared to lower latitudes (Smol et al., 2005; Smol and Douglas, 2007).

The shifts in diatom assemblage delineating zones 4 and 5 as well as 1 and 2 have similarities (Aulacoseira decline while planktonic cyclotelloid taxa such as S. alpinus or D. stelligera increase and diatom richness increase). However, we also see clear differences (cold habitat preferences from S. alpinus vs. warm and stratified preferences of D. stelligera, increase in F.radians and S. pinnata in recent years only) which suggests that the increase in air temperatures in recent decades notably influenced the lake ecosystem and its diatom assemblages.

Figure 6 Diatom assemblage changes and selected biogeochemical proxies in relation to hydroclimatic variability: Aulacoseira vs. D.stelligera abundances of the short core EN18232-1, diatom accumulation rates (DARs), mass accumulation rates (MARs) and carbon accumulation rates (OCARs), δ³⁰Si_diatom record, δ¹⁸O_diatom record from the same short core of Lake Khamra (Stieg et al., 2024b), and anomalies of annual mean temperature and annual precipitation from the weather station in Vitim (59.45° N, 112.58° E; 186 m a.s.l.). Data from ECA (European Climate Assessment), accessible via the KNMI Climate Explorer https://climexp.knmi.nl (last access: 15 July 2024) (Klein Tank et al., 2002). Vertical lines indicate the individual mean values. Horizontal lines indicate diatom zones (DZs), and the dashed line marks the time before and after 1950 CE.

4.3 Biogeochemical proxies as indicators of bioproductivity and organic matter sources at Lake Khamra

4.3.1 Biogeochemical shifts after ca. 1950 CE

Carbon and nitrogen indicators (TOC, TN, C $/$ N, OCAR, δ¹³C; Fig. 5) show a notable shift at ∼1950 CE, highlighting substantial limnological changes that precede the major shifts in diatom communities observed starting at ca. 1970 CE (onset of diatom zone 5). We propose that hydroclimatic variability substantially influences limnological conditions at Lake Khamra, with some diatom species responding immediately. However, the most ecologically important shift in the diatom assemblage appears to occur at ∼1970 CE and is more likely driven by rising temperatures and their associated limnological effects, as discussed above.

At ∼1950 CE, the δ¹⁸O_diatom values reach a maximum (Fig. 6), interpreted as a dry anomaly with a reduced input of isotopically light winter precipitation into the lake (Stieg et al., 2024b). This reconstructed anomaly aligns with the overall dry period inferred from persistently below-average annual precipitation starting from the 1930s, at the onset of meteorological recordings, and continuing until the 1970s (Fig. 6; details in Stieg et al., 2024b). Additionally, a charcoal record from the same lake identifies enhanced fire activity during the 1950s (Glückler et al., 2021), coinciding with the δ¹⁸O_diatom maximum, providing further evidence of this dry anomaly. At the same time, MARs drop (Fig. 6), which could be linked to a reduction in erosional input, further supported by the flattening of the age–depth relationship from the 1950s onward (Fig. A1), indicating lower sedimentation rates. The concurrent decrease in OCAR to a minimum of 11.5 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{a}}^{-\mathrm{1}}$ between ca. 1950 and 1970 (Figs. 5 and 6) suggests lower erosional input of organic material, a reduced primary productivity or a combination of both. In addition to the dry phase, declining annual temperatures, reaching a minimum between ∼1950 and 1970 CE (Fig. 6), indicate simultaneously cool conditions. While Aulacoseira remains dominant, it decreases especially between the 1950s and the 1970s. DARs decline after the peak in the 1930s, indicating rather a reduction in primary productivity around the 1950s. This is likely linked to cold annual temperature anomalies and less nutrient refreshment due to reduced turbulence and inflow, as Aulacoseira begins to decline. In contrast, TOC and TN, which represent the fraction of organic matter and are closely related to primary productivity (Meyers and Teranes, 2001), both show a peak between ca. 1950 and 1970 CE (Fig. 5), indicating higher bioproductivity. It is plausible that low sedimentation rates could result in a relative increase in TOC and TN concentrations within the sediment, explaining the contrasting tendency.

During the dry period between the 1950s and 1970s, the erosional input from the catchment was probably reduced. However, a differentiation between autochthonous (in-lake) or allochthonous origin of organic material since then remains indistinct. The C $/$ N atomic ratio and δ¹³C help to identify the organic matter origin (Meyers, 2009). C $/$ N drops from ca. 13 to a slightly lower level of around 12 from the 1950s onward, reflecting a mixed signal in both cases. According to the literature, lacustrine algae have a C $/$ N ratio ranging between 4 and 10, whereas organic matter from land plants have a ratio>20 (Meyers and Teranes, 2001). Also together with δ¹³C, a differentiation still remains inexplicit (Fig. 7), as both algae and C₃ land plants have similar δ¹³C signatures, making them indistinguishable (Meyers, 2009). The δ¹⁵N data underline a mixed organic origin, as the value of 3 ‰ lies in between land plants (around 0 ‰ ) and freshwater algae (5‰–10‰; Meyers, 2009). However, the rather low C $/$ N values tend towards an autochthonous organic matter source from the 1950s onward, underlining a reduced inflow from the catchment.

Figure 7 C $/$ N and δ¹³C plot showing sediment samples of the short core EN18232-1. In general, lacustrine algae have a C $/$ N ratio between 4–10, and C₃ land plants have a ratio of >20 (Meyers and Teranes, 2001). Figure is based on the values given by Meyers and Teranes (2001).

After ∼1970 CE, during diatom zone 5, annual precipitation in the region increases (Fig. 6) and offers the possibility for an increased erosional input from the catchment. Simultaneously, we observe an increase in both TOC and TN (Fig. 5), pointing towards a higher bioproductivity (Meyers and Teranes, 2001). Since OCAR increases at the same time (Fig. 5), it supports an enhanced accumulation of TOC from ca. 1970 CE onward. However, the C $/$ N ratio in combination with δ¹³C still mirrors a mixed signal of its origin (Fig. 7). The DAR lags slightly behind and increases after the 1980s (Fig. 6), indicating a tendency towards enhanced diatom productivity likely related to warming (Douglas and Smol, 2010; Biskaborn et al., 2012, 2021b).

The continuous decrease in δ¹³C especially since the 1950s (Fig. 5) initially appears to contradict the interpretation of increased bioproductivity in the upper part of the record, even though algae and land plants have similar signatures (Meyers, 2009). In general, phytoplankton preferentially incorporates lighter ¹²C than ¹³C, which enriches the remaining reservoir over time and leads to increasing δ¹³C values of organic matter during enhanced productivity (Meyers, 2009). The observed δ¹³C decline might, at least partly, be related to fossil fuel combustion and biomass burning, which release CO₂ depleted in ¹³C into the atmosphere, known as the Suess effect (Keeling, 1979). As the CO₂ in the atmosphere and in the water tend to reach equilibrium, this also affects the δ¹³C values of lake sediment organic matter (Verburg, 2007), as observed e.g. in a lake study in west Greenland (Stevenson et al., 2021). Most of the decline in δ¹³C of about 1.7 ‰ at Lake Khamra could, thus, be explained by the decline in atmospheric δ¹³C, which has accelerated particularly after 1950 and has fallen by about 1 ‰ since then (Verburg, 2007). Aside from that, an increased algae productivity might not contradict the decreasing δ¹³C values, as light CO₂ can be constantly assimilated during increased growth, as observed in an Arctic lake (Jiang et al., 2011).

As with δ¹³C, δ³⁰Si_diatom also shows a decline but slightly later at around 1970 CE (Fig. 5 and 6), coinciding with the shift in the diatom assemblages between diatom zones 4 and 5. We assume that the decrease in δ³⁰Si_diatom is influenced by multiple factors, as the interpretation of silicon isotopes is complex and involves various processes (Sutton et al., 2018). Overall, the δ³⁰Si_diatom values are low (mean: +0.24 ‰), comparable to e.g. a study at a groundwater-dominated Swedish lake (Zahajská et al., 2021a). Groundwater inflow could lead to low isotope values of dissolved silica (DSi, δ³⁰Si_DSi) and hence to low δ³⁰Si_diatom (Opfergelt et al., 2011; Frings et al., 2016; Zahajská et al., 2021b, a). Unfortunately, there is no information on lake water δ³⁰Si_DSi or groundwater inflow at Lake Khamra. It is likely that groundwater input at Lake Khamra contributes to the generally low δ³⁰Si_diatom. However, a sudden change in groundwater seems unlikely and cannot explain the distinct decrease at around 1970 CE.

The fractionation of silicon isotopes occurs through both inorganic chemical weathering and biological activities, mainly diatom growth in marine and terrestrial environments (Opfergelt et al., 2011). In marine environments, diatom silicon utilisation is linked with the silicon isotope composition (De la Rocha et al., 2000; Varela et al., 2004; Cardinal et al., 2005; De la Rocha, 2006), whereas studies in lacustrine environments are sparse. There are studies linking diatom blooms in lakes with silicon isotope fractionation by seasonal in situ measurements to validate δ³⁰Si_diatom as a bioproductivity proxy (Alleman et al., 2005; Opfergelt et al., 2011; Panizzo et al., 2016). In this context, δ³⁰Si_diatom is used as a paleoenvironmental proxy in only a few paleolimnological studies in Asia, for example at Lake El'gygytgyn in north-east Siberia (Swann et al., 2010) or in south China (Chen et al., 2012).

Diatoms preferably incorporate light ²⁸Si instead of heavier ³⁰Si (and ²⁹Si) during biomineralisation (De la Rocha et al., 1997; Leng et al., 2009). Assuming this fractionation dominates the isotopic signal of the diatoms at Lake Khamra, the decrease in δ³⁰Si_diatom to values of −0.53 ‰ in recent decades contradicts the discussed tendency towards higher algae productivity in diatom zone 5. Higher bioproductivity increases DSi utilisation and would result in an enrichment of both the DSi signature of the reservoir (lake water) and the diatom frustules therein.

In unlimited lake systems, such as Lake Khamra, alteration of DSi sources can influence δ³⁰Si_diatom (Zahajská et al., 2021b) and may overprint the isotopic effect of bioproductivity. We suggest that substantial changes in the catchment have led to isotopic changes in the reservoir (lake water), which influenced the δ³⁰Si_diatom record at Lake Khamra. From ∼1970 CE onward, both temperature and precipitation increase in the region (Fig. 6), forcing enhanced weathering in the catchment and e.g. leading to a higher nutrient inflow to the lake (Douglas and Smol, 2010), as supported by the increase in OCAR (Fig. 5). In addition, the recent wet and warm period at Lake Khamra follows an earlier cold and dry period around the 1950s with increased wildfire activity in the region (Glückler et al., 2021; Stieg et al., 2024b). Fires in the catchment area have likely altered vegetation and soil conditions, consequently affecting the source of dissolved silica. Soils and terrestrial plants, which tend to have low δ³⁰Si_DSi values (Frings et al., 2016; Sun et al., 2018), may be subject to increased erosion following the (partial) loss of the vegetation cover. Additionally, Lake Khamra lies within a sporadic permafrost zone (Fedorov et al., 2018). Observed warming could increase the depth of the active layer, which, as it thaws, would also preferably transport low δ³⁰Si_DSi from litter and soil material into the lake.

Besides the change in catchment supply, an alteration of species composition driven by variations in summer temperatures and the duration of ice-free periods (Zahajská et al., 2021b) could further influence δ³⁰Si_diatom. It remains uncertain if a species effect leads to a different silicon isotope fractionation. While Sutton et al. (2013) provided evidence of such effects in marine diatoms, Schmidtbauer et al. (2022) observed varying δ³⁰Si_diatom values in lacustrine samples, which were associated with different preferred habitat regimes. We hypothesise that the δ³⁰Si_diatom decrease at Lake Khamra could partially mirror the major shift in diatom assemblage from the 1970s onward, from heavy and compact Aulacoseira taxa to less silicified D. stelligera. Overall, interpreting δ³⁰Si_diatom in lacustrine systems is more complex than in marine environments. At Lake Khamra, δ³⁰Si_diatom appears to be influenced by multiple factors, including shifts in diatom assemblages and catchment alterations that affect silica sources and utilisation.

4.3.2 Biogeochemical trends before ca. 1950 CE

Before ∼1950 CE we observe a prolonged, rather stable period in most biogeochemical proxies (TOC, TN, C $/$ N, OCAR, δ¹³C, δ³⁰Si_diatom, δ¹⁵N; Fig. 5) characterised by a dominance of Aulacoseira in diatom zones 2, 3 and 4. The slightly elevated C $/$ N ratio of around 13 in combination with elevated δ¹³C values (Fig. 7) still indicates a mixed source signal of organic material, as outlined above. Compared to the period after ∼1950 CE, the older section of the record shows a tendency towards a more allochthonous origin of the organic material (Meyers, 2009), which could indicate either lower lake-internal productivity or increased erosional inflow of organic matter. Only a few species dominated the record at this time (lowest Hill's N2; Fig. 3), mainly represented by the two Aulacoseira taxa. TOC and TN values show their lowest concentration between ∼1830 and 1950 CE (Fig. 5), suggesting a reduced bioproductivity (Meyers and Teranes, 2001). OCAR values are less indicative and remain relatively constant until the 1950s (Fig. 5), with DAR showing similar stability, though with some dips indicating periods of lower accumulation (Fig. 6). These trends point to a relatively constant or slightly decreased bioproductivity. As determined from the δ¹⁸O_diatom record (Fig. 6) (Stieg et al., 2024b), the hydroclimatic conditions were rather stable from ∼1830 CE onward, likely with sufficient inflow from the catchment by snowmelt. This rather wet period could explain the tendency towards an allochthonous origin of organic material.

Instead, elevated δ³⁰Si_diatom would suggest enhanced bioproductivity based on the classical interpretation of this proxy (De la Rocha et al., 1997). Aulacoseira has a high demand of dissolved silica for their heavy and thick frustules, which could have led to a moderate δ³⁰Si enrichment of the reservoir during phases of diatom bloom. Furthermore, we assume that, in contrast to the post-1950s period, the catchment was not yet affected by wildfire disturbances, leading to a difference in the isotopic signature of the pre- and post-1950s. Less influence of isotopic light δ³⁰Si plant material from soil erosion potentially resulted in a higher δ³⁰Si level of dissolved silica, which is reflected in the overall elevated δ³⁰Si_diatom (Fig. 6). This highlights the complexity of δ³⁰Si_diatom as a proxy, linked likely to both in-lake biogeochemical processes and catchment dynamics at Lake Khamra.

4.4 Anthropogenic influence on Lake Khamra linked to atmospheric pollution

There are indications of human pollution at Lake Khamra despite its remote location and the absence of industrial activity in its vicinity. The decline in δ¹³C values since the 1950s (Fig. 5) is likely linked to the release of lighter carbon isotopes from fossil fuel combustion (Verburg, 2007), as reflected in the lake's sediments.

A previous study (Stieg et al., 2024b) identified a marked increase in mercury concentrations in the same lake sediment samples, tripling since ∼1930 CE, with fluxes rising more than 4-fold from the 1990s onward. Mercury in the lake can originate from both natural sources, such as permafrost (Rutkowski et al., 2021) and biomass burning (Burke et al., 2010; Driscoll et al., 2013), as well as from anthropogenic sources, such as emissions (Wang et al., 2004; Streets et al., 2011). The observed increase at Lake Khamra is at least partly related to industrial growth and associated atmospheric mercury emissions in Asia and Russia (Pacyna et al., 2016; Sundseth et al., 2017; Eckhardt et al., 2023), comparable to mercury accumulation rates observed at Lake Baikal post-1850 CE (Roberts et al., 2020). Here, atmospheric fallout from the 1980s onward is also suggested to have likely contributed to rising mercury accumulation rates, in addition to mining activities in the vicinity of Lake Baikal (Roberts et al., 2020). Similarly, a rapid mercury increase since 1850 CE has been linked to atmospheric pollution at another remote boreal lake in eastern Siberia (Biskaborn et al., 2021b), suggesting that Lake Khamra is likely also influenced by anthropogenic mercury deposition.

In contrast to the recent marked increase in mercury levels and δ¹³C depletion, δ¹⁵N in Lake Khamra sediments has shown only minimal variation over the last ∼220 years, fluctuating by ±0.5 ‰, with a slight decrease of ∼0.3 ‰ since the 1970s (Fig. 5). Human activities, such as fossil fuel combustion and fertiliser production, are relevant sources of reactive nitrogen (Nr) that contribute to the deposition of δ¹⁵N-depleted nitrogen in lake sediments, typically in a range from 1‰–3‰ (Gruber and Galloway, 2008; Holtgrieve et al., 2011; Wolfe et al., 2013). Despite the possible influence of atmospheric pollution, no substantial δ¹⁵N depletion is observed in Lake Khamra. However, not all lakes display δ¹⁵N depletion, as nitrogen cycling in lakes is complex and factors such as nitrogen inputs, water residence time and aquatic activity play crucial roles (Meyers and Teranes, 2001; Galloway et al., 2003; Anderson et al., 2018). For example, increased abundance of A. formosa in oligotrophic alpine lakes in North America has been linked to atmospheric nitrogen deposition (Saros et al., 2005, 2010). At Lake Khamra, we observe only a slight increase in planktonic A. formosa in the most recent diatom zone 5 (Fig. 2). We suggest that this slight increase in abundance is more likely a response to climate warming and related changes in lake-water mixing and thermal stability rather than nitrogen enrichment by atmospheric deposition (Sivarajah et al., 2016). This further supports the argument that Lake Khamra is primarily influenced by recent climate warming, which is altering the lake's properties, rather than by atmospheric nitrogen deposition from human sources.

5 Conclusions

Our study of diatom assemblages and biogeochemical proxies over the last ∼220 years from Lake Khamra, eastern Siberia, reveals distinct ecosystem changes, likely driven by hydroclimatic dynamics. This research addresses a spatial gap in lake studies over shorter time periods in Yakutia.

Lake Khamra's diatom assemblage is dominated by planktonic diatom species, primarily Aulacoseira taxa. At ca. 1970 CE, we observe a major shift in the diatom assemblages, characterised mainly by an increase in planktonic D. stelligera, while the two dominant Aulacoseira taxa decrease. We attribute this change to warming, a probable earlier ice-out and a potential increase in summer thermal stratification in recent decades. Rapidly increasing chrysophyte scales (Mallomonas) since the 1990s support the warming trend and a related increasing stratification of the lake. These findings are temporally consistent with similar changes observed in temperate lake ecosystems globally.

Carbon and nitrogen proxies exhibit changes starting in the 1950s, indicating the hydroclimatic influence with a delayed shift in the diatom assemblages, which we attribute primarily to rising temperatures and associated effects.

At Lake Khamra, the interpretation of δ³⁰Si_diatom is not straightforward, as it is likely influenced by a combination of factors, including weathering, shifts in diatom assemblages and diatom productivity. We suggest that wildfire activity in the catchment in the 1950s, along with increased precipitation and temperatures since ∼1970 CE, altered silica sources by transporting isotopically light dissolved silica into the lake, which likely led to the observed decrease. Furthermore, the lower silica uptake of D. stelligera compared to Aulacoseira taxa might have further influenced this decline.

In addition to elevated mercury levels determined in a previous study, the human influence on Lake Khamra is evident from δ¹³C depletion, likely derived from fossil fuel combustion. No substantial δ¹⁵N depletion from atmospheric nitrogen deposition is observed in Lake Khamra sediments.

These findings underline the effects of both climate warming and anthropogenic air pollution on remote lake ecosystems, emphasising the need for comprehensive research to mitigate these impacts and protect natural water resources.

Appendix A

Figure A1 (a) Bayesian accumulation model based on the ²¹⁰Pb and ¹³⁷Cs dating (green lines) of the short core EN18232-1 (grey lines indicate the 2σ range; red line indicates the median). (b) Model iteration log. (c, d) Prior (green line) and posterior (grey area) distributions for accumulation rate and memory, respectively. Age–depth model was first published in Stieg et al. (2024b).