Spatial distribution of environmental indicators in surface sediments of Lake Bolshoe Toko, Yakutia, Russia

Rapidly changing climate in the Northern Hemisphere and associated socio-economic impacts require reliable understanding of lake systems as important freshwater resources and sensitive sentinels of environmental change. To better understand time-series data in lake sediment cores, it is necessary to gain information on within-lake spatial variabilities of environmental indicator data. Therefore, we retrieved a set of 38 samples from the sediment surface along spatial habitat gradients in the boreal, deep, and yet pristine Lake Bolshoe Toko in southern Yakutia, Russia. Our methods comprise laboratory analyses of the sediments for multiple proxy parameters, including diatom and chironomid taxonomy, oxygen isotopes from diatom silica, grain-size distributions, elemental compositions (XRF), organic carbon content, and mineralogy (XRD). We analysed the lake water for cations, anions, and isotopes. Our results show that the diatom assemblages are strongly influenced by water depth and dominated by planktonic species, i.e. Pliocaenicus bolshetokoensis. Species richness and diversity are higher in the northern part of the lake basin, associated with the availability of benthic, i.e. periphytic, niches in shallower waters. δ¹⁸O_diatom values are higher in the deeper south-western part of the lake, probably related to water temperature differences. The highest amount of the chironomid taxa underrepresented in the training set used for palaeoclimate inference was found close to the Utuk River and at southern littoral and profundal sites. Abiotic sediment components are not symmetrically distributed in the lake basin, but vary along restricted areas of differential environmental forcing. Grain size and organic matter are mainly controlled by both river input and water depth. Mineral (XRD) data distributions are influenced by the methamorphic lithology of the Stanovoy mountain range, while elements (XRF) are intermingled due to catchment and diagenetic differences. We conclude that the lake represents a valuable archive for multiproxy environmental reconstruction based on diatoms (including oxygen isotopes), chironomids, and sediment–geochemical parameters. Our analyses suggest multiple coring locations preferably at intermediate depth in the northern basin and the deep part in the central basin, to account for representative bioindicator distributions and higher temporal resolution, respectively.

1 Introduction

Over the past few decades, the atmosphere in boreal and high-elevation regions has warmed faster than anywhere else on Earth (Mountain Research Initiative EDW Working Group, 2015; Huang et al., 2017). Dramatic socio-economic and ecological consequences are expected (AMAP, 2017) as well as substantial feedbacks from thawing permafrost and the associated release of greenhouse gas into the global climate system (Schuur et al., 2015). Boreal Russia is identified as a global hotspot where surface air temperature increases have led to substantial ground warming over the past decade (Biskaborn et al., 2019a). Accurate estimates of the amplitude of environmental impacts are compounded by an imprecise understanding of ecological indicators of past environmental conditions (Miller et al., 2010). Lake ecosystems, whose development is archived in their sediments, act as sensitive sentinels of environmental changes (Adrian et al., 2009), while even small changes in climate can profoundly deteriorate ecosystem services (Saulnier-Talbot et al., 2014). Assessments of the impact of climate change to lake systems rely on careful interpretation of suitable proxy data. Proxy information on present and past ecological conditions is provided by various biological and physicochemical properties of the sediment components (Meyer et al., 2015; Solovieva et al., 2015; Nazarova et al., 2017a). However, the spatial within-lake distributions of preserved remnants of ecosystem inhabitants and associated sediment–geochemical properties depend on habitat differences between the epilimnion and the hypolimnion (Raposeiro et al., 2018) and are therefore expected to be non-uniform. Accordingly, precise palaeolimnological reconstruction of past environmental variability requires a detailed, quantitative understanding of the modern (21st century) within-lake heterogeneity.

Here, we employ a multi-proxy approach to attain a holistic view of a lake's depositional history in boreal Russia. Variables include diatom and chironomid taxonomy, stable oxygen isotopes in diatom silica, grain-size distributions, elemental compositions, organic carbon content, and mineralogy. Abiotic sediment properties may represent signals resulting from either the external input of material and lake-internal conditions during deposition or post-sedimentary diagenetic processes near the sediment surface (Biskaborn et al., 2013b; Bouchard et al., 2016). Hence, our integrated approach enables the identification and distinction between internal lake processes and external forcing (Cohen, 2003).

Diatoms (unicellular, siliceous microalgae) are major aquatic primary producers. They appear ubiquitously and their opaline frustules (SiO₂•nH₂O) are well preserved in the sedimentary record, allowing exact identification down to sub-species level by high-resolution light microscope analysis (Battarbee et al., 2001). Diatoms are widely applied bioindicators for past and present ecosystem changes in boreal environments (Miller et al., 2010; Pestryakova et al., 2012; Hoff et al., 2015; Herzschuh et al., 2013; Biskaborn et al., 2012, 2016; Palagushkina et al., 2017; Douglas and Smol, 2010). Widespread responses of planktonic diatoms to recent climate change indicate that lakes in the Northern Hemisphere have already crossed important ecological thresholds (Smol and Douglas, 2007; Rühland et al., 2008). The very rapid cell lifecycles of days to weeks (Round et al., 1990) enable changes in diatom assemblages on very short timescales in response to changes in environmental circumstances, e.g. cooling or warming (Anderson, 1990). The link between climate change and diatoms, however, cannot easily be addressed via simple temperature-inference models and instead requires a more complete understanding of the interactions between the aquatic ecosystem with lake habitat preferences, hydrodynamics, and catchment properties (Anderson, 2000; Palagushkina et al., 2012; Biskaborn et al., 2016; Bracht-Flyr and Fritz, 2012; Hoff et al., 2015). It is thus necessary to identify the relationship between diatom species occurrence, the isotopic composition of their opaline valves, and internal physico-limnological factors (Heinecke et al., 2017) within spatial heterogenic lake systems before drawing direct inferences about external climatic-driven factors from single core studies.

Chironomid larvae (Insecta: Diptera) comprise up to 90 % of the aquatic secondary production (Herren et al., 2017; Nazarova et al., 2004) and hence their preserved head capsules well represent the aquatic heterotrophic bottom-dwelling ecosystem component (Nazarova et al., 2008; Syrykh et al., 2017; Brooks et al., 2007). Furthermore, the literature reports a net mutualism of chironomids and benthic algae between the primary consumer and primary producer trophic levels in benthic ecosystems (Specziar et al., 2018; Zinchenko et al., 2014). Factors influencing the spatial distribution of chironomids within single lakes are water temperature (Nazarova et al., 2011; Luoto and Ojala, 2018), sedimentological habitat characteristics (Heling et al., 2018) and/or water depth and nutrients (Yang et al., 2017), as well as hypolimnetic oxygen (Stief et al., 2005) and the availability of water plants (Raposeiro et al., 2018; Wang et al., 2012b).

As previous studies described, pollen distributions in lake sediments are less influenced by lake zonation than aquatic communities (Zhao et al., 2006). Accordingly, our study does not consider spatial pollen distributions.

Secondary factors influencing the spatial distribution of subfossil assemblages are selective transitions from living communities to accumulation of dead remains. Both biological remains and physico-chemical properties are influenced by sediment resuspension and redistribution processes described as sediment focusing (Hilton et al., 1986). These are primarily dependent on slope steepness (Hakanson, 1977) or, in shallow areas, wind-induced bottom shear stress (Bennion et al., 2010; Yang et al., 2009). Nevertheless, it already has been proven for other lake sites that within-lake bioindicator distributions are laterally non-uniform, contradicting the assumption that mixing processes cause homogenous microfossil assemblages before deposition (Anderson, 1990; Wolfe, 1996; Anderson et al., 1994; Earle et al., 1988; Kingston et al., 1983; Puusepp and Punning, 2011; Stewart and Lamoureux, 2012; Yang et al., 2009). However, many palaeolimnological studies employ single-site approaches using only one sediment core and hence do not encompass the full spatial extent and natural variability of the entire lake sediment archive. Heggen et al. (2012) report that sediment cores from the deep centre of small and shallow lakes with high spatial proxy variability in the littoral zones contain representative bioindicator assemblages. The authors also conclude that in larger and deeper lakes similar multi-site studies are necessary to make recommendations about the “ideal” coring positions for multi-proxy palaeolimnological studies.

In this respect, our broad research question is how spatially reliable palaeolimnological proxy data in a complex lake system are. To answer this question, we set up our research hypothesis: bioindicators and abiotic sediment properties will respond to different habitat conditions and lake zonation, including water depth, proximity to the main inflow in the south, and old moraines in the north of Lake Bolshoe Toko.

An analysis of spatio-temporal within-lake bioindicator distribution requires a suitable and large lake system with an anthropogenically untouched ecosystem and sufficient variability in water depth, catchment setting, and sedimentological regime. These demands are met by Lake Bolshoe Toko, the deepest lake in Yakutia, Russia (Zhirkov et al., 2016) (Fig. 1). Our study aims to gain a better local understanding of proxy data for future palaeoenvironmental analyses of long sediment cores from Bolshoe Toko. Therefore, our objectives are to (1) detect the spatial variability of abiotic (elements, minerals, grain size) and biotic (diatoms, chironomids, organic carbon) components of the lake's surface sediments, (2) reveal the causal relationship between the distribution of aquatic microfossils, lake basin features, and sedimentary parameters, and (3) attribute proxy variability to specific environmental factors.

2 Study site

Lake Bolshoe Toko (56^∘15^′ N, 130^∘30^′ E; 903 m.a.s.l.) is an oligotrophic, freshwater lake located in the Sakha Republic, Russia (Fig. 1). The lake surface area is 84.3 km², with a mean water depth of 29.5 m (maximum, 72.5 m) and a Secchi depth of 9.8 m (Zhirkov et al., 2016). The Utuk River runs through Lake Maloe Toko and brings water from the southern igneous catchment. Lake Maloe Toko (called “small Toko”, size 2.7 × 0.9 km, 168 m depth, tectonic origin) is located between high mountains south of Bolshoe Toko. The river inflow south of Bolshoe Toko forms deltaic sediments. The bay in the south-east is called Zaliv Rybachiy (“Fishing bay”). It is partly separated from the main basin and supplied with water by a small creek that itself is connected to a small lake (Fig. 1). The bay is reported to have somewhat different fauna as compared to the Bolshoe Toko main basin, i.e. occurrence of fish that are typical for small lakes and not found out of the basin (Semenov, 2018). The “Banya lake” in the north-east is isolated from Bolshoe Toko and is not considered in this study. The Mulam River is the lake's predominant outflow towards the north along the south-eastern border of Yakutia flowing into the Uchur, Aldan, and finally Lena rivers.

Figure 1 Lake Bolshoe Toko study site. (a) Geological map, bathymetry, and moraines. Map compiled using data from Konstantinov (2000) and Kornilov (1962). (b) Overview map of Siberia. World Borders data are derived from http://thematicmapping.org/downloads/world_borders.php (last access: April 2016) and licensed under CC BY-SA 3.0. (c) Catchment area around Bolshoe Toko delineated from the ASTER GDEM V2 model between the latitudes 54 and 56^∘ N and longitudes 130 to 131^∘ E (1) (Tachikawa et al., 2011) and a corresponding multispectral Landsat 8 OLI TIRS satellite image using QGIS (QGIS-Team, 2016). Most of the river catchment is located in the igneous Precambrian Stanovoy mountain range supplying the southern part of the lake with water and sediment. The shallower north-eastern part of the lake is influenced by the surrounding moraines and Mesozoic sandstones and mudstones.

There are no permanent settlements in the study area. During the time of field work there was a temporary mining settlement (built in 2011) located 17 km north-west of Bolshoe Toko in the upper course of the Elga River. This settlement was accessible by off-road vehicles we used to reach the lake, partly along temporary winter roads (frozen rivers and lakes) in March 2013. The exploitation of the El'ginsky coal deposits, planned for a productivity of 15–20 million t yr⁻¹ (Konstantinov, 2000), will strongly affect the lake and its catchment. The territory of the watershed will increasingly be damaged and contaminated by off-road vehicles, and rainfall will produce muddy water which potentially can cause lake pollution (Sobakina and Solomonov, 2013).

The lake basin is adjoined to the northern slope of the eastern Stanovoy mountain range in a depression of tectonic and glacial origin between two north-west-trending right-lateral strike-slip faults (Imaeva et al., 2009). A southward thrust fault runs along the southern border of the lake separating the Precambrian igneous rocks in the south from sandstones and mudstones of the Mesozoic Tokinski Plateau in the north. The Stanovoy mountain range in the southern catchment of the lake consists mainly of highly mafic granulites and other high-pressure metamorphic rock types (Rundqvist and Mitrofanov, 1993). At its north-eastern margins the lake is bordered by moraines of three different glacial sub-periods (Kornilov, 1962) (Fig. 1).

The study area is situated within the East Siberian continental temperate climate zone exhibiting taiga vegetation (boreal forests) and fragments of steppes and a predominant westerly wind system (Shahgedanova, 2002). The meteorological station in Yakutsk has recorded historical climate data (Gavrilova, 1993). In the 19th century the mean annual temperature (January–December) was circa −11 to −11.5 ^∘C and during the 20th century these temperatures increased to around −10.2 ^∘C, in parallel with an increase in precipitation from 205 to 250 mm yr⁻¹ (Konstantinov, 2000). The meteorological station “Toko” located approximately 10 km north-east of the lake, however, recorded an increase in air temperature of ca. 0.48 ^∘C per decade from the 1970s to 2010 (calculated from NOAA data, only those years involved that have average air T data in 12 months). Measurements taken directly at the lake were lower, indicating the influence of cold meltwater from the Stanovoy mountain range in summer and the high volume of ice during wintertime. Since the average air temperature in southern Yakutia increases with height (temperature inversion of ∼2 ^∘C 100 m⁻¹), permafrost can be locally discontinuous where taliks (unfrozen zones) underneath topographically high and deep lakes penetrate the permafrost zone (Konstantinov, 1986). As observed in 1971 (Konstantinov, 2000), ice cover lasts at least partly until mid-July.

3 Materials and methods

3.1 Field work

Field work was conducted during the “Yakutia 2013” German–Russian expedition between 19 March and 14 April 2013 by the Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research (AWI) and North-Eastern Federal State University in Yakutsk (NEFU). Vertical holes were drilled in the lake ice cover using a Jiffy ice auger with a diameter of 250 mm. Lake basin bathymetry was measured using a portable echo sounder. Ice cores were retrieved by drilling multiple holes around a central part. Water samples for hydrochemical analyses were collected prior to sediment coring using a UWITEC water sampler. Water samples were analysed in situ using a WTW Multilab 340i for pH, conductivity, and oxygen values at the day of retrieval during field work. A sub-sample of the original water was passed through a 0.45 µm cellulose–acetate filter, stored, and transported in 60 mL Nalgene polyethylene bottles for subsequent anion and cation analyses in AWI laboratories in autumn 2013. Cation samples were acidified during field work with HNO₃, suprapure (65 %), to prevent microbial conversion processes and adsorptive accretion.

At 42 sites within the lake, short cores containing intact sediment surface material were retrieved using an UWITEC gravity corer. Water depth at sampling sites was measured using either a hand-held HONDEX PS-7 LCD digital sounder and/or the cord of the coring device when the lake ice cover disturbed the signal. The sediment was identified as clayish silt deposits with a predominantly dark (black) colour and a weak smell of hydrogen sulfide, a sticky and viscous mud mixed with plant and other organic residues. The uppermost ca. 2 cm at some sites had a dark red colouring indicating the redox boundary between oxygenated and anoxic sediments. We identified the uppermost 0.5 cm of short cores as surface sediments and subsampled these layers on site during fieldwork to avoid sediment mixture during transport. Sediment samples were transported in sterile Whirl-Pak bags and sediment cores were transported in plastic liners to the AWI laboratories in Potsdam, Germany, and stored at 4 ^∘C in a dark room for further analyses and as back-up.

During this expedition long core material was also retrieved from multiple sites, including the northern and central parts of the lake, and is planned for publication in a separate paper.

According to the amount of the uppermost 0–0.5 cm layer in the short cores available, the sample size n for different sediment properties measured in this study varies.

3.2 Laboratory analyses

3.2.1 Hydrochemistry

Water depth profiles were taken during the March 2013 expedition from the deepest part of the lake (PG2208, water depth 70 m) and in the lagoon (PG2122, 18 m) as well as in August 2012 (sample site near the western shoreline, 37 m). The temperature was determined in the field and the samples analysed for isotopes (δ¹⁸O, δD; see Fig. 6). From the water samples anions were analysed using ion chromatography (Dionex DX 320) and cations were determined using inductively coupled plasma-optical emission spectrometry (ICP-OES, Perkin-Elmer Optima 8300DV Perkin-Elmer – Optical Emission Spectrometer). Alkalinity was measured by titration with 0.01 M HCl using an automatic titrator (Metrohm 794 Basic Titrino).

Stable hydrogen and oxygen isotope analyses were carried out with Finnigan MAT Delta-S mass spectrometers with two equilibration units using common equilibration techniques (Meyer et al., 2000) and given as δ¹⁸O and δD in ‰ vs. VSMOW (Vienna Standard Mean Ocean Water) with respective analytical errors of better than ±0.1 % and 0.8 ‰. The secondary parameter d excess (d) is calculated as d = δD-8δ¹⁸O (Dansgaard, 1964; Merlivat and Jouzel, 1979).

3.2.2 X-ray fluorescence and X-ray diffractometry

To gain information on the variability of the elemental sediment composition, 20 freeze-dried and milled surface samples were semi-quantitatively analysed by X-ray fluorescence (XRF) using a novel single-sample modification for the AVAATECH XRF core scanner at AWI Bremerhaven. A rhodium X-ray tube was warmed up to 1.75 and 3 mA with a detector count time of 10 and 15 s for elemental analysis at 10 kV (no filter) and 30 kV (Pd thin filter), respectively. The average modelled chi square values (χ²) of the measured peak intensity curve fitting for the relevant elements were variable but generally low (Zr = 0.92, Mn = 1.49, Fe = 2.32, Ti = 1.53, Br = 3.65, Sr = 4.79, Rb = 4.98, Si = 16.11). Values above 3 were ascribed to suspiciously high count rates from sample PG2133, which was subsequently excluded from XRF interpretation. The relatively low amount of total sample material available did not facilitate the removal of organic matter prior to sample measurement and may have contributed to the variable modelled chi square values.

As interpretation of raw device obtained element intensities (in counts per second, cps) is problematic due to non-linear matrix effects and variations in sample density, water content, and grain size (Tjallingii et al., 2007), cps values were transformed using a centred log-ratio transformation (CLR). Element ratios were calculated from raw cps values and transformed using an additive log-ratio transformation (ALR) (Weltje and Tjallingii, 2008).

The mineralogical composition of 32 freeze-dried and milled samples was analysed by standard X-ray diffractometry (XRD) using a Philips PW1820 goniometer at AWI Bremerhaven applying cobalt–potassium alpha (CoKα) radiation (40 kV, 40 mA) as outlined in Petschick et al. (1996). The intensity of diffracted radiation was calculated as counts of peak areas using XRD processing software MacDiff 4.0.7 (freeware developed by Rainer Petschick in 1999). Individual mineral content was expressed as percentages of bulk sediment XRD counts (Voigt, 2009). Mineral inspection focused on quartz, plagioclase and K-feldspar, hornblende, mica, and pyrite. Clay minerals involved kaolinite, smectite and chlorite. Accuracy of the semi-quantitative XRD method is estimated to be between 5 % and 10 % (Gingele et al., 2001).

3.2.3 Grain-size, carbon, and nitrogen analyses

In order to gain high-resolution information on the spatial variability of particle sizes and related water energy in the lake, we analysed the grain-size distribution using a laser technique. Organic material was removed from 32 surface sediment samples by hydrogen peroxide oxidation over 4 weeks on a platform shaker. Two homogenized subsamples were weighted and 93 subclasses between 0.375 and 2000 µm were measured using a Coulter LS 200 Laser Diffraction Particle Analyser. Grain-size fractions coarser than 2 mm were sieved out, weighted, and added to the volume percentage data afterwards to indicate the proportion of gravel.

To assess the accumulation of organic matter in the lake, we analysed total carbon (TC) and total nitrogen (TN) of 35 freeze-dried and milled samples. For TC and TN we quantified bulk samples by heating the material in small tin capsules using a Vario EL III CNS analyser. Total organic carbon (TOC) was measured using a Vario MAX C in per cent by weight (wt %). The measurement accuracy was 0.1 wt % for TOC and TN, and 0.05 wt % for TC. TOC and TN were compared to calculate the TOC∕TN_atomic ratio by multiplying with the ratio of atomic weights of nitrogen and carbon following Meyers and Teranes (2002).

To gain additional bioproductivity information we analysed the stable carbon isotope composition δ¹³C of the total organic carbon fraction in 15 samples using a Finnigan Delta-S mass spectrometer. Dried, milled, and carbonate-free (HCl treated) samples were combusted in tin capsules to CO₂. Results are expressed as δ¹³C values relative to the PDB standard in parts per thousand (‰) with an error of ±0.15 %.

Radiocarbon dating of two bulk sediment surface sample from short cores, each ranging from 0 to 0.5 cm depth below the sediment surface, was performed in the Poznan Radiocarbon Laboratory on the soluble (SOL) fraction using an accelerator mass spectrometer.

3.2.4 Diatoms

Twenty-three samples were prepared for diatom analysis following the standard procedure (Battarbee et al., 2001). To calculate the diatom valve concentration (DVC), 5×10⁶ microspheres were added to each sample following organic removal with hydrogen peroxide. Diatom slides were prepared on a hot plate using Naphrax mounting medium. For the identification of diatoms to the lowest possible taxonomic level we used several diatom flora including Lange-Bertalot et al. (2011), Lange-Bertalot and Metzeltin (1996), Krammer and Lange-Bertalot (1986–1991), and Lange-Bertalot and Genkal (1999). For rare taxa (i.e. Pliocaenicus) literature research was applied in scientific papers, including Cremer and Van de Vijver (2006) and Genkal et al. (2018). A minimum of 300 (and up to 400) diatom valves were counted in each sample using a Zeiss AXIO Scope.A1 light microscope with a Plan-Apochromat 100×/1.4 Oil Ph3 objective at 1000× magnification. Identification of small diatom species was verified using a scanning electron microscope (SEM) at the GeoForschungsZentrum Potsdam.

During counting of diatom valves, chrysophycean stomatocysts and Mallomonas were counted but not further taxonomically identified. Count numbers were used to estimate the chrysophyte cyst to diatom index (C:D) and Mallomonas to diatom index (M:D) relative to counted diatom cells (Smol, 1984; Smol and Boucherle, 1985). Diatom valve preservation was measured and calculated as the f-index (Ryves et al., 2001). Diatom valve concentration was estimated as the number of valves per gram dry sediment following Battarbee and Kneen (1982).

3.2.5 Oxygen isotopes of diatom silica

To analyse the oxygen isotope composition from diatom silica (δ¹⁸O_diatom) from nine representative surface samples, a purification procedure including wet chemistry (to remove organic matter and carbonates) and heavy liquid separation was applied for the fraction <10 µm following the method described in Chapligin et al. (2012a). After freeze-drying the samples were treated with H₂O₂ (32 %) and HCl (10 %) to remove organic matter and carbonates and wet sieved into <10 and >10 µm fractions. Four multiple heavy liquid separation (HLS) steps with varying densities (from 2.25 to 2.15 g cm⁻³) were then applied using a sodium polytungstate (SPT) solution before being exposed to a mixture of HClO₄ (65 %) and HNO₃ (65 %) to remove any remaining micro-organics.

To remove exchangeable hydrous groups from the diatom valve structure (amorphous silica SiO₂⋅nH₂O), inert gas flow dehydration was performed (Chapligin et al., 2010). Oxygen isotope analyses were performed on dehydrated samples using a laser fluorination technique (with BrF₅ as a reagent to liberate O₂) and then directly measured against an oxygen reference of a known isotopic composition using a PDZ Europa 2020 mass spectrometer (MS2020, now supplied by Sercon Ltd., UK). The long-term analytical reproducibility (1σ) is ±0.25 ‰ (Chapligin et al., 2010).

Every fifth sample was a biogenic working standard to verify the quality of the analyses. For this, the BFC biogenic working standard calibrated within an inter-laboratory comparison was used (Chapligin, 2011). With a δ¹⁸O value of $+\mathrm{29.0}\pm \mathrm{0.3}$ ‰ (1σ), BFC (this study: $+\mathrm{28.7}\pm \mathrm{0.17}$ ‰, n=49) is the closest diatom working standard to the Bolshoe Toko samples (δ¹⁸O values range between +22 and +24 ‰) available. A contamination correction was applied to δ¹⁸O_diatom using a geochemical mass-balance approach (Chapligin et al., 2012a; Swann et al., 2007) determining the contamination end-member by analysing the heavy fractions after the first heavy liquid separation resulting in ${\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}=\mathrm{16.2}\pm \mathrm{1.3}$ % (via EDX; n=9) and δ¹⁸O = 8.5±0.8 ‰ (n=6).

3.2.6 Chironomids

Treatment of 18 sediment samples for chironomid analysis followed standard techniques described in Brooks et al. (2007). Subsamples of wet sediments were deflocculated in 10 % KOH, heated to 70 ^∘C for up to 10 min, to which boiling water was added, and left to stand for up to another 20 min. The sediment was passed through stacked 225 and 90 µm sieves. Chironomid larval head capsules were picked out of a grooved Bogorov sorting tray under a stereomicroscope at 25–40× magnifications and were mounted in Hydromatrix two at a time, ventral side up, under a 6 mm diameter cover slip; 48 to 117 chironomid larval head capsules were extracted from each sample to capture the maximum diversity of the chironomid population. Chironomids were identified to the highest taxonomic resolution possible with reference to Wiederholm (1983) and Brooks et al. (2007). Information on the ecology of chironomid taxa and groups was taken from Brooks et al. (2007), Moller Pillot (2009), and Nazarova et al. (2011, 2015, 2008, 2017b). Ecological information of the taxa associated with biotopes (littoral, profundal), water velocity (standing, running water), and relation to presence of macrophytes were taken from Brooks et al. (2007) and Moller Pillot (2009). T July optima of chironomids were taken from Far East (FE) chironomid-based temperature inference model (Nazarova et al., 2015). The FE chironomid-based temperature inference model (WA-PLS, two components; r² boot = 0.81; RMSEP boot = 1.43 ^∘C) was established from a modern calibration data set of 88 lakes and 135 taxa from the Russian Far East (53–75^∘ N, 141–163^∘ E; T July range 1.8–13.3 ^∘C). Mean July air temperature for the lakes from the calibration data set was derived from New et al. (2002). All modern and chironomid-inferred temperatures were corrected to 0 m a.s.l. using a modern July air temperature lapse rate of 6 ^∘C km⁻¹ (Livingstone et al., 1999; Heiri et al., 2014).

3.3 Statistical analyses

Detrended correspondence analysis (DCA) with detrending by segments was performed on the chironomid and diatom data (rare taxa down-weighted) to determine the lengths of the sampled environmental gradients, from which we decided whether unimodal or linear statistical techniques would be the most appropriate for the data analysis (Birks, 1995). For diatom data the gradient lengths of the species scores were 2.07 and 1.49 standard deviation units (SDU) for DCA 1 and 2, respectively, suggesting that lineal numerical methods should be used. A principal component analysis (PCA) was used to explore the main taxonomic variation of the data (ter Braak and Prentice, 1988). The gradient lengths of chironomid species scores were 3.78 and 4.12 SDU, indicating that numerical methods based on a unimodal response model should be more appropriate to assess the variation structure of the chironomid assemblages (ter Braak, 1995). However, test PCA performed on chironomid data showed that the lineal method captures more variance of species data (ESM, Table 2); therefore, we further applied lineal methods for both chironomid and diatom data. In order to summarize the response of lacustrine biota to abiotic, physicochemical explanatory variables, a redundancy analysis (RDA) was performed on diatom and chironomid data in comparison to environmental variables (Figs. 2 and 3).

Figure 2 RDA biplots of diatoms in the surface sediments of Lake Bolshoe Toko. (a) Common diatom taxa and significant environmental variables. (b) Diatom sampling sites and significant environmental variables. Abbreviated species names: AchBio – Psammothidium bioretti; AchDid – Achnanthes cf. didyma; AchHel – Psammothidium helveticum; AchLae – Eucocconeis laevis; AchLan – Planothidium lanceolata; AchLat – Karayevia laterostrata; AchLev – Psammothidium levanderi; AchOes – Planothidium oestrupii; AchSp – Achnanthes sp.; AchSub – Psammothidium subatomoides; AchSuc – Karayevia suchlandtii; AmpPed – Amphora pediculus; AulAmb – Aulacoseira ambigua; AulSub – Aulacoseira subarctica; AulVal – Aulacoseira valida; CocPla – Cocconeis placentula; CycAto_ – Cyclotella cf. atomus; CycCom – Cyclotella comensis; CycCyc_ – Cyclotella cyclopuncta; CycIri – Cyclotella iris; CycSp – Cyclotella comensis-tripartita-complex; CymMin – Encyonema minutum; CymSin – Cymbella sinuata; DiaMes_ – Diatoma mesodon; EncLun_ – Encyonema lunatum; FraCap – Fragilaria capucina; FraConVe – Staurosira venter; FraPin – Fragilaria pinnata; GomAcu – Gomphonema acuminatum; GomIns – Gomphonema insigne; HipCos – Hippodonta costulata; NavCoc – Cavinula cocconeiformis; NavJae – Cavinula jaernefeltii; NavPsc – Cavinula pseudoscutiformis; NitSp – Nitzschia sp.; PlioBol – Pliocaenicus bolshetokoensis; StaCon_ – Staurosira construens.

Figure 3 RDA biplots of chironomids in the surface sediments of Lake Bolshoe Toko. (a) Common chironomid taxa and significant environmental variables. (b) Chironomid sampling sites and significant environmental variables. Abbreviated species names: ChaetPi – Chaetocladius piger-type; ChirLar – Chironomini larvula; CladMan – Cladotanytarsus mancus-type; Cory arc – Corynoneura arctica-type; CoryCor – Corynoneura coronata-type; Cricbi – Cricotopus bicinctus-type; CricInt – Cricotopus intersectus-type; CricoP – Cricotopus type P; Cryptoch – Cryptochironomus; Diamesi – Diamesa sp.; Eukiefe – Eukiefferiella; EukiCla – Eukiefferiella claripennis-type; HeteGri – Heterotrissocladius grimshawi-type; HetMaeo – Heterotrissocladius maeaeri-type 1; HetMaet – Heterotrissocladius maeaeri-type 2; HeteMar – Heterotrissocladius marcidus-type; Limnophy – Limnophyes – Paralimnophyes; MicroIn – Micropsectra insignilobus-type; MicrPed – Microtendipes pedellus-type; Orthocla – Orthocladius/Cricotopus; OrthOli – Orthocladius oliveri-type; OrthoS – Orthocladius type S; Paraclop – Paracladopelma; Paracri – Paracricotopus; ParaBat – Parakiefferiella bathophila-type; ParaTri – Parakiefferiella triquetra-type; Procladi – Procladius; Prodiam – Prodiamesa; Propsil – Propsilocerus type N; Protanyp – Protanypus; Psectro – Psectrocladius narrow; Pseudoch – Pseudochironomus; SergCor – Sergentia coracina-type; Smittia – Smittia – Parasmittia; Zavrelia – Stempellinella – Zavrelia; Synortho – Synorthocladius; TanyMen – Tanytarsus mendax-type; Tanytar – Tanytarsus pallidicornis-type 2; Tveteni – Tvetenia bavarica-type; ZaluMuc – Zalutschia mucronata-type; Chirono – Chironomini unidentified; Unid Tan – Tanytarsini unidentified; Unid Pen – Tanypodinae  unidentified.

Initially, all environmental variables shown in this paper were tested in a RDA to assess the relationships between the distribution of bioindicator taxa and abiotic habitat parameters. Apart from the chemical and physical parameters of the lake and sediments (Fig. 5), we include in the analysis the presence/absence of the submerged vegetation and distances of the sampling stations from the shore and from the inflowing rivers. All explanatory variables were tested for normality prior to the analyses. Skewness reflects the degree of asymmetry of a distribution around its mean. Normal distributions produce a skewness statistic of about zero. Values that exceeded 2 standard errors of skewness were identified as significantly skewed (Sokal and Rohlf, 1995). Environmental variables with skewed distributions (gravel, grain-size EM2, smectite–chlorite, mica, K-feldspar) were log transformed and remaining parameters were left untransformed. To reveal intercorrelated parameters, we performed a variance inflation factor (VIF) analysis prior to ordination techniques to only retain non-correlated parameters in further multivariate analysis. Environmental variables with a VIF greater than 20 were eliminated, beginning with the variable with the largest inflation factor, until all remaining variables had values <20 (ter Braak and Šmilauer, 2012). A set of RDAs was performed on chironomid and diatom data with each environmental variable as the sole constraining variable. The percentage of the variance explained by each variable was calculated and the statistical significance of each variable was tested by a Monte Carlo permutation test with 999 unrestricted permutations. Significant variables (P≤0.05) were retained for further analysis. DCA, PCA, and RDA were performed using CANOCO 5.04 (ter Braak and Šmilauer, 2012).

Percentage abundances of the chironomid taxa that are absent or rare in the modern calibration data set were calculated at each sampling site in order to see the distribution of the taxa that could potentially hamper a T July reconstruction in the case of a palaeoclimatic study that could be done at each of the sampling sites. It is known that less reliability should be placed on the samples in which more than 5 % of the taxa are not represented in the modern calibration data or more than 5 % of the taxa are rare in the modern calibration data set (i.e. if the effective number of occurrences in the training set, Hill's N2, is less than five) (Heiri and Lotter, 2001; Hill, 1973; Self et al., 2011).

Species richness and the Simpson diversity on diatom and chironomid data were estimated after sample-size normalization using a rarefaction analysis of Hill numbers in the iNEXT package in R.

To assess the relative contribution of different sedimentary processes to the bulk sediment, such as fluvial or aeolian transport (Wang et al., 2015; Biskaborn et al., 2013b), a statistical end-member analysis on grain-size data was performed using the MATLAB modelling algorithm of Dietze et al. (2012). In this method, individual grain-size populations identified as end-member loadings (vol %, Fig. 4) as well as their contributions to the bulk composition identified as scores (%) were derived by eigenspace analysis, weight transformation, varimax rotations, and different scaling procedures.

Figure 4 End-member analysis grain-size distributions in 33 samples from Lake Bolshoe Toko.

Figure 5 Correlation matrix of selected environmental parameters. (a) Pearson correlation. Positive correlations indicated in red, negative correlations in blue. To keep the false discovery rate below 5 %, only p values of <0.001 were used to assign colours (Colquhoun, 2014). (b) Spatial autocorrelation associated with coordinates of sample sites and shown as p values generated by Moran's autocorrelation coefficient (R package “ape”).

A Pearson correlation matrix of the main important variables (Fig. 5a) was calculated using the basic R core (R Core Team, 2012) and plotted using corrplot. To keep the false discovery rate below 5 %, a p-value adjustment was applied prior to assignment of colours using only values that revealed p<0.001 (Colquhoun, 2014). To identify the pattern, the correlation matrix was reordered according to the correlation coefficient. Exceptional sites within the heterogenic lake system lead to disturbance of good correlation coefficients within areas along natural borders, e.g. water depth isobaths. Spatial autocorrelation of variables was estimated using latitudes and longitudes recorded of each sample site and displayed as p values generated by Moran's autocorrelation coefficient (R package “ape”).

To guarantee the sustained availability of our research (Elger et al., 2016), the data used in this study are freely accessible at the PANGAEA data repository (https://doi.org/10.1594/PANGAEA.906317).

4 Results

4.1 Water chemistry

Sampled surface waters of Bolshoe Toko (Table 1, ESM) were well saturated in O₂ (101 %–113 %) with a pH value in the neutral range (6.8–7.2). Electrical conductivity was very low for all waters (35.1–39.1 µS cm⁻¹), with slightly higher levels in the lagoon (67.8 µS cm⁻¹). Traces of Al (mean 72 µg L⁻¹), Fe (mean 46.6 µg L⁻¹), and Sr (mean 37.1 µg L⁻¹) were present, but there is no evidence of Pb, Cr, V, Co, Ni, or Cu. Mean sulfate concentrations (${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) were 2.35 mg L⁻¹ on average, with lower values in the lagoon (0.51 mg L⁻¹). The concentrations of nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) were 0.76 mg L⁻¹, but were lower in the lagoon (0.29 mg L⁻¹). ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ was 37.5 mg L⁻¹ in the lagoon and 14.9 mg L⁻¹ on average in the remaining samples. There was no phosphorus in any sample. Overall the water can be characterized as water of the Ca–Mg–${\mathrm{HCO}}_{\mathrm{3}}^{-}$ type.

Surface waters were characterized by mean isotope values of −18.7 ‰, −140.2 ‰ and 9.5 ‰ for δ¹⁸O, δD and d excess, respectively (n=6). The isotopic composition was relatively uniform in the main lake basin (δ¹⁸O = $-\mathrm{18.58}\pm \mathrm{0.15}$ ‰, δD = $-\mathrm{139}\pm \mathrm{0.7}$ ‰), while the lagoon (PG2122) exhibited slightly lower δ¹⁸O (δD) values of −19.2 ‰ (−145 ‰) (Fig. 6).

Figure 6 Hydrochemical situation between 2012 and 2013 in Lake Bolshoe Toko. (a) Profiles of water isotopes (δ¹⁸O) and temperature from different locations taken in August 2012 and March 2013. (b) δ¹⁸O/δD diagram for Lake Bolshoe Toko water samples. GMWL is the Global Meteoric Water Line (black line) and LMWL is the Local Meteoric Water Line for Yakutsk (dashed line; δD = 7.59⋅δ¹⁸O − 6.8) based on own data (monthly mean precipitation values between 1997 and 2006; n=106; Kloss, 2008).

In March 2013, isotope-depth profiles at PG2208 exhibited a slight isotopic enrichment trend from the surface to ∼5 m depth ($\sim +\mathrm{0.35}$ ‰ for δ¹⁸O), with a relatively uniform isotopic composition (δ¹⁸O = $-\mathrm{18.2}\pm \mathrm{0.2}$ ‰) below 10 m (Fig. 6a). These subtle variations likely reflect minor isotopic fractionation of surface waters during ice formation in spring and a well-mixed water column below. Conversely, the August 2012 depth profile at the western shoreline exhibited a gradual depleting isotope trend below ∼6 m depth, with marked variability that closely tracks water temperature changes (Fig. 6a). Meteorological data from the nearby weather station (Toko RS, 10 km northward) recorded heavy rainfall for August 2012 (25 mm above the long-term mean of 83 mm). Such precipitation events could cause temporary isotopic stratification or a variation in the isotopic signal throughout the water column. Due to ongoing mixing, these variations were then evened. In conclusion, variations in the isotopic composition throughout the August profile rather represent a temporal phenomenon and are not characteristic of Bolshoe Toko. In contrast, the lagoon showed a lighter isotope composition (δ¹⁸O = $-\mathrm{18.9}\pm \mathrm{0.2}$ ‰) than the main lake basin. All samples were positioned close to the Global Meteoric Water Line (GMWL, Fig. 6), indicating negligible evaporative effects on lake water isotope composition and a dominant influence of meteoric inputs both directly (i.e. precipitation) and indirectly (i.e. river inflows). The Local Meteoric Water Line for Yakutsk (dashed line; δD = 7.59⋅δ¹⁸O − 6.8), based on own data (monthly mean precipitation values between 1997 and 2006; n=106; from Kloss, 2008), is given for comparison and is indicative for more continental climate conditions.

4.2 Physicochemical sediment composition

The typical surficial lake bottom sediments consisted of either brown organic-enriched gyttja or sandy, organic-poor siliciclastic material. Sand contents ranged between 10.2 % and 96.2 % (mean 45.9 %, Fig. 7); silt contents ranged from 3.6 % to 83.3 % (mean 47.1 %); clay contents ranged from 0.2 % to 11.3 % (mean 5.8 %). Gravel was found only in four samples in the north-eastern near-shore areas, with contents of up to 13.1 %. The mean grain size ranged from 12 to 479 µm (mean 72 µm). The mean grain size generally correlated negatively with water depth (r −0.45). Mineral grains are composed mainly of quartz (32.7 %–76.2 %, mean 55.4 %), plagioclase (13.4 %–39.5 %, mean 26.2 %), K-feldspar (0.0 %–9.8 %, mean 5.6 %), and, to a smaller degree, pyrite (0.2 %–5.5 %, mean 3.3 %), hornblende (0.5 %–10.8 %, mean 3.1 %), mica (0.3 %–2.4 %, mean 1.1 %), and the clay minerals smectite, kaolinite, and chlorite (together 0.0 %–4.6 %, mean 2.0 %). The spatial distribution of minerals (Fig. 7) revealed a generally decreasing gradient of minerals relative to quartz starting from the Utuk River delta (proximal) towards the northern areas (distal).

Figure 7 Spatial distribution of the grain-size and mineral compositions of the surface sediments of Lake Bolshoe Toko. Maps compiled in ArcGIS 10.4. Scales chosen as 10 classes with equal intervals.

The CLR-transformed XRF data (Fig. 8) revealed high proportions of Zr and intermediate to high Ti near the Utuk River inflow and at the northern and eastern shore proximal areas. Zr values decreased with increasing water depth towards the lake centre with the exception of the shallow lagoon, where low values were observed. Mn values were highest in the lake centre and at the very deep site at the western steep subaquatic slope, and intermediate at shallow areas close to the shore. A minimum in Mn was found in the lagoon. Fe tends to be highest in the southern part of the lake basin, in the very shallow site in the north, and in the lagoon. Br showed a variable distribution; however, high values were found at two sites within the eastern lagoon and correspond to high TOC contents.

Figure 8 Spatial distribution of elements obtained from XRF measurements of surface sediments of Lake Bolshoe Toko. Maps compiled in ArcGIS 10.4. Scales chosen as 10 classes with equal intervals.

Additive log ratios (ALR) of Mn∕Fe were variable, with intermediate values found at sites surrounding the Utuk River inflow and low values within the lagoon and at basin central sites. High values were located at the deepest lake site as well as in the shallow north-eastern region. Both Sr∕Rb and Zr∕Rb ratios showed high values directly in front of the Utuk River inflow and decreased with distance toward the basin centre. Both Sr∕Rb and Zr∕Rb exhibited intermediate to high values in the north-eastern lake region and lower values in the lagoon. Si∕Ti ratio values demonstrated an increasing trend from the southern lake region and lagoon to the northern lake region.

The contents of total organic carbon (TOC, Fig. 9) range from 0.1 % to 12.3 % (mean 4.9 %). Maximum values occurred in the eastern area, intermediate values in the central basin, and lowest values in the northern shallow water areas. The difference between TOC and total carbon is within the error of the devices, and hence no inorganic carbon was detected. TOC contents and the TOC∕TN ratios were highest near the Utuk River inflow in the southern part of the lake, in the lagoon, and in proximity to the eastern shoreline. δ¹³C was measured in 15 samples and showed maximum values at the eastern shore (−25.7 ‰) and minimum values elsewhere (−27.8 ‰).

Figure 9 Spatial distribution of organic properties and statistical parameters inferred from diatom assemblages in the surface sediments of Lake Bolshoe Toko. Maps compiled in ArcGIS 10.4. Scales chosen as 10 classes with equal intervals.

Radiocarbon dating of a surface sample at site PG2139 (0–0.5 cm) indicated an age of 720±30 ¹⁴C yr BP (Lab-ID: Poz-105350, NaOH-SOL), while PG2207 (0–0.5 cm) suggested 1790±130 ¹⁴C yr BP (Lab-ID: Poz-105355, NaOH-SOL). Considering that the carbon concentration dissolved in sample PG2207 was too low (0.03 mgC), we use sample PG2139 as an estimated reservoir effect on the lake caused by the input of old carbon. Given that a hypothetical sediment surface is just a momentum only collectable as a range of past surfaces and there was more time available for radioactive decay at 0.5 cm depth than at 0 cm, the actual reservoir effect will be a little bit lower and should be confirmed by ²¹⁰Pb and ¹³⁷Cs measurements of downcore material before establishing an age–depth model for sediment cores.

4.3 Diatoms

The Bolshoe Toko diatom assemblages were characterized by boreal and arcto-alpine types and exhibited distinct spatial variations across the lake. In total, 142 different diatom taxa were found at 23 sites, dominated by planktonic species Pliocaenicus bolshetokoensis (Genkal et al., 2018) (0.0 %–27.9 %, mean 14.7 %), Cyclotella comensis (0.0 %–23.1 %, mean 10.9 %), and benthic species Achnanthidium minutissimum (0.0 %–38.0 %, mean 11.8 %). The relative content of planktonic species (Fig. 9) was 2.0 %–73.7 % (mean 54.2 %), epiphytic species 19.2 %–83.9 % (mean 36.4 %), and epibenthic species 2.6 %–23.0 % (mean 9.3 %). The spatial distributions of the main taxa are presented in Fig. 10. Small benthic fragilarioid species were represented by 0.0 %–27.6 % (mean 7.4 %), Naviculoid species ranged from 3.3 % to 12.9 % (mean 7.2 %), and Aulacoseira species ranged from 0.0 % to 10.8 % (mean 4.5 %). Pliocaenicus bolshetokoensis maximum abundance occurred in areas of deepest water such as the southern part of the lake and in the eastern lagoon. Cyclotella species were more abundant in the central lake and were not as strictly bound to water depth as Pliocaenicus. Aulacoseira species displayed no clear spatial pattern, but were less abundant in the northern shallow water areas. Tabellaria species were more abundant in shallow near-shore areas than in central and deep-water areas.

Achnanthoid (monoraphid) species were most abundant in near-shore areas, especially near the eastern lake terrace. Fragilarioid (araphid) species were common in the southernmost part near the inflow, as well as the lagoon. Other benthic species, i.e. Navicula, Cymbella, and Eunotia, were generally more abundant in shallow near-shore areas than in deeper water areas.

In pelagic areas planktonic diatoms were generally more abundant than epiphytic and epibenthic species. Epiphytic species, however, predominated in some shallow areas in the northern and eastern parts of the lake. Epibenthic species occurred in smaller abundances in shallow lake littorals. Together with an increased amount of non-planktonic species, the Simpson diatom species diversity was higher in the northern and eastern parts of the lake. The chrysophyte index was high near the river inflow in the south and along the river-like bathymetrical structure, as well as the lagoon where another small river inflowed into the lake. The Mallomonas index, reported for high nutrients and low pH (Smol et al., 1984), was highest near the inflow and in the central part and lowest in near-shore areas in the north and east. The maximum f-index value, representing the highest valve preservation, was found in the near-shore areas, whereas lower values were found at the shallow bathymetrical structure in the central part of the lake. Maximum valve concentrations were observed in the central and northern lake basin.

The initial RDA with all environmental variables indicated that axes 1 and 2 explain 39.6 % of variance in diatom species data. After deleting all intercorrelated variables, 13 parameters with VIFs <20 were left for manual selection with Monte Carlo test. The analysis revealed eight statistically significant (p≤0.05) explanatory variables: TOC∕TN, TOC, water depth, distance from river, distance from the shore, presence of vegetation, sand, and EM3, (ESM diatoms, Fig. 2). Eigenvalues for RDA axes 1 and 2 constrained by eight significant environmental variables constitute 81 % and 59 %, respectively, of the initial RDA, suggesting that the selected significant variables explain the major variance in the diatoms data. The RDA biplots of the species scores and sample scores (Fig. 2) show that diatom species and sites are grouped according to the main environmental forcing responsible for their spatial distribution. The clearest environmental signals in the RDA are related to water depth, habitat preferences and river influence. The upper left quarter of the biplot is strongly influenced by water depth, grain size (EM3), and the ratio between TOC and TN. The species found next to water depth are planktonic Cyclotella taxa, whereas Aulacoseira is closer to TOC∕TN and the total carbon content. In the lower right quarter epiphytic and benthic taxa prevail, i.e. achnanthoid, naviculoid, and cymbelloid taxa, associated with the presence of vegetation and coarser (sand) substrate conditions. The distances to river and to shore are crossing the lower left quarter and are associated with different planktonic Cyclotella and achnanthoid taxa, while in the opposite direction, with increasing Utuk River influence, fragilarioid taxa, Eunotia, Tabellaria, and Gomophonema prevail, next to the high influence of TOC∕TN.

Mean surface sediment δ¹⁸O_diatom was 22.8 ‰ (minimum 21.9 ‰, maximum 23.6 ‰, n=9, Fig. 9) with a standard deviation of ±0.6 ‰ (1σ). The spatial distribution indicated higher values ∼23.3 ‰ in the deeper south-western part of lake (PG2113, 2115, 2005) and lower values ∼22.3 ‰ in the shallower northern lake basin (PG2139, 2140, 2147, 2209). The two samples from the lagoon exhibited values of 22.2 ‰ in the shallower northern area and 23.6 ‰ in the deeper part. Four samples from the southern part could not be purified well enough and had contamination corrections >2 ‰.

Figure 10 Spatial distribution of the main diatom taxa in the surface sediments of Lake Bolshoe Toko. Maps compiled in ArcGIS 10.4. Scales chosen as 10 classes with equal intervals. Maps e and h had exceptionally high values of achnanthoid and cymbelloid taxa only in the very shallow (0.5 m) site PG2142. These values are shown in purple, indicated separately on the right-hand side of the scales.

4.4 Chironomids

A total of 79 different chironomid taxa were present in the surface sediment samples, of which 48 belong to the subfamily Orthocladiinae, 25 to Chironominae (15 from the triba Tanitarsini and 10 from the triba Chironomini), 4 to subfamily Diamesinae, and 2 to Tanypodinae.

Figure 11 Spatial distribution of chironomid taxa and inferred statistical parameters in the surface sediments of Lake Bolshoe Toko. Maps compiled in ArcGIS 10.4. Scales chosen as 10 classes with equal intervals.

The initial RDA with all environmental variables shows that axes 1 and 2 explain 46.7 % of variance in the taxon data. Most of the environmental parameters were intercorrelated, and following sequential deletion of all redundant variables, eight parameters with VIFs <20 remained for the further analysis. The manual Monte Carlo test selection demonstrates four statistically significant (p≤0.05) explanatory variables: TOC∕N, water depth (WD), distance from river, and presence of vegetation (Table 2). Distance from the river and presence of vegetation showed lower than TOC∕N and WD level of significance. However, we still use these parameters for interpretation of the chironomid data, as there was a clear gap between the four chosen parameters (p=0.001 to 0.059) and much higher p values (>0.25) of the following tested parameters (TC, distance to the shore, silt, clay). Eigenvalues for RDA axes 1 and 2 constrained by four significant environmental variables were 0.200 and 0.150, respectively, and constituted 70 % and 85 % of the RDA performed on all environmental variables (0.289 and 0.177, respectively). This minor difference suggests that the four selected variables sufficiently explain the major gradients in the chironomid data.

The RDA biplot of the sample scores shows that sites are grouped by their location in relation to the major environmental variables (Fig. 11), and distribution of chironomid taxa along the RDA axes reflects their ecological spectra. Figure 11 and Table 6 in the ESM show median values of eco-taxonomical chironomid groups and their relation to environmental parameters.

Sites most strongly influenced by the inflowing rivers grouped in the lower left quadrant of the biplot, as the vector in the upper right quarter shows increasing distance from the river mouth. In total 64 chironomid taxa were found in this group of sites, and of these 33 were only found here. Chironomid fauna were chiefly represented by phytophilic littoral taxa from the Orthocladiinae genera Cricotopus, Orthocladius, Eukiefferiella, and Parakiefferiella etc. (Fig. 11). Another important feature is the presence of a relatively high amount of lotic environmental taxa, among which are several Diamesa taxa, Rheocricotopus effusus-type, Synorthocladius, Brillia, and for lotic–lentic environments Parakiefferiella bathophila-type, P. triguetra-type, Nanocladius rectinervis-type, N. branchicolus-type, several Eukiefferiella taxa, and Stictochironomus.

The group in the opposite upper right quadrant represents the northern part of the lake situated far from the inflowing rivers. Here, mainly profundal taxa prevail, i.e. Procladius, Polypedilum nubeculosum-type, Cryptochironomus (eurytopic), and Heterotrissocladius maeaeri-type 1 (acidophilic).

The lower right group of sites represent eastern shallow littoral with presence of macrophytes. Species richness and proportion of semiterrestrial and littoral taxa in this group is generally low. Littoral taxa were generally phytophilic: Cricotopus intersectus-type, C. cylindraceus-type, Dicrotendipes nervosus-type (mesotrophic), and Cladotanytarsus mancus-type and Psectrocladius sordidellus-type (acid-tolerant mesotrophic). Most abundant profundal taxa here belong to the acid-tolerant Heterotrissocladius genera represented by H. macridus-type, H. maeaeri-types 1 and 2, H. grimschawi-type (acidophilic), and to the subfamily Tanypodinae represented by Procladius. The sites grouped in the opposing upper left quadrant represent lotic and lotic–lentic taxa (Diamesinae, Thenimaniella clavicornis-type, Eukiefferiella claripennis-type, Eukiefferiella fittkaui-type, several Orthocladius taxa).

5 Discussion

5.1 Spatial control of abiotic and biogeochemical sediment components

Sediment-geochemical and physical properties of the uppermost surface of the sediment basin in Bolshoe Toko are spatially variable. Physical properties of particles within the surface sediments depend chiefly on transportation processes and the characteristics and availability of clastic compounds in the lake catchment. The main catchment comprises the Stanovoy mountain range in the south channelled through the Utuk River into Bolshoe Toko. Accordingly, the lake experiences annual input of suspended material through a single source at the Utuk River mouth that likely is at its maximum during spring snowmelt (Bouchard et al., 2013). The grain-size data and their end-members (Figs. 4 and 7) indicate that the relative proportions of sand, silt, and clay are somewhat constant in proximity to the Utuk River inflow but change towards the north and at the lake shoreline. Whereas in the central northern lake basin the amount of silt increases, the proportions of sand increase along the northern shoreline on top of the drowned moraine. Gravel is only present in samples near the lake terraces in the east. The constant distribution in the southern–central lake basin reflects the river input. Decreasing river influence and hence decreasing water transport energy with increasing distance from the river mouth leads to the observed predominance of finer grain-size (silt dominated) samples in the northern central parts of the lake. Sandy samples along the shoreline reflect direct input from the moraines around the northern part of the lake. Other relevant within-lake sedimentary processes include shore erosion and inwash and winnowing of fine sediment grains by surface currents as well as alluvial processes and debris flows which continue basin ward as subaquatic flows. The restriction of gravel at the eastern shore can be attributed to the availability of source material and suitable transport pathways of coarser clasts from the third moraine. In consequence to the described lateral transport trajectories and local control factors within the lake, there is only weak negative correlation between mean grain size and water depth (r −0.45, Figs. 7 and 12).

The modelled end-member loadings of the observed grain-size classes (Figs. 4 and 7) indicate an EM1 major peak in fine silt that represents fluvial sediment input. EM2 has peak values in fine to medium sandy grain-size fractions and in the northern part of the lake indicative of depositional processes associated with the erosion of moraines distal from the river inflow, where the hydrological dynamics in the lake basin are weak. The weak positive correlation between EM3 and the concentration of diatom valves (r 0.44) likely represents both in situ diatom valves that could not be removed from allochthonous sediment particles during sample processing, and possibly redistributed ice-rafted debris (Wang et al., 2015).

Intermediate concentrations of TOC and high ratios of TOC∕TN in the south as compared to the north suggest differences in catchment characteristics, i.e. a considerable allochthonous contribution of terrestrial plant material from the Utuk River. This assumption is supported by previous findings that show non-vascular plants, i.e. phytoplankton and other algae, with TOC∕TN ratios between ca. 5 and 10, while organic matter from vascular land plants has higher values of about 20 (Meyers and Teranes, 2002). High values of TOC∕TN in lake sediment surfaces at river inflows have also been observed in other studies (Vogel et al., 2010). δ¹³C is generally low on average (−26.8 ‰) and only slightly higher at the eastern shore (−25.7 ‰), suggesting a strong overall dominance of C₃ plants and phytoplankton in the bulk organic matter fraction (Meyers, 2003). It remains unclear as to the degree of old and reworked organic carbon, e.g. from charcoal deposits, transported to the lake.

The distribution of elements from the XRF scanning data suggests strong abiotic relationships with grain-size and mineral distributions. We focus on heavier elements because lighter elements, even though commonly in higher concentrations, show potential contribution from multiple sources. Sr∕Rb ratios and Zr are negatively correlated with kaolinite and chlorite (r −0.73 and −0.85, respectively). As described in Kalugin et al. (2007), Rb substitutes for K in clay minerals. The Sr∕Rb ratios do not however show a significant correlation with grain-size parameters, as found in other studies (Biskaborn et al., 2013b). We assume therefore that Sr, as a substituent for Ca, is influenced by multiple minerals represented in different grain-size fractions, i.e. K-feldspar (r 0.45) and Hornblende (r 0.24). Associated with high metamorphic grades in the Stanovoy mountains, Sr is preferentially taken into the K-feldspar phase (Virgo, 1968). Conversely, the Zr∕Rb ratio correlates well with the sand fraction (r 0.50) and with the mean grain size (r 0.49), but negatively with silt (r −0.54) and clay (r −0.39). We account for this effect by a higher diversity of minerals in the input of the Utuk River supplying the lake basin with mafic Ca-rich metamorphic rocks from the Stanovoy mountains. The strong influence of the Utuk River in the spatial distribution of physicochemical sediment components is further demonstrated by the decreasing gradient of minerals relative to quartz starting from the Utuk River towards the northern lake basin (Fig. 7). The most representative indicator of grain-size variations in surface sediments is given by CLR-transformed values of Ti, which correlate well with the sand fraction (r 0.74) and the mean grain size (r 0.88).

Si∕Ti ratios have traditionally been used as a proxy for the biogenic silica content of sediments (Melles et al., 2012). This stems from the fact that Ti is generally attributed to detrital influx and Si to both detrital and biogenic (diatom) origins. At Bolshoe Toko positive correlations between Si∕Ti ratios, diatom valve concentrations (r 0.36) and the ratio of planktonic to benthic diatoms (r 0.42) suggests that Si∕Ti may be useful to trace the relative portion of diatom valves in intermediate grain-size fractions. Moreover, the Si∕Ti ratio correlates significantly with silt (r 0.81).

Mn∕Fe ratios have been ascribed to redox dynamics associated with bottom water oxygenation processes (Naeher et al., 2013). In Bolshoe Toko, however, the detrital input of ferrous minerals, i.e. pyrite, suggests that Mn∕Fe ratios cannot be directly attributed to redox processes in the surface sediments. This is supported by the correlation of Fe with the sand fraction (r 0.6) and grain size (r 0.59). Accordingly, we found no significant correlations between Mn∕Fe and other abiotic or biotic proxies.

Lastly, there is an uncertainty in the spatial distribution of elements measured by XRF techniques. We attribute this lack of clear patterns to (1) methodological hurdles to apply XRF techniques to surface sediments commonly rich in water and organic material, and (2) multiple sources of the same elements coming from minerogenic input, grain-size differences in individual samples and different intensities of redox processes at different habitat settings. The high variance of elements are therefore representative of the high complexity of this lake system, rather than unequivocal validations or falsifications of the applicability of XRF scanner data as an environmental proxy at Bolshoe Toko.

5.2 Factors explaining the spatial diatom distribution

Diatom communities in Yakutia respond rapidly to environmental changes including hydrochemical parameters, water depth, nutrients, and catchment vegetation type (Pestryakova et al., 2018). Planktonic diatom species are ubiquitous across Bolshoe Toko, with a distinct tendency of the ratio between planktonic and benthic species to greater water depths (r 0.77, Figs. 5 and 12), due to the limited availability of light for benthic species (Gushulak et al., 2017; Raposeiro et al., 2018). Especially Aulacoseira species were never abundant along the shallower northern and eastern shorelines. The primary difference between the two most abundant genera in the lake is that Pliocaenicus exhibits the highest abundances proximal to the inflow and in the south-eastern lagoon, whereas Cyclotella are more abundant in the lake centre and absent in the lagoon. Little is as yet known about the new species Pliocaenicus bolshetokoensis (Genkal et al., 2018). Our findings suggest factors other than water depth (r 0.39), such as proximity (e.g. nutrient supply) to the Utuk River and small streams, as controlling parameters for bloom intensities of this species. Cyclotella, however, are restricted to stratification of the water column and hence are more abundant at distance from the river mouth, where incoming water causes turbulence (Rühland et al., 2003; Smol et al., 2005). Cyclotella are therefore also believed to benefit from recent air temperature warming trends and will likely increase in abundance (Paul et al., 2010). Aulacoseira is a dense, rapidly sinking tychoplanktonic group of species requiring water turbulence to remain in the photic zone (Rühland et al., 2008, 2015), which explains the lower abundances in the northern and hydrologically less dynamic zones within the lake. Lightly silicified Tabellaria species are known to occur in zigzag planktonic colonies, yet they also appear as short-valved populations in the benthos (Lange-Bertalot et al., 2011; Biskaborn et al., 2013a; Krammer and Lange-Bertalot, 1986–1991). In Bolshoe Toko, the spatial distribution of Tabellaria indicates benthic habitats are more favourable than planktonic.

Figure 12 Distribution of grain size, organic carbon and nitrogen indices, diatom and chironomid parameters, and selected elements and minerals in dependence on water depth in Lake Bolshoe Toko.

The most common non-planktonic species in Bolshoe Toko belong to achnanthoid (monoraphid) genera, of which most species are epiphytic. Epiphytic species exhibit a stronger negative correlation with water depth (r −0.68) than epibenthic species (r −0.4), indicating that aquatic plants, in turn controlled by water transparency, pH, water depth and nutrient status (Valiranta et al., 2011), have an important function in the lake ecosystem (Fig. 12). The highest abundance of achnanthoid and cymbelloid valves occurs at 400 m distance to the northern shore at a water depth of 0.5 m.

Fragilarioid species are adapted to rapidly changing environments and are thus good indicators of ecosystem variability (Wischnewski et al., 2011). The peak occurrences of Staurosira species, which are pioneering small benthic fragilarioids (Biskaborn et al., 2012), therefore indicates the formation of a new ecosystem habitat type in the lagoon at the south-eastern lake basin. We assume this basin is successively separated from the main basin and will eventually form a small isolated remnant lake, similar to “Banya” lake (Fig. 1). High productivity of epiphytic species and low detrital input suggested by elemental and grain-size data, together with higher organic content (high TOC and Br), indicate a calm sedimentological regime with high bioproductivity. Similar neutral pH values measured in water samples from the central basin and the lagoon (Table 1) questions pH as a main driving factor of the Eunotia peak in the lagoon. However, Barinova et al. (2011) suggest a 5.0–5.8 pH range for the identified Eunotia species, which rather indicates that the pH values obtained during April in 2013 are not representative for the annual average and the specific catchment of the lagoon, which likely will differ from this point measurement. The ice break-up during spring and transport of water from the catchment restricted to the lagoon likely leads to milieu differences in the lagoon relative to the main basin.

High autocorrelation coefficients (Moran's I p values) for species richness and valve concentration indicate strong local influence of biotic processes, i.e. reproduction, leading to spatial autocorrelation (Legendre et al., 2005). The lowest observed autocorrelation for the diatom planktonic ∕ benthic ratio confirms the strong relationship between diatom species assemblage composition and water depth. A strong relationship between diatom diversity and water depth is supported by a study comparing morphological count data and phylogenetic species data gained by next-generation sequencing DNA analysis (Stoof-Leichsenring et al., 2019).

The RDA biplot of diatoms (Fig. 2) suggests that both water depth and distance to river are important lake attributes accounting for the species distributions across the lake. Especially Eunotia, fragilarioids, Tabellaria, and also Aulacoseira subarctica appear more frequently at sites that are close to the Utuk River mouth (e.g. PG2113, PG2115, PG2117, and PG2118). The high TOC∕TN ratios in these samples illustrates the strong riverine input of allochthonous material. In the biplots, high water depth is primarily associated with Cyclotella species (and Aulacoseria), while Aulacoseira species tend to be additionally influenced by incoming rivers and also thrive closer to the shorelines. Areas close to river mouths are usually dominated by river taxa and species that prefer higher nutrient content related to river input and associated early ice-cover melting (Kienel and Kumke, 2002). Accordingly, the influx of diatoms from wetlands in the lake catchment is an important additional factor influencing the spatial diatom distribution (Earle et al., 1988). Compared to direct conductivity, water depth and nutrient controls, the link between temperature and diatom species is poorly understood in Yakutian lake systems (Pestryakova et al., 2018) and should be avoided.

Our RDA also shows that a high diversity of benthic, and particularly epiphytic diatom species, i.e. several achnanthoid species and some naviculoid taxa, plot in the opposite direction from water depth together with vegetation and the coarse grain-size fraction. Kingston et al. (1983) revealed spatial diatom variability in the Laurentian Great Lakes, where the stability of diatom assemblages increased with water depth. In shallower marginal waters of the Great Lakes, the availability of diverse habitats, including benthic and periphytic niches, leads to high species diversity. According to our data in Bolshoe Toko, the Simpson diversity index suggests higher effective numbers of dominant species associated with increased habitat complexity (Kovalenko et al., 2012), i.e. availability of water plants and benthic substrates in shallower depths along the eastern and northern shores. Thus, higher diversity in this area is facilitated by differential catchment preferences. However, it can be assumed that due to lesser water supply rates from the small northern part of the catchment (Fig. 1), a single location at the north-eastern lake margin will likely not receive significantly higher loadings of nutrients as compared to the Utuk River coming from the igneous mountain range. Nevertheless, moraine deposits typically contain high amounts of silt and clay which can more easily be weathered and altered to fertilizing substances that are transported into the calm and shallower northern part of the basin.

The indices of chrysophyte cysts and Mallomonas relative to diatom cells exhibit indistinct patterns in spatial distributions but a slight tendency towards proximity to river input and high water depths. Although chrysophyte cysts commonly represent planktonic algae (Smol, 1988b), periphytic taxa are also common in boreal regions (Douglas and Smol, 1995) with cool and oligotrophic conditions (Gavin et al., 2011). Mallomonas was reported as an indicator of lake eutrophication and acidification (Smol et al., 1984).

Taphonomic effects on the preservation of subfossil assemblages are generally influenced by clastic transport mechanisms depending on the lake morphology (Raposeiro et al., 2018). The preservation of diatom valves in Bolshoe Toko is found to be lowest in samples from a plateau-like feature in the central part of the lake bottom, which indicates increased re-working associated with bottom currents and/or increased dissolution of diatom valves due to lesser accumulation rates and/or increased grazing activity of herbivorous organisms (Flower and Ryves, 2009; Ryves et al., 2001).

The spatial distribution of δ¹⁸O_diatom from the sediment surface indicates higher δ¹⁸O_diatom values at the deeper, south-western part of the lake with a difference of approximately 1 ‰ compared to lower δ¹⁸O_diatom values in the shallower northern part. This could reflect a combination of spatial δ¹⁸O_water variations, water temperatures, and/or a potential species-driven fractionation effect. However, existing studies demonstrate no apparent species composition effects on lacustrine δ¹⁸O_diatom (Bailey et al., 2014; Chapligin et al., 2012b). Additionally, the sieving step reduces the assemblage before the isotope analysis to a small size interval, resulting in a similar species composition. Furthermore, dissolution effects in nature and during sample preparation could have an impact on δ¹⁸O_diatom. However, we suppose differential dissolution to have a minor influence on the spatial variability of δ¹⁸O_diatom at BT samples tackled in our study as these are (1) of similar age, (2) have been treated with wet chemistry at low temperatures, and (3) after preparation do not show any microscopical signs of dissolution effects, i.e. a low diatom dissolution index (Smith et al., 2016).

Regarding δ¹⁸O_water variability, waters sampled at the same time in different parts of the lake show a uniform isotopic composition (within ±0.15 ‰) and indicate an isotopically well-mixed lake. Considering this is a one-time recording, slight seasonal variation between shallower and deeper parts (for example due to evaporation) cannot be excluded and could account for some differences in ¹⁸O. However, lake surface evaporation would result in isotopic enrichment and overall higher δ¹⁸O_diatom values.

Alternatively, the lake temperature in which the diatoms grow has an impact of ca. −0.2 ‰ ^∘C⁻¹ on δ¹⁸O_diatom (Brandriss et al., 1998; Dodd et al., 2012; Moschen et al., 2005). Shallower areas heat up faster, especially in the photic zone. The temperature profile near to the western shoreline taken in August 2012 (Fig. 6) shows 12 ^∘C at the surface, with an average of approximately 10 ^∘C in the first 15 m of the water column decreasing to approximately 6 ^∘C in 30 m depth. Although a spatial difference of 5 ^∘C in the photic zone for causing a 1 ‰ shift is rather unlikely, this could account for part of the variation in surface δ¹⁸O_diatom.

5.3 Factors explaining the spatial chironomid distribution

The chironomid RDA indicates that spatial variations are primarily influenced by the distribution of tributary rivers. For example, high species diversity is found adjacent to the Utuk River inflow (2117) and in the south-eastern lagoon fed from a small inflowing stream (PG2122). Semiterrestrial taxa, like Smittia-Parasmittia, Pseudosmittia, and Limnophies-Paralimnophies, have been found only here, with the highest abundances of 6 % and 3.2 % at the sites opposite of the inflowing rivers (PG2117 and PG2122), suggesting these taxa were transported from marshy river deltas.

Species at lentic sites with no tributary influence are primarily controlled by water depth. Deep profundal sites of the lake have much lower taxonomic richness in chironomid communities. Higher taxonomic richness at site PG2118 can be explained by an enriching riverine influence. High proportions of lotic and lotic–lentic taxa lead to a high taxonomic similarity of this profundal site to littoral sites in the south and south-east. Similarly, in relation to temperature, sublittoral and profundal sites both have much higher representation of the taxa characteristic of semi-warm conditions and lower abundances of the taxa preferring warm and cold conditions. However, high depths of the sublittoral and profundal sites lead to the development of a poor chironomid fauna at these sites. High distance from the shore and presumably only weak transportation of chironomid remains of littoral fauna to the profundal zone could be another limiting factor for diversity of chironomid communities in the profundal zone.

Eastern relatively shallow littorals are inhabited by more diverse, phytophilic, mesotrophic, and partly acidophilic fauna with absence of lotic taxa, related to a less disturbed and turbulent environment and presence of macrophytes. This fauna has higher abundance of the semi-warm and warm taxa. The presence of mesotrophic to eutrophic and acidophilic taxa can be attributed to paludification of the shore zone and decomposition of macrophytes and submerged vegetation in the shallow littoral (Nazarova et al., 2017b).

It is still debated how spatial and local environmental processes influence the distribution of chironomids at a small spatial scale in a lake (Luoto and Ojala, 2018; Yang et al., 2017). It is known that within one water body the concentration of chironomid head capsules can vary from zero to several thousand per 1 cm³ of sediments (Kalinkina and Belkina, 2018; Walker et al., 1997) depending on factors such as water depth, rate of sediment accumulation, the hydrological conditions, or anthropogenic influence. Water depth in particular is a major driving factor of chironomid assemblages (Ali et al., 2002; Luoto, 2012; Vemeaux and Aleya, 1998), with depth optima of several species consistent across broad spatial scales (Nazarova et al., 2011). Chironomid remains from the deepest zones of Bolshoe Toko represent an assemblage of elements of profundal necrocenosis (Hofmann, 1971) mixed with secondary components of littoral fauna transported with in-lake hydrological and sedimentary processes into the profundal zone from outside. Thus, the re-deposition of littoral taxa into the profundal zone is an important factor that affects the final composition and abundance of subfossil assemblages. While in small lakes, subfossil assemblages from the profundal zone quite adequately reflect the fauna of the entire water body (Brooks and Birks, 2001; Walker and Mathewes, 1990), our findings support the hypothesis that in large lakes the taphonomy of chironomid communities seems to be more complex (Yang et al., 2017; Árva et al., 2015).

5.4 Lake Bolshoe Toko as a site for palaeoclimate reconstructions

Compared to small lowland lakes of central and northern Yakutia, sedimentary processes are quite different in Bolshoe Toko. One reason is the lack of thaw slumps, subsidence, and other permafrost-related phenomena (Biskaborn et al., 2013b) that are typical of shallow thermokarst lake settings across northern permafrost regions (Biskaborn et al., 2016, 2012, 2013a, b; Bouchard et al., 2016; Schleusner et al., 2015; Subetto et al., 2017).

The Bolshoe Toko mineral composition is primarily influenced by the Utuk River, and only samples in extremely shallow areas are influenced by direct shoreline input. The grain-size signal is influenced by dissolution effects associated with organic matter and in situ growth of diatom valves. Conversely, the coarser fractions parallel minerogenic compositions and water depth. Accordingly, the grain-size distribution originated from multiple processes and should only be considered an environmental proxy when combined with biotic indicators.

Diatoms are spatially distributed according to their preferred habitat. Aside from the spatial habitat conditions associated with basin morphology, an additional consideration is the annual duration and thickness of lake ice cover (Keatley et al., 2008; Smol, 1988a). For instance, planktonic communities in Lake Baikal, including Aulacoseira species, are found to grow under the ice if the surface snow properties (i.e. thickness, density) allow sufficient light penetration (Jewson et al., 2009; Mackay et al., 2005). Generally, planktonic and benthic diatom species have strategies to survive in ice-covered lakes by growing in benthic mode, forming resting spores, or attaching to the ice-cover substrate (D'souza, 2012). Hence, the duration and presence of ice cover can significantly impact both changes in assemblage composition and spatial distribution, particularly including the ratio of planktonic to benthic diatoms (Wang et al., 2012a; Bailey et al., 2018).

The applicability of chironomids for temperature reconstructions reveals clear spatial constraints; 22 % of the taxa in sites with riverine influence are absent or rare from the FE mean July chironomid-based temperature inference model (Nazarova et al., 2015), whereas fewer of these rare/absent taxa occur in the central and northern littoral, sublittoral, and profundal parts of the lake (Fig. 4). However, low taxonomic richness of the profundal zone also hampers palaeoclimatic inferences. Also, the number of chironomid head capsules are generally lower here relative to littoral sites. Maximum taxonomic diversity in areas influenced by lake tributaries can be explained by both a taxonomic enrichment from the lake catchment as well by more favourable oxygen and nutrient conditions.

The applicability of δ¹⁸O_diatom as a proxy of past hydroclimate conditions at Bolshoe Toko is facilitated by the main controls influencing δ¹⁸O_diatom, which are here found to be (1) lake water temperature (T_lake) and (2) lake water isotope composition (δ¹⁸O_lake) (Dodd and Sharp, 2010; Leng and Barker, 2006; Labeyrie, 1974; Leclerc and Labeyrie, 1987). The fractionation between lake water and biogenic opal can be calculated when comparing δ¹⁸O_lake (mean: −18.7 ‰) with recent surface sediments of Lake Bolshoe Toko and their respective mean δ¹⁸O_diatom (of +22.8 ‰) using this isotope fractionation correlation between sedimentary diatom silica and water as determined by Leclerc and Labeyrie (1987). The mean T_lake can be estimated to ca. 6 ^∘C for the photic zone/diatom bloom. This estimate is at the lower end of summer temperatures between 4.8 and 12 ^∘C. The corresponding derived mean isotope fractionation factor for the system diatom silica–water α=1.0424 matches the fractionation factor for sediments proposed by Leclerc and Labeyrie (1987) well (α_{(silica–water)}=1.0432).

Additionally, as lacustrine δ¹⁸O_diatom also reflects the isotopic composition of the water where the diatoms grow (δ¹⁸O_lake), δ¹⁸O_diatom typically reflects meteoric inputs associated with precipitation and riverine inflows (Fig. 6b). For example, existing studies have used lacustrine δ¹⁸O_diatom to reconstruct past changes in precipitation amount and seasonality, the precipitation–evaporation balance, spring snowmelt inputs, and synoptic-scale shifts in atmospheric circulation (Bailey et al., 2015, 2018; Meyer et al., 2015; Kostrova et al., 2013; Mackay et al., 2013). It is envisaged that changes in δ¹⁸O_diatom through time at a single site in Bolshoe Toko will yield insights into the long-term air temperature and palaeohydrological history of the region.

Positive feedback mechanisms between benthic algae and chironomid larvae in benthic ecosystems are well documented (Herren et al., 2017). Chironomids in Bolshoe Toko, however, showed less significant correlations with benthic diatom species but weak correlations with planktonic species and lake attributes associated with benthic habitats and water depth, highlighting the potential of chironomids for independent water depth and temperature reconstruction in future sediment core studies (Nazarova et al., 2011).

High correlation coefficients between organic carbon and Pliocaenicus bolshetokoensis (0.66) and silt (0.65) suggest that the accumulation of organic matter and intermediate grain-size fraction is, to a certain degree, controlled by the productivity of siliceous microalgae (Biskaborn et al., 2012). A strong contribution of plankton indicates that TOC∕TN ratios can provide insights into the relative influx between land and water plants (Meyers and Teranes, 2002). The relatively weak correlation between TOC∕TN ratios and water depth (r 0.51) demonstrates the accuracy limits of TOC∕TN as a proxy for relative lake-level changes. This is caused by transport and accumulation of allochthonous organic matter in proximity to the Utuk River. Furthermore, correlations between TOC∕TN and TOC, as well as negative correlations with grain-size indicators, suggest diagenetic alteration (i.e. loss) of nitrogen in the surface sediments (Galman et al., 2008).

The distinct difference between two samples along the subaquatic slope near the western shore (diatoms, minerals, organics) indicates redistribution of sediment. Downslope transport of surface layers over the time could lead to redistribution of old material into the deepest parts of the basin. Due to higher accumulation rates, a sediment core from the deepest part of the basin would potentially provide a higher temporal resolution but also a higher risk of repositioned sediment layers. On top of redistribution processes, hump-shaped relations between lake depth and species diversity observed in other studies suggest that the total subfossil species assemblages are better represented at intermediate depths than at the maximum depth (Raposeiro et al., 2018). A coring site at intermediate depth in the shallow northern and sedimentologically calm sector of the basin would enable the tracking of different river and glacial influences and offers greater chances of undisturbed successions of bioindicator time series.

6 Conclusions

Our study on the within-lake variance of environmental indicator data and its attribution to habitat factors improves the understanding of lake-internal filters between environmental forcing and the resulting sediment parameters of Lake Bolshoe Toko and comparable boreal, cold, and deep lakes. We found that the spatial variabilities of biotic ecosystem components are mainly explained by static habitat preferences such as water depth and river distance. Abiotic sediment features are not symmetrically distributed in the basin, but vary along restricted areas of differential environmental forcings (e.g. river input, rocky shore, steep shore, and shallow shore). They depend, in addition to water depth and riverine activity, on multiple interacting factors, such as catchment characteristics, geochemical sediment diagenesis, and hydrochemical dynamics. Our main findings can be highlighted as follows.

The lake water of Bolshoe Toko can be characterized as Ca–Mg–HCO₃-type water. It is well saturated in O₂, neutral to slightly acidic, showing a low conductivity and corresponding ion concentrations suggesting unpolluted freshwater conditions. Lake Bolshoe Toko is a cold, polymictic, oligotrophic, open through-flow lake system and can be regarded as an undisturbed ecosystem.

Water depth is a strong factor explaining the spatial variability of diatoms and chironomids. The proportions of planktonic to benthic diatoms and profundal to littoral chironomids serve as a reliable lake-level proxy.

The diatom assemblage is dominated by planktonic species, i.e. Pliocaenicus bolshetokoensis, which is unique for this lake, and more common plankton such as Cyclotella and Aulacoseira, as well as non-planktonic taxa, such as Achnanthidium. Diatom species richness and diversity are higher in surface sediments in the northern part of the basin, associated with shallower waters and the availability of benthic and periphytic niches.

The δ¹⁸O_diatom values (22.8±0.6 ‰) show slight spatial variations with higher values in the deeper south-western part of the lake probably related to water temperature differences in the photic zone during the main diatom bloom. The silica–water isotope fractionation is suitable for further downcore investigations for assessing palaeohydrological information and potential air-temperature changes in the region.

The water of Bolshoe Toko is well mixed and does not show significant isotopic stratification apart from lake ice-cover formation where thermal stratification prevents mixing. The isotopic lake water composition (δ¹⁸O = 18.2±0.2 ‰) corresponds to the GMWL and does not show evaporative enrichment. Both isotopic and hydrochemical data indicate atmospheric precipitation (and meltwater run-off) as the main water source. Accordingly, δ¹⁸O_lakewater is directly linked to δ¹⁸O_{precipitation}.

The highest amount of the chironomid taxa underrepresented in the FE training set used for regional palaeoclimate inference was found close to the Utuk River and at southern littoral and profundal sites. Poor chironomid communities from the deep profundal zone would also hamper palaeoclimate reconstruction. Cold-stenotherm chironomid taxa were influenced by river proximity, while taxa preferring warm conditions were more frequent at shallow littorals of the lake.

Weak negative correlation between mean grain size and water depth is explained by end-members revealing influences of river input and diatom valves in the grain-size composition.

Observed TOC values (mean 4.9 %) and TOC∕TN ratios indicate strong allochthonous supply of organic matter from the Utuk River. δ¹³C (mean −26.8 ‰) indicates dominance of C₃ plants and phytoplankton in the bulk organic matter fraction. Radiocarbon dating suggests that there is a reservoir effect caused by input of old organic carbon by a maximum of 720±30 ¹⁴C yr BP.

Elemental (XRF) data and mineral (XRD) distribution is influenced by the methamorphic lithology of the Stanovoy mountain range. Ratios of minerals relative to quartz decrease from the Utuk River towards the northern lake basin. Ti correlates well with mean grain size. There is no clear pattern in Mn∕Fe ratios, due to mixture of allochthonous elements and differential intensities of redox processes in the lake basin.

The observed proxy variabilities in the surface sediments suggest at least two locations for sediment coring: (1) at intermediate depth in the northern basin to account for representative bioindicator distributions and (2) the deep part in the central basin to potentially receive higher temporal resolution in the sedimentary record.