In May 2019, we surveyed the lake with a close-spaced grid of echosounding lines, with the aim to compile a morphobathymetric map of its floor and identify a suitable position for sediment coring. We used techniques and instruments like those employed by Gasperini et al. (2007) at Cheko Lake. The bathymetric measurements were performed with a 200 kHz echosounder tied to a wooden stick aside of a boat while a transducer was located at about 30 cm below the water surface. Instead of acquiring only the depth value provided by the echosounder through the serial port, we collected the entire echogram by means of a 16 bit A/D converter driven by SwanPRO acquisition software. The acquisition parameters of the echosounder signal were as follows: sample rate 0.856 µs, trace length 25 ms, and 2,140 trace samples (seven traces per second). The collected data were stored in XTF format, converted to SEGY format, consistency was checked using the Segy-change preprocessing software program (Stanghellini and Carrara 2017), and processed using SeisPrho software program (Gasperini and Stanghellini 2009) to display echograms and semi-automatically digitize the sediment-water interface. The lake-floor profiles were subsequently corrected for the vertical offset between the water surface and the transducer vertical position, and then exported in ASCII (lon, lat, z) format. The file was then processed with GMT to produce both a regular grid file and the final map. The original data were affected by spatial noise due to GPS precision being limited to about 1 m. To reduce this spatial noise, the final morphobathymetric map was compiled using a 2D median filter with a searching radius of 10 m.

During the same field survey, two short lake sediment cores, SUZ1 and SUZ3, were collected using a Kajak gravity corer (sampling tube diameter of 5.8 cm and length of 50 cm) from the deepest middle part of Suzdalevo Lake. The lengths of the obtained cores SUZ1 and SUZ3 were 42 and 46 cm, respectively. Both cores were immediately sliced in steps of 1 cm, with the uppermost 1 cm (0–1 cm) divided into two layers—0–0.5 cm and 0.5–1 cm. The color, grain size, and amount of plant macroremains (qualitative data) were recorded during slicing and during analyses of aquatic biota remains.

Both sediment cores were measured using gamma spectroscopy for the specific activity of ²¹⁰Pb, ¹³⁷Cs, and ²²⁶Ra isotopes in the Radiometry Laboratory (Institute of Geochemistry, Mineralogy and Mineral Resources Charles University, Prague, Czechia). Individual samples (1-cm-thick layers; the two uppermost layer were merged) were measured by well-in-well geometry in a SILAR® low-background anti-Compton-anti-coincidence gamma spectrometer with a specially designed 40 × 40 mm Na(Tl) well-type LEADMETER® detector with a total efficiency of 46.2% (for ²¹⁰Pb line of 47 keV) placed in a well-type guard NaI(Tl) detector 160 × 125 mm in 10-cm-low-background lead shielding (Hamrová et al., 2010). A Canberra DSA 2000 multichannel analyzer controlled by GENIE 2000 software was used to determine the specific activities of the above-mentioned isotopes. The measuring time for individual samples was 2 days, and 6 days for background. Special in-house standards with a light matrix were used for radionuclide quantification in the same geometry of the 8 ml vials. The IAEA-447 standard (moss-soil) was used as ²¹⁰Pb reference material with a result of 340 ± 6 Bq kg⁻¹ for the recommended value of 338 Bq kg⁻¹ and corrected for decay from the reference date (15/11/2009). Due to low weights of the SUZ1 samples affecting the results, we decided to measure isotope activities also for the SUZ3 core by merging four adjacent 1-cm-thick layers in a row. The output from the heavier 4-cm-merged samples (∼4 g of dry sediment) was expected to be more accurate. The age of individual layers was obtained using the Constant Rate of Supply (CRS), the Constant Initial Concentration (CIC), and the Constant Flux Constant Sedimentation rate (CFCS) age-depth models that are implemented in an R package “serac” (Bruel and Sabatier 2020 and references therein). These alternative models were used because none of them can be considered as a universal method of age-depth model construction.

All sediment samples from both cores were placed into plastic cups and magnetic susceptibility was measured using a SM30 magnetic susceptibility meter (ZHinstruments Inc., Czechia) at an oscillation frequency of 8 kHz and generated magnetic field amplitude of 40 A/m. The magnetic susceptibility values [X (SI), average of three measurements] were mass normalized.

The total content of selected elements was determined by means of X-ray fluorescence (XRF) spectrometry using a handheld ED-XRF analyzer VANTA VMR with a Silicon Drift Detector (Olympus, United States). The ED-XRF analyzer was coupled with a programmable moving core holder and autonomously run by a computer (for details Kletetschka et al., 2018; Kletetschka et al., 2019b). National Institute of Standards and Technology (NIST) standard reference materials 2711a Montana II Soil and 2710a Montana I Soil were used for quality control. Finally, we focused on an interpretation of concentrations and ratios of individual elements that are likely relevant to tracking the history of the lake’s productivity as well as the erosional activity in the lake catchment. The strontium to rubidium ratio (Sr/Rb) and titanium (Ti) concentration were selected as proxies characterizing the input of material from the catchment (Cuven et al., 2010; Xu et al., 2010), phosphorus (P) and the silicon to zirconium ratio (Si/Zr) as proxies of lake productivity (Cuven et al., 2011), and the iron to manganese ratio (Fe/Mn) with iron (Fe) concentration control as proxy of near bottom redox conditions (Mackereth 1966; Davison 1993).

For separation of magnetic microspherules from the sediment samples (all SUZ1 layers), we performed a standard separation technique described by Israde-Alcántara et al. (2012). The separated material was analyzed under a dissecting optical microscope at 25–70× magnification. Suspected spheric objects (diameter >3 µm) were manually placed on aluminum stubs with nonmagnetic tools (sharpened wooden sticks). Finally, the spheric objects were identified with a scanning electron microscope (SEM) TESCAN Vega in back-scattered electron (BSE) and secondary electron (SE) modes, and subsequent elemental microanalysis was conducted using an energy dispersive X-ray spectroscope (EDS, detector X-Max 50; Oxford Instruments, United Kingdom) at the Laboratory of Scanning Electron Microscopy (Institute of Petrology and Structural Geology, Faculty of Science, Charles University, Prague, Czechia).

Diatom analysis and the analysis of freshwater fauna remains were conducted at 1 cm resolution throughout the SUZ1 sediment core. Diatom samples were prepared following the method described by van der Werff (1955). The material was cleaned by adding 37% H₂O₂ and heating to 80°C for about 1 h, then by the addition of KMnO₄. Following digestion and centrifugation, the resulting clean material was diluted with distilled water to avoid excessive concentrations of diatom valves, which may hinder reliable observations. Known quantities of Lycopodium spores were added to estimate diatom concentrations. Cleaned diatom valves were mounted in Naphrax^®, a high-refractive index medium. In each sample, 400 diatom valves were identified and enumerated on random transects at 1,000× magnification (oil-immersion), using an Olympus BX53 microscope equipped with Differential Interference Contrast (Nomarski) optics and an Olympus UC30 Imaging System. Diatoms were categorized into five ecological groups according to their life form (modified according to Buczkó et al., 2013): 1) aerophyton (living in subaerial and terrestrial habitats), 2) benthos (living at the bottom and near shore, mainly in soft sediments), 3) periphyton (attached to submerged substrata or organisms; more sessile than benthic diatoms), 4) plankton (living in the water column), and 5) tychoplankton (taxa frequently encountered both in the water column and in surface sediments). In addition, the category of benthos was further subdivided into three subcategories: epipsamon (living attached to sand grains), benthos sensu stricto (lake benthos besides epipsamon), and river benthos. Furthermore, the D to C ratio of diatom valves to chrysophyte cysts was calculated. The chrysophyte cyst sums are expressed relative to the number of diatom frustules (two valves = one frustule) plus chrysophyte cysts, with the following formula: D:C = [number of diatom frustules/(number of chrysophyte cysts + number of diatom frustules)] × 100 (Douglas and Smol 1995).

Samples for zoological indicators were washed with distilled water through a 100 µm sieve and transferred into a modified plastic Sedgewick-Rafter counting cell. All identified animal remains, mainly chironomid (non-biting midge) head capsules (HCs), were picked with either fine forceps or a steel needle using a dissecting optical microscope (15× magnification). Samples were then dehydrated in ethyl alcohol (∼96%) and mounted on glass microscope slides in Euparal® mounting medium. For taxa identification, we used an optical microscope (100–400× magnification) and followed several identification keys—mainly Wiederholm (1983), Rieradevall and Brooks (2001), Brooks et al. (2007), and Szeroczyńska and Sarmaja-Korjonen (2007). Due to low abundances of all zoological indicators, including chironomid HCs, in most of the samples, we show only concentrations of individual taxa (i.e., fossils per 1 g of dry sediment).

Statistically significant changes in the diatom record were determined using the optimum sum of squares partitioning on percentage data implementing the ZONE program (Birks and Gordon 1985; Lotter and Juggins 1991). The results from optimum sum of squares partitioning on percentage data were tested for significance using a broken-stick model to derive a stratigraphical zonation (Bennett 1996). Due to the presence of zero values in the chironomid record in some of the analyzed layers (which do not allow independent zonation), we use the diatom zonation when the zoological indicator data and other proxies are discussed.

Figure 1C displays the morphobathymetric map of Suzdalevo Lake obtained after data processing and compilation, while Supplementary Figure S2 shows echographic (depth) profiles with the depth magnified by a factor of 10. The lake is very shallow, with a maximum depth of about 2.3 m and a surface area of ∼1.8 ha. The depth-to-diameter ratio is about 0.015, and its irregular inverted-cone morphology is slightly elongated in the N-S direction. The depocenter is located close to the N shore, while the E and W shores show morphological irregularities. Both sediment cores (SUZ1 and SUZ3) were recovered from the depocenter (Figure 1C) and consisted of homogeneous dark brown gyttja with no distinct lamination or a clayey layer.

In the SUZ1 core, ²¹⁰Pb, ²²⁶Ra, and ¹³⁷Cs activities were measured only for the upper 26 cm. Due to their low values and large error bars (Supplementary Figure S3, Supplementary Table S1), the same measurements were not performed for the deeper layers. Excess (“unsupported”) ²¹⁰Pb (²¹⁰Pb_ex) was detected even in the deepest analyzed sample (depth of 25–26 cm), therefore the ²¹⁰Pb_ex inventory could not be calculated. The ¹³⁷Cs specific activities are relatively low and irregular at depths below 16 cm. The maximum ¹³⁷Cs activity of 93 Bq kg⁻¹ was detected at the depth of 9.5 cm. Weights of some SUZ1 samples were <1 g of dry sediment (average 0.9 g).

The problem of low sample weights in the SUZ1 core was resolved using the parallel core SUZ3, where four adjacent layers were merged to produce larger samples (3.1–5.6 g of dry sediment). In nine merged (4-cm-thick) layers, the specific activities of ²¹⁰Pb decrease exponentially from the topmost layers to background (“supported”) activities in the lower part of the core, showing a relatively smooth decreasing curve (Supplementary Figure S3). ²¹⁰Pb_ex was successfully detected in upper eight merged samples (Figure 2, Supplementary Table S1). The trend in the SUZ3 relation between ²¹⁰Pb_ex and depth in core is similar to the SUZ1 record, but the curve is smoother and shows lower error bars (Supplementary Figure S4). The ¹³⁷Cs specific activities are relatively low (maximum activity of 61 Bq kg⁻¹), demonstrating low contamination by radioactive fallout. An increase in ¹³⁷Cs activity, that should be equivalent with the global fallout due to nuclear weapon tests before its maximum in 1963 CE, was found between depths of 14 and 18 cm (Figure 2), with corresponding calculated ages of 1958–1963 years. CE (CRS model), 1988–1993 years. CE (CIC model), and 1963–1971 years. CE (CFCS model) for the average depth of 16 cm (Supplementary Table S2). However, this ¹³⁷Cs peak is flat and low ¹³⁷Cs activities are dispersed in all samples of the upper part of the core, probably due to continuous flushing out from the surrounding contaminated area, occasional sediment-water interface mixing (redeposition), and younger atmospheric fallout (later nuclear tests and power plant accidents).

Calculated sediment accumulation rates for the SUZ3 core differ among the age-depth models (Supplementary Table S2). Based on the CRS model, the accumulation rates decrease downwards from 3.27 (depth of 4 cm) to 1.24 (depth of 32 cm) mm year⁻¹. The uppermost samples between 0 and 3 cm show intermediate accumulation rates increasing from 2.44 to 3.10 mm year⁻¹. The CIC model shows two periods of high accumulation rate around 14 and 30 cm, and the CFCS model expects stable accumulation rate of 2.88 mm year⁻¹. The age of 1908 CE (best fit) is placed to depths of 27.75, 32.84, and 32.00 cm according to the CRS, CIC, and CFCS models, respectively.

Mass normalized magnetic susceptibility values [X (SI)] in the SUZ1 core show a peak at the depth of 27.5 cm and an increasing trend between the depths 15 and 0 cm (Figure 3). For the same core, the concentrations of selected elements (ppm) and their ratios are plotted with depth and the diatom zonation (see details in Section 3.5 and Figure 3) to better interpret the lake sedimentation history. Two proxies characterizing catchment material input to the lake are used. The Sr to Rb ratio has a similar trend as magnetic susceptibility (a peak at the depth of 26.5 cm and an increasing trend from 16 cm), except for a drop in the two uppermost layers that correspond to an increased concentration of Ti. Ti concentrations are otherwise relatively stable in the profile. Both proxies of primary productivity in the lake environment, P concentration and the Si to Zr ratio, show peaks at 37.5 cm, low values around 26.5 cm, and increased values between 24 and 12 cm (max. P concentration of ∼2,250 ppm and max. Si/Zr of ∼5,800) followed by low values (min. P concentration of ∼1,070 ppm and max. Si/Zr of ∼1850). Trends of the Fe to Mn ratio curve correspond well with Fe concentrations, showing increased values around 27.5 cm and a sharp increase with maxima in the two uppermost layers (shifts from 85 to 133 for Fe/Mn and from ∼50,500 to ∼137,200 ppm for the Fe concentration). The Fe to Mn ratio, a potential proxy for changes in the near-bottom redox conditions, was therefore determined by fluctuations in the Fe concentration.

A similar magnetic susceptibility and geochemical record was found for the SUZ3 core (Figure 3). Nevertheless, it shows two noticeable differences: 1) absence of the magnetic susceptibility peak at the depth of 27.5 cm followed by the Sr/Rb peak at the depth of 26.5 cm, and 2) absence of the magnetic susceptibility, Ti, Fe/Mn, and Fe peaks in the two uppermost layers. Based on this comparison, the SUZ1 core provided more dynamic records and was selected for extraction of iron-rich microspherules and analyses of biological indicator remains.

During the SEM observation of the magnetic extract from the Suzdalevo Lake sediment, we found reduced (oxygen content <50%) titanomagnetite (TM) grains in the depth interval of 26–30 cm (Supplementary Figure S5). However, their concentration is very low, not exceeding three grains per 1 g of dry sediment. Only one magnetic microspherule larger than 3 µm in diameter was found in the same samples with an original wet volume of 28 cm³ (14 cm³ for the two youngest layers). This iron oxide microspherule (10 µm in diameter; Fe 67.6%, O 30.4%, Cr <1%, Cu <1%) was separated from the sample 3–4 cm (∼2011 CE). Its outer surface exhibits a distinctive skeletal texture characteristic of rapid melting-quenching processes (Supplementary Figure S6).

Altogether, 156 diatom taxa representing 55 genera were identified in the sediment studied. The most abundant taxa include Staurosira costruens var. venter (Ehrenberg) P. B. Hamilton, Aulacoseira laevissima (Grunow) Krammer, Pseudostaurosira brevistriata (Grunow) D. M. Williams and Round, Aulacoseira alpigena (Grunow) Krammer, and Tabellaria flocculosa (Roth) Kützing (Supplementary Figure S7). Most of the taxa belong to benthos sensu stricto, tychoplankton, and periphyton, whereas the epipsamon, river benthos, and aerophyton ecological groups are rare (Figure 4). Although Tabellaria flocculosa is usually considered to be a planktonic diatom, periphytic forms can be also common. Therefore, we decided to keep this taxon as a separate category.

The diatom record was divided into four significant diatom assemblage zones (DAZ 1–4) based on the relative abundances. DAZ 1 was further divided into two subzones (DAZ 1a, DAZ 1b), because of striking changes in the diatom valve concentration per 1 g and diversity (number of species per 400 valves) at the depth of 28 cm (Figure 4).

The average concentration of chironomid head capsules (HCs) is 9.5 HCs per 1 g of dry sediment. While HCs in nine samples are absent, the maximum concentration of HCs in a single sample reaches 131 per 1 g of dry sediment. Altogether 33 chironomid taxa were identified (poorly preserved HCs identified to subfamily level were not counted). The most common components of the chironomid assemblages (at least 5 HCs per 1 g in one sample) are Chironomus anthracinus-type, Cladotanytarsus mancus-type, Cricotopus intersectus-type, Dicrotendipes nervosus-type, Endochironomus impar-type, Psectrocladius sordidellus-type, Tanytarsus lugens-type, and T. mendax-type (Supplementary Figure S8). In addition to chironomids (Diptera: Chironomidae), we also observed head capsules of biting midges (Diptera: Ceratopogonidae); mayfly (Ephemeroptera) mandibles; ephippia (genus Ceriodaphnia and Daphnia), shells, and head shields (genus Bosmina) of planktonic cladocerans; bryozoan statoblasts (Cristatella mucedo and genus Plumatella); and microturbelarian (Turbellaria: Rhabdocoela) resting-egg capsules (Figure 5). However, these other zoological indicators are rather complementary to the chironomid record.

The cores SUZ1 and SUZ3 were collected at the same location within the Suzdalevo Lake basin and show relatively similar geochemical, magnetic, and short-lived isotope activity records (Figure 3, Supplementary Figure S3). Due to this similarity, although partial differences are apparent, we used the ²¹⁰Pb dating of the SUZ3 core for interpretation of SUZ1 proxy data. In both cores, the observed peaks in ¹³⁷Cs (onset of the peaks at the depth of ∼16 cm; Figure 2, Supplementary Figure S4) were interpreted as signals of the global fallout from atmospheric nuclear weapons with the ¹³⁷Cs maximum in 1963 CE (e.g., Appleby 2008; Bruel and Sabatier 2020), including the tests carried out by the Soviet Union at Novaya Zemlya from 1957 to 1962 (Khalturin et al., 2004), which is in good agreement with the results of the CRS age-depth model (1958–1963 CE at the depth of 16 cm) for the SUZ3 core. This model is likely the most relevant one for dating of the younger half of the sedimentary record as it reflects well sediment compaction in middle parts of the cores. The signal of the Chernobyl accident in 1986, that was reported from many regions of the Northern Hemisphere, is not known from lake sediments investigated in the area of the TE (Rogozin et al., 2017; Darin et al., 2020). However, the maximum ¹³⁷Cs peaks in Suzdalevo Lake cores are not very prominent. In SUZ1, low ¹³⁷Cs activities were detectable also below the depth of 16 cm, likely due to partial mixing and/or remobilization and downward diffusion within the sediments (Matisoff 2017). Despite the fact that the “no sediment mixing” assumption for construction of the ²¹⁰Pb-based age-depth models is not fully met, all three models (CRS, CIC, CFCS) show that the base samples of the studied sediment cores are older than 1908 CE (Figure 2, Supplementary Table S2). Since even the oldest layers contain valves of lake diatoms and remains of lake invertebrates, sedimentation in a lake environment is apparent throughout the whole record (Figures 4, 5). Moreover, the diatom and zoological indicator assemblages below the first layer with no ²¹⁰Pb_ex activity (>36 cm, calculated ages >1908 CE) are not uniform (Figures 4, 5, Supplementary Figures S7, S8) which contradicts potential very high accumulation rate in a lake formed by the TE. Therefore, Suzdalevo Lake was evidently formed before the TE. Due to the nature of the studied material (short cores sampled by a gravity corer) it is not possible to determine the time or the process of lake formation. It can only be speculated that it may have been associated with the selective melting of permafrost that is still present in the region (Gasperini et al., 2007; Ponomarev et al., 2019). Despite the oral testimony of the Evenki people, Suzdalevo Lake is not a TE impact crater filled with water, as has been proposed for the basin of Cheko Lake (Collins et al., 2008; Rogozin et al., 2017). Note that this conclusion is supported by the depth to diameter ratio being lower than 1/100 (Figure 1C and Supplementary Figure S2), while typical simple impact craters have values between 1/5 and 1/7 (Melosh 1989; Włodarski et al., 2017).

On the other hand, the existence of the lake before the TE is not supported by identification of a specific “TE layer”. No distinct clayey TE layer was observed, in contrast to those documented in the sediments of Cheko Lake and Zapovednoe Lake (thickness from 2 to 8 mm) by Kletetschka et al. (2019a) and Darin et al. (2020), respectively (Supplementary Figure S10). This is probably due to the fact that Zapovednoe and Cheko lakes are deep riverine lakes with large catchments, while Suzdalevo is a small shallow water body separated from the Chamba River (Figure 1C, Supplementary Figures S1, S2). Sediments of riverine lakes and lakes with large tributaries are often very sensitive to erosion in a wider area compared to lakes with small catchments (e.g., Xu et al., 2010; Wilhelm et al., 2017). The most conspicuous, but not strong, clastic material input event in Suzdalevo Lake occurred in DAZ 4, based on the simultaneously increased values of mass normalized magnetic susceptibility, the Sr to Rb ratio, and Ti concentration in the SUZ1 core (Rapuc et al., 2020). However, DAZ 4 represents the youngest part of the sedimentary record and no similar erosion signal was found in the SUZ3 core. Another evidence of disturbance to the lake catchment was observed between the depths 25 and 28 cm (DAZ 1b), i.e., 14 cm above the base of the SUZ1 core (Figure 3). The DAZ 1b (SUZ1 core) contains abundant fine plant remains (Supplementary Figure S9), and its onset is characterized by increases in mass normalized magnetic susceptibility and the Sr to Rb ratio. The mass normalized magnetic susceptibility values reflect the concentration of potential magnetic carriers in the sediment (magnetite, pyrrhotite, maghemite, titanomagnetite, titanomaghemite, greigite); therefore, their enhancement in DAZ 1b reflects the presence of higher concentrations of at least one of these carriers. Since the area was exposed to fires due to the explosion and extreme radiation from the sky (Svetsov 2002; Jenniskens et al., 2019; Johnston and Stern 2019), the burning of organic material (e.g., grass, herbs, tree leaves and branches) has been associated with the conversion of iron-bearing organic matter to magnetite (Kletetschka and Banerjee 1995). The increase in the Sr to Rb ratio also reflects the higher content of plant remains in the sediment samples rather than a decline in terrigenous clastic input (Xu et al., 2010). This is supported by the relatively stable Ti concentration values, indicating no distinct change in the proportion or grain size of clastic material (Cuven et al., 2010). The above mentioned geochemical and magnetic anomalies in DAZ 1b were, however, observed only in the SUZ1 core (Figure 3). The absence of similar shifts in the core SUZ3 implies variability of the disturbance signal strength within the lake basin. The DAZ 1b onset at the depth of 28 cm is in agreement with the CRS-model-calculated depth for 1908 CE (27.75 cm), however, this “TE-depth” is placed ca. 5 and 4.3 cm deeper according to the CIC and CFCS models, respectively. It is important to notice that the time of the TE (i.e., 111 years before the sediment coring) is close to the limit of the ²¹⁰Pb dating method and the related models’ applicability. Moreover, when the CRS model is applied, a “too-old” age error is always present for the deeper core layers due to underestimation of ²¹⁰Pb_ex (Binford 1990).

We found no increase in concentration of melted magnetic microspherules in layers of potential TE age. Possible existence of a layer containing these objects would have been used for validating the age-depth models. The lack of melted magnetic microspherules in the SUZ1 sediments means that the TE explosion produced less than 354 magnetic microspherules per 1 m², assuming the diameter of our gravity corer. This agrees with the known maximum concentration of these objects near the TE epicenter (90 magnetic microspherules per 1 m², when objects with diameter from 10 to several hundred µm were counted) as reported by Badyukov et al. (2011). They studied material collected in the early 1960s by an expedition led by K. P. Florenskii, and distinguished magnetic iron-oxide microspherules and magnetic silicate microspherules with metallic droplets, both subdivided based on the presence/absence of a Ni admixture. The single microspherule we found in the SUZ1 core (Supplementary Figure S6) corresponds to the category of iron-rich spherules without a Ni admixture that was, according to Badyukov et al. (2011), likely created by a high-temperature process from terrestrial rocks (i.e., melted and evaporated due to sudden heat pulse, and later nucleated and/or quenched due to subsequent cooling). The same type of microspherules can be, however, of anthropogenic origin and transported for long distances if sufficiently small (<150 µm in diameter) (Zhang et al., 2011). The maximum calculated age of the extracted microspherule of 2003 CE (CFCS model, depth of 4 cm; Supplementary Table S2) indicates that it likely originated from the background micrometeorite flux, redeposition of an older TE sediment or soil, or anthropogenic activities. It seems that, due to their low concentration, these objects (diameter >3 µm) cannot be effectively used as markers of the TE in lake sediment cores from the Tunguska region. Also, the somewhat reduced TM grains that were found at depths of 30–26 cm (Supplementary Figure S5) cannot be interpreted as clear evidence of the TE explosion. The concentration of these grains is low (<3 grains per 1 g of dry sediment) and could have originated from the erosion of local volcanic rocks (Kamo et al., 2003). Such erosion in the lake catchment and subsequent inwash of the TM grains to the lake depocenter could have rather been associated with the TE-related forest disturbance. In future studies focused on lake sediments in the Tunguska region, the detection of increased concentrations of Ir and other Pt group elements may be more promising as follows from the “Northern” peat bog record presented by Hou et al. (2004).

In this study, we present the first diatom and freshwater zoological indicator records from the area exposed to the TE explosion in 1908 CE (Figures 4, 5). Throughout its recent history, Suzdalevo Lake lacked a significant planktonic diatom community and was dominated by tychoplanktonic Aulacoseira species together with small benthic alkaliphilic species of the Fragilariaceae family, commonly found in nutrient limited lakes with long periods of ice cover in alpine (Lotter et al., 1997; Karst-Riddoch et al., 2005), arctic (Douglas and Smol 1999; Laing and Smol 2000), and tundra regions (Westover et al., 2006). Similarly, the remains of planktonic cladocerans are rare in the sediments, demonstrating a greater suitability of the shallow lake habitats for benthic and littoral taxa (Szeroczyńska and Sarmaja-Korjonen 2007). This is supported by frequent findings of various benthic zoological indicators, such as chironomids, ceratopogonids, mayflies, bryozoans, or microturbelarians. All invertebrate taxa that were identified are known from the lakes in southern Siberia (42–65°N) or have a Holarctic distribution (e.g., Massard and Geimer 2008; Nazarova et al., 2008; Kotov 2016; Biskaborn et al., 2019).

According to the alternative age-depth models (Figure 2, Supplementary Table S2), the “TE-layer” and related signals of ecosystem changes should be located in DAZ 1a around the depth interval of 32–33 cm (CIC and CFCS models) or at the DAZ 1a/DAZ 1b transition, i.e., near the depth of 28 cm (CRS model). The first scenario is more likely in the case of the underestimation sediment accumulation rate in deeper layers mentioned by Binford (1990), however, the second scenario is also realistic when considering the error bars and thickness of the sediment layers used for the dating (4 cm). Moreover, only the second scenario corresponds to a distinct abrupt shift in lake biota (Figures 4, 5). The reality of the first scenario would imply an absence of a unique signal in the diatom and zoological indicator records. Such a conclusion would, for example, be consistent to the observations from sediments of Lake Strzeszyńskie (W Poland) that experienced the Morasko iron meteorite shower 4.9–5.3 kyr.cal. BP (Pleskot et al., 2018). The estimated entry mass of the Morasko meteoroid was between 600 and 1,100 tons and fragments of this cosmic body formed seven impact craters with diameters ranging from 20 to 90 m (Bronikowska et al., 2017).

If the age of the DAZ 1b onset is ∼1908 CE as suggested by the CRS model, the pre-TE status of Suzdalevo Lake is represented by the bottom zone (DAZ 1a) that is dominated by benthic colonial Fragilariaceae taxa (Figure 4, Supplementary Figure S7). These ecological generalists are common in lakes at high elevations and high latitudes, where ice-free intervals are short and nutrient supplies are low, although they also occur across a wide ecological gradient (Wolfe 2003). Furthermore, they inhabit the littoral zone, which is the first area to become ice-free (Westover et al., 2006). Thus, these taxa may quickly form blooms and outcompete larger diatoms during the short growing season (Smol 1988; Lotter et al., 1999; Grönlund and Kaupilla 2002). Fluctuations in abundances of all benthic zoological indicators (Figure 5 and Supplementary Figure S8) imply the alternating presence of at least one adverse factor. Such severe suppression of benthic fauna in lakes is often caused by oxygen depletion (hypoxia) or even anoxia (e.g., Moravcová et al., 2021). Also, the presence of Ceriodaphnia ephippia and Plumatella and Cristatella mucedo statoblasts indicates limiting near-bottom oxygen concentrations caused by the intensive decomposition of organic matter accompanied by pronounced thermal stratification, limiting oxygen replenishment (Ursenbacher et al., 2020). However, the stratification and oxygen depletion were not permanent, as documented by the high abundance and high taxonomic diversity of chironomid HCs at depths from 39 to 32 cm, where both taxa that are tolerant (e.g., Chironomus anthracinus-type) and sensitive (Tanytarsus lugens-type) to low oxygen concentrations (Brooks et al., 2007; Tichá et al., 2019) were observed. Short periods of more intense water column mixing (affecting the sediment-water interface) could have occurred during the recent lake history as suggested by detectable low ¹³⁷Cs activities bellow the depth of 16 cm in both cores (Figure 2, Supplementary Figure S3), but such short-term changes cannot be reconstructed due to the time resolution of our records. The hypothesized changes in the mixing regime during DAZ 1a were likely triggered by an external factor, such as a climate change or a forest disturbance. A potential impact of floods is not supported by any erosion proxies (Figure 3) or presence of taxa indicative of running water environment (Figures 4, 5). Therefore, if a flood occurred, it did not have an erosive effect on the lake shore. Another possible explanation comes from the alternative age-depth models. If the “TE-depth” calculation using CIC and CFCS models is more realistic then in the case of the CRS model, an effect of the TE explosion cannot be ruled out as well.

As mentioned above, the assemblages in DAZ 1b could represent a direct response of the lake-catchment ecosystem to the TE (Figures 4, 5). In diatoms, the occurrence of aerophyton indicates the influence of the substantial transport of terrestrial material from the lake catchment, whereas the high abundances of periphyton and benthos suggest greater substrate variability. These changes therefore correspond to the observed increase in the amount of fine plant remains and the peak in mass normalized magnetic susceptibility (SUZ1 core), documenting a disturbance of the surrounding taiga. Diatom assemblages are known to also be structured by organic matter content and/or availability (Grimes et al., 1980), which significantly increased during this period. The diversified conditions, together with the input of organic matter, most probably led to increased diatom diversity and valve concentration (Adrian et al., 1999). Our data do not imply any conditions promoting the growth of acidophilous taxa, possibly connected with acidification or acid rains, as previously suggested by Kolesnikov et al. (1998). This also applies to the diatom assemblages in DAZ 1a. However, the peak in diatom valve concentration in DAZ 1b could have been supported by mild eutrophication after the increased deposition of nitrogen (N), estimated by Kolesnikov et al. (2003) to 200,000 tons across the 2,000 km² of tree fall area. This extra N supply was interpreted as a result of the high-temperature oxidation of nitrogen in the atmosphere with formation of nitrogen oxides during the impact (Kolesnikov et al., 1998, 2003); however, it could also have arisen from another unconventional process. Boslough and Crawford (2008) modelled the TE explosion in the atmosphere and calculated an overpressure reaching 1 GPa at 6–8 km above the Earth’s surface. Such a pressure wave allows the compression of N in air to reach near the liquid phase (Alkhaldi and Kroll 2019), and thus may have allowed the oxidation and creation of N species over the impacted area. The deposition value published by Kolesnikov et al. (2003) (1 ton per 1 ha) may be overestimated or apply only to sites located close the TE epicenter (distance <10 km); however, even a smaller deposition could have affected Suzdalevo Lake, as by shown by Bergström and Jansson (2006). Moreover, an input of P-rich dust after the TE explosion-related thermal radiation and shockwave could be expected, transporting this key limiting nutrient to the studied lake and its catchment. No increase in P concentration was observed in the studied cores near the DAZ 1b onset or around the depth interval 32–32.84 cm (Figure 3), but trends in total P concentration sometimes differ from trends in biologically available P as shown by Norton et al. (2011). The taxonomic diversity and abundance of zoological indicators increase in DAZ 1b, with the dominance of Chironomus anthracinus-type and the absence of Tanytarsus lugens-type (Supplementary Figure S8) suggesting an improvement in near-bottom oxygen conditions but to a lesser extent than in the case of the period around the chironomid abundance and diversity peak at 35.5 cm (Johnson and Wiederholm 1989; Gąsiorowski and Sienkiewicz 2010). The potential destruction of the shoreline forest by the TE explosion could have caused increased wind speeds and shortened the periods of thermal stratification via wind mixing of the water column, resulting in improved dissolved oxygen concentrations (Klaus et al., 2021). On the other hand, the concurrent inwash of organic material could have accelerated respiration and oxygen consumption rates, which had the opposite effect and prevented the reappearance of the sensitive Tanytarsus lugens-type (Luoto 2013). In any case, the observed changes in DAZ 1b do not indicate any major disturbance or lake ecosystem collapse.

After the DAZ 1b, the ecosystem did not return to the pre-TE state. A new equilibrium, represented by DAZ 2, was established, and lasted until the 1990s (i.e., ca. 50–70 years) (Figure 4). We interpret this zone as a period of low water transparency and pronounced oxygen depletion, indicated by the substantial increase in the tychoplanktonic diatom taxa and extremely low abundances of all benthic zoological indicators, respectively. We assume that both processes were driven by the presence of additional dead plant biomass in the catchment after the TE. This led to a subsequent increase in dissolved organic carbon (DOC) production, and later to DOC leaching, DOC transport to the lake, and lake water brownification (Kopáček et al., 2018). The DOC leaching may have been amplified by reduced evapotranspiration in the disturbed forest, leading to increased soil moisture (Brothers et al., 2014; Kopáček et al., 2018). Brownification is a fundamental factor influencing lake ecosystem structure and function through more stable thermal stratification caused by increased absorption of solar radiation and limited light availability to benthic primary producers due to decreased water column transparency (Brothers et al., 2014; Solomon et al., 2015; Vasconcelos et al., 2016). These conditions likely caused the decline in benthic diatoms and especially in all benthic invertebrates, such as chironomids, ceratopogonids, bryozoans, and microturbelarians (Figures 4, 5). The chironomid abundances are so low in DAZ 2 (0–3 HCs per 1 g of dry sediment) that even episodic near-bottom anoxic conditions at the deepest part in the lake cannot be ruled out (Johnson and Wiederholm 1989; Houfková et al., 2017; Ursenbacher et al., 2020). These conditions would have supported Aulacoseira species (Supplementary Figure S7), which are adapted to low light availability (Kilham 1990; Bradbury et al., 1994), but most of them also indicate turbulent, unstable environments, likely caused by mixing of the water column. Such ecological characteristics are related to the production of heavy resting cells during the Aulacoseira life cycle that requires resuspension from the sediments into the lake water to establish (tycho)planktonic populations (Round et al., 1990). In addition, Aulacoseira species also create heavy silicified vegetative cells and therefore require an increased degree of turbulence to prevent their sedimentation (Saunders et al., 2009; Buczkó et al., 2013). In view of these facts, the stratification could not be long-term, but occurred frequently enough to disadvantage the benthic fauna. A similar record showing a long-term Aulacoseira presence in a shallow lake affected by anoxia was published by Houfková et al. (2017). We can speculate that a severe winter stratification due to long ice cover period was the key factor that suppressed the benthic invertebrates but enabled a summer bloom of Aulacoseira.

While the highest diatom valve concentration per 1 g of dry sediment was observed in DAZ 1b, the highest P concentrations and the highest values of the Si to Zr ratio, the other two potential proxies for algal primary production (e.g., Cuven et al., 2011), were observed in DAZ 2 (Figure 3). We explain this discrepancy by the above-mentioned presence of large and heavy silicified Aulacoseira valves in DAZ 2, which could have increased the total diatom biomass and the concentration of biogenic silica (i.e., the Si to Zr ratio) in the sediment samples. Alternatively, the increased P content in DAZ 2 could be explained by an increase in P immobilization in the sediments caused by a higher aluminum hydroxide to iron hydroxide ratio in highly organic sediments during the brownwater phase, as shown by Kopáček et al. (2007). In the latter case, P concentration values would not reflect the availability of this limiting nutrient for aquatic organisms.

Since about the 1990s (DAZ 3), the lake biota in Suzdalevo Lake began to resemble those of the DAZ 1a (Figures 4, 5). The overall dominance of benthic diatoms (around 70% of the total diatom assemblage) may be attributed to processes linked to increased water transparency and with more stable conditions after recovery of the taiga forest, when DOC concentration in the lake water decreased and mature trees on the lake shore started to attenuate water mixing by wind (Kilham, 1990; Westover et al., 2006). The higher diatom diversity could be explained by the decreased dominance of Aulacoseira taxa, leading to more even diatom communities. These new conditions were also more favorable for benthic invertebrates, including taxa associated with aquatic macrophytes, that were absent in DAZ 2.

We also attempted to document changes in oxygen availability and related redox conditions in the sedimentary record using the Fe to Mn ratio (Figure 3). This proxy is based on the fact that Mn is more easily mobilized and passes into solution more readily than Fe under reducing conditions (Davison, 1993). However, in Suzdalevo Lake, trends in the Fe to Mn ratio coincide with the Fe concentrations, indicating that the ratio was mainly affected by changes in Fe supply from the lake catchment (Mackereth 1966), and thus could not be used to reconstruct changes in redox conditions.

The SUZ1 core provided an anomalous multi-proxy record for the uppermost two layers (DAZ 4), representing approximately the years 2016–2019 CE. In spite of the absence of these anomalies in the SUZ3 core that indicates variability of the signal strength within the lake basin or partial loss of the surface layer during the coring, the shifts observed in most of the proxies in the SUZ1 core are conspicuous and deserve a separate comment (Figure 3). The low Sr to Rb ratio and the sharp increase in Ti concentrations in DAZ 4 document an increased input of clastic material (Cuven et al., 2010; Xu et al., 2010), which was not observed in the other zones, including the depth interval where “TE-depth” is expected according to the age-depth models. Such a change in the proportion of inorganic particles could have been caused by increased erosion in the catchment or by flood activity of the nearby Chamba River (Supplementary Figure S1). Also, river benthic diatom taxa, such as P. frequentissimum (Figure 4, Supplementary Figure S7), reach their maximum abundance in DAZ 4, indicating a possible riverine influx and, together with higher abundances of aerophyton, may reflect the increased input of allochthonous material. In addition, the low diatom diversity and valve concentration and the dominance of S. pinnata and T. flocculosa suggest the persistence of relatively unproductive conditions throughout this period (Dam et al., 1994; Hofmann et al., 2011).

The same zone is characterized by the presence of abundant remains of two planktonic cladocerans, Daphnia and Bosmina, while the abundances of the benthic zoological indicators remain rather low (Figure 5). Among these benthic invertebrates, the only exception is the slight increase in chironomids associated with aquatic macrophytes and the solitary finding of the rheophilic chironomid Eukiefferiella claripennis-type (Brooks et al., 2007) (Supplementary Figure S8). Thus, the zoological record also suggests the influence of an external factor that changed conditions in the lake. DOC concentrations in boreal streams are often the lowest after high discharge events arising from snowmelt (e.g., Jonsson et al., 2007). If the factor responsible for the observed changes in lake biota was a snowmelt flood (or several floods during the period 2016–2019 CE) on the Chamba River, a drop in DOC accompanied with an increase in water transparency, along with other changes in water chemistry, could have occurred in Suzdalevo Lake. Such a shift would, for example, support the growth of aquatic plants and increase the abundances of the associated chironomid taxa (Figure 5). Due to the presence of Bosmina, a cladoceran well adapted to fish predation, short-time colonization of the lake by fish after a flood event cannot be ruled out either (Johnsen and Raddum 1987); however, we did not observe any fish during our field survey.

We do not have data on floods on the Chamba River and we do not know of any observations of changes in the connectivity between Suzdalevo Lake and the river, as the site is very remote. Nevertheless, based on available data from the nearest meteorological station in Vanavara (Weather and Climate 2021), severe floods caused by rapid snowmelt or an ice jam near the study site could have occurred during the given period. At this station, the average annual temperature and precipitation between the years 1933 and 2020 were -5.5°C and 421 mm, respectively, while between the years 2016 and 2019 the average annual temperature increased to −4.1°C and average annual precipitation increased to 625 mm. This change in the local climate would have been favorable for stronger snowmelt floods and corresponds to trends observed in the wider area of Central Siberia (e.g., Tchebakova et al., 2011).