Late Quaternary Climate Reconstruction and Lead-Lag Relationships of Biotic and Sediment-Geochemical Indicators at Lake Bolshoe Toko, Siberia

Introduction

Over the past few decades, boreal and high elevation regions warmed faster than elsewhere (Mountain-Research-Initiative et al., 2015) and was strongest in northern Russia (Biskaborn B. K. et al., 2019). Although the Siberian region is generally regarded as being underrepresented in paleoenvironmental research (Sundqvist et al., 2014; Mckay et al., 2018; Kaufman et al., 2020), there is an ever-expanding interest in the area due to expected environmental changes and associated strong socio-environmental implications (Amap, 2017). In particular, remote Arctic lakes in permafrost regions represent pristine early-warning systems of environmental change, as they rapidly respond to climatic perturbations and changes in geo-ecological boundary conditions in their catchments (Subetto et al., 2017; Biskaborn et al., 2021; Nazarova L. B. et al., 2021). In terms of environmental sustainability, they still provide an infinite resource of uncontaminated fresh water, which make them unique compared to lake systems in the highly populated mid-latitude areas (Serreze et al., 2006). In order to assess the dimensions of modern landscape changes in Siberia, information about ecosystem changes during comparable climate warming episodes of the past provide valuable clues. In the terrestrial part of Siberia, such information, for instance, is stored in lake sedimentary records. The last strong warming of the atmosphere in a comparable setting took place after the Pleistocene-Holocene transition, referred to as the Holocene Thermal Maximum (HTM) that progressed in a non-uniform spatiotemporal event (Kaufman et al., 2004; Renssen et al., 2012). Knowledge about trends within this event was successively gained from paleolimnological records across a north-south transect through Siberia (Biskaborn et al., 2012; Biskaborn et al., 2016). These and our new record rely on radiocarbon dating and sediment-geochemical parameters that provide the background information on the sedimentation history necessary to reliably compile and evaluate paleoenvironmental bioindicator records (Biskaborn et al., 2013b; Wang et al., 2015; Heinecke et al., 2017).

The response of lake ecosystems to climate forcing is well represented by biotic indicators (Cohen, 2003). One of the main lake ecosystem components is given by diatoms, which are unicellular, siliceous microalgae that appear ubiquitous with opaline valves being well preserved in the sedimentary record allowing identification up to sub-species level by high-resolution light microscope analysis (Battarbee et al., 2001). Diatom assemblages (species communities) respond to hydrochemical changes (pH, Oxygen, nutrients) as well as seasonal climate changes via changes in the length of the ice cover, but also to physical habitat changes (water depth, sediment input, turbulence) (Biskaborn et al., 2012; Pestryakova et al., 2012; Biskaborn et al., 2013a; Herzschuh et al., 2013; Hoff et al., 2015; Biskaborn et al., 2016; Palagushkina et al., 2017; Pestryakova et al., 2018).

The simplest straightforward measure to describe diversity in lakes is alpha diversity (species richness). Alpha diversities in wetland ecosystem are affected by past long-term and ongoing recent climate change (Engels et al., 2020). In particular increasing summer temperature at northern sites, as a main driver, leads to increases in diatom diversity (Hobbs et al., 2010). As diatoms represent the main aquatic primary producers that are used as bioindicators in boreal lakes (Smol and Stoermer, 2010), the understanding of diatom diversity is important because it can stabilize primary production and thus overall biomass production (Marzetz et al., 2017). This in turn relevant for the global carbon budget as lakes act as sink for organic matter (Vyse et al., 2021 accepted and close to publication).

Larval head capsules of Chironomidae preserved in lake sediments have proven to be among the best indicators of environmental variations (Battarbee, 2000; Smol et al., 2005). Inference models (Barley et al., 2006; Nazarova et al., 2015; Luoto et al., 2016), which link the present distribution and abundance of chironomids to contemporary climate enable quantitative reconstructions of fundamental environmental parameters, among which are air paleotemperature (Nazarova et al., 2015), water depth (Nazarova et al., 2011a), continentality (Self et al., 2011), late-winter hypolimnetic oxygen (Luoto and Salonen, 2010).

Sedimentary pollen and their modern vegetation analogues have widely been used to reconstruct terrestrial summer temperatures (Biskaborn et al., 2016; Rudaya et al., 2016; Zhilich et al., 2017; Andreev et al., 2021). In our study design, we include a published pollen record from a parallel sediment core (Courtin et al., 2021) as additional independent proxy and counterpart to chironomid based air temperature reconstruction. Sediment-geochemical proxies that were analyzed from Bolshoe Toko surface sediments (XRF-gained elements, XRD-gained minerals, laser-derived grain-size distributions, organic carbon characteristics) are also used in support to biotic indicators. As novelty in context to previous research on Bolshoe Toko our study presents time series of aquatic ecosystem indicators and numerical paleoclimate reconstructions that can be correlated to sediment-geochemical downcore data. The resolution of bioindicator samples allows new insights into environmental variability and relationships up to centennial scale over the Holocene, while in the Pleistocene the age model (due to necessary projection of modern reservoir assumptions) only allows reliable conclusions at millennial scales.

The overall objectives of this paper are to 1) reconstruct millennial-scale climate change history in the yet understudied region of southern Yakutia and, by setting multi-proxy data in a paleoclimate context, to 2) infer relationships between diatom diversity, paleotemperature, and sediment-geochemical lake system trajectories.

Study Site

Lake Bolshoe Toko (56°15′N, 130°30′E, 903 m.a.s.l) was reported as possibly the deepest lake in Yakutia (Zhirkov et al., 2016). The freshwater lake is oligotrophic and circumneutral. Biskaborn B. K. et al. (2019) provided a detailed description of the lake system and reported pH 6.8 and conductivity 35.1 µS cm⁻¹ at the coring site PG2208 (Figure 1). According to their findings, the characteristics of the lake as coring locality for paleolimnological study can be summarized as follows: 1) pristine lake (far from settlements), 2) potentially old (multiple surrounding moraines and fault systems), deep (max. ca. 72.5 m), lake ecosystem represented in sediments (diatoms and chironomids preserved among multiple other proxies).

FIGURE 1. Study site of Lake Bolshoe Toko. (A) Overview map. (B) Lake Bolshoe Toko with coring sites. Water depth isobaths based on Konstantinov (2000) slightly modified from Biskaborn B. K. et al. (2019). Moraine cartography and topographic features based on Kornilov (1962).

Today, the sedimentary basin of Lake Bolshoe Toko is 15 km long and 8 km wide, situated at the northern foot of the Stanovoy Mountains representing most of the lake’s catchment (Biskaborn et al., 2021). The Bolshoe Toko basin forms a surficial landform on the southeastern Siberian Platform, which is built up of cratonic basement rocks of Precambrian age and overlain by Paleozoic to Mesozoic sedimentary rock strata and Quaternary cover sediments, in places with intercalated volcanic rock units (Nikishin et al., 2010; Barnet and Steiner, 2021). During late Jurassic to early Cretaceous time, the area formed part of an intramontane depression within the Aldan-Stanovoy Block (Imaeva et al., 2012) filled up with continental sedimentary rocks and coal deposits of the southern Yakutian coalfields (Golovin et al., 2007). Today a fault line separates the mountain range of the Stanovoy Mountains from Lake Bolshoe Toko, which is deepened into the coal deposits of southern Yakutia. Regional mountain glaciation during the Quaternary shaped the lake basin and left behind three generations of morainic arcs and ancient lake terraces of yet unknown ages (Kornilov, 1962; Konstantinov, 2000) (Figure 1).

Modern climate in the study area reveals strong seasonal contrasts under the regime of the Siberian High during winter and low-pressure conditions during the summer months (Mock et al., 1998; Shahgedanova, 2002). Regional climate data are available from the meterological station Toko, situated about 8 km north of Lake Bolshoe Toko (www.worldmeteo.info/en/europe/russia/toko/weather-101939). There, warmest temperatures and highest precipitation occur in July, ranging on average between 6 and 22°C and rain falls around 95 mm per month. Coldest temperatures between −49 and −29°C and low precipitation around 10 mm per month characterize the winter conditions in January. Climate conditions control the extent of lake ice. The study of modern satellite images (https://earthdata.nasa.gov) revealed that annual ice cover usually began in the first week of November and lasted until the end of May to mid-June during the last decade. In earlier times, lake ice cover may have existed longer, as for example observed in the year 1971, when lake ice did not disappear before the mid-July (Konstantinov, 2000).

The study site belongs to the Aldan Floristic Region (Kuznetsova et al., 2010) dominated by northern taiga consisting mainly of Larix accompanied by Picea and Pinus (Courtin et al., 2021).

Materials and Methods

Expedition and Core Retrieval

We conducted field work at Lake Bolshoe Toko in the first week of April 2013 as part of the expedition “Yakutia 2013” as a joint campaign of the North Eastern Federal State University (NEFU, Russia, Yakutsk) and the Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research (AWI, Germany, Potsdam). Coring equipment was brought by trucks onto the lake. We used a Jiffy ice auger (250 mm) to open holes through the 80 cm thick ice cover as conduits for scientific instruments and coring equipment to the unfrozen water column and bottom sediments. A portable Echo sounder (HONDEX PS-7 LCD) gave estimates of water depth and verified ambiguous values with a calibrated rope. We used UWITEC piston and gravity corer to retrieve sediment cores. Core PG2208 was retrieved by repeated coring of sub-sections at a water depth of 68.3 m. To avoid core loss due to the piston coring technique, parallel 3 m core sections were retrieved with an overlap of 0.5 m. Compaction was reduced by leaving the cores for tension release for several hours after retrieval. Sediment cores were transported in PVC tubes to the AWI laboratories in Potsdam and stored at 4°C in a dark room. Hydrochemical measurements performed during fieldwork are already published in detail in Biskaborn B. K. et al. (2019). The study in hand is based on core PG2208 and additionally involves pollen time series published by Courtin et al. (2021) from a core at 26 m water depth.

X-Ray Fluorescence Scanning

Core sections were opened alongside and cut into halves of which one was used for non-destructive XRF scanning and back-up material, while the other one was subsampled for single sample analyses. We used an AVAATECH XRF Core Scanner at AWI (Bremerhaven), equipped with Rh x-ray tube at 1.2 mA, 12 s count time for 10 kV (no filter) and 1,5 mA, 15 s count time for 30 kV (Pd filter) to scan the core in 10 mm steps for elements ranging from Al to Bi within a He atmosphere. The mean error (expressed as chi square χ²) for elements of interest was Al 2.0, Si 3.5, Ti 1.6, Mn 14.4, Fe 12.1, Br 0.6, Rb 0.7, Sr 1.0.

Core Splicing and Dating

To gain a composite core record for core PG2208 (Figure 1) we used XRF data, i.e., Sr, to correlate overlapping parts of retrieved sections resulting in 367 cm in total. For dating we scanned the sediment core for macro fossils and retrieved bulk sediment samples at depths without any reliable macro remains. A prior attempt to date samples at the Poznan radiocarbon laboratory resulted in a large deviation between humic acid fraction (soluble in NaOH) and the humin fraction (residual, i.e., alkali-insoluble fraction). Therefore, we kept only the Poznan date of the single macro fossil (terrestrial wood remain) that was present and analyzed with standard (reliable) dating routines. For creating a reliable model, we dated another set of bulk sediment samples at MICADAS radiocarbon laboratory at AWI using an acid treatment method outlined in Vyse et al. (2020). All dating results are shown in detail in Courtin et al. (2021) including the age-model for PG2133 used for pollen in their study.

To estimate a reservoir (old carbon) effect we dated an undisturbed surface sample PG2119 that was 1.6 km SW of PG2208. To estimate recent accumulation rates, we developed chronologies for a short core PG2203 that was 935 m SE from coring location PG2208 using the Lead-210 (²¹⁰Pb) and ¹³⁷Cs technique by direct gamma assay in the Liverpool University Environmental Radioactivity Laboratory, using Ortec HPGe GWL series well-type coaxial low background intrinsic germanium detectors (Appleby et al., 1986). Details on the methods including analyses of short cores from Bolshoe Toko are provided by Biskaborn et al. (2021). All MICADAS datings performed on core PG2208, the Poznan dating of the plant remain, and the Lead-210 samples are shown in Table 1. We applied the R package “Bacon” (Blaauw and Christen, 2011) and the IntCal20 calibration curve (Reimer et al., 2020) to model the age-depth relationship based on 11 radiocarbon dates from PG2208 and nine calibrated ages (Lead-210 and ¹³⁷Cs technique) from the short core PG2203 reported onto composite depths from PG2208 (Figure 2). We used a reservoir value of 762 ± 66 years based on PG2119 to correct bulk sample radiocarbon estimates (Table 1).

TABLE 1. List of samples and dating methods used for age-depth modelling. All dates are converted to C14 ages relative to 1950 CE. All dating results from MICADAS and Liverpool are shown. Just one Poznan date on wood residuals (Courtin et al., 2021) was used in this study. The composite depth indicated as “Com. Depth (cm)” represents the total depth below the sediment surface; the specific depth (cm) within core sections is indicated in “Core sample ID”.

FIGURE 2. Bacon age-depth model based on radiocarbon data from sediment core PG2208 (367 cm) and ²¹⁰Pb and ¹³⁷Cs dating from the uppermost 4 cm of short core PG2203.

Grain Size

Before conducting grain-size analysis, we treated 36 samples with hydrogen peroxide (30%), set on a platform shaker for 4–5 weeks to remove organic matter, while keeping the pH neutral using ammonia and acetic acid. After washing, sediment samples were treated with sodium pyrophosphate for dispersion and homogenization of the sediment suspension for a day. At least three measurements per sample were carried out using a laser diffraction particle size analyzer MASTERSIZER 3000. Their mean was used to estimate sand, silt, and clay fractions as sums of volume percentages between 2 and 63 mm, 63 and 2 mm, and <2 mm grain-size classes, respectively.

Carbon and Nitrogen

To gain information on organic matter variability we analyzed 204 samples for total carbon (TC) and nitrogen (N), 190 samples for total organic carbon (TOC), 45 samples for stable carbon isotopes (δ¹³C), and determined the water content in 417 samples. We used an Elementar Vario EL III to estimate TC and N from two repeated measurements of each freeze-dried and milled sample. The measurement accuracy is 0.05% for TC and 0.1% for TN. We used a Vario MAX C analyser to determine TOC. Measurement accuracy was 0.1%. We analyzed δ¹³C from carbonate-freed samples using a Finnigan Delta-S mass spectrometer. Results are presented as δ¹³C values relative to the PDB reference standard (‰) with an error of ±0.15%.

X-Ray Diffractometry

To analyze the mineralogical composition of sediments we freeze-dried and milled 19 samples for X-ray diffractometry (XRD) using a PHILIPS PW1820 goniometer (40 kV, 40 mA, from 3 to 100°, step-rate 0.05°, Co ka radiation). Data processing was performed using the free software MacDiff 4.0.7. To avoid influence of opal and organic matter in the samples (Voigt, 2009), we used the ratio of the peak area intensities of each mineral to the sum of the peak area intensities of all detected minerals, and used this value to estimate percentages.

Diatoms

A total of 37 samples from sediment core PG2208 were analyses for diatoms, of which we found valves in 32 samples. Average resolutions for diatoms and other indicators over the Holocene and Pleistocene was generally centennial and millennial, respectively (listed in the Supplementary Table; Holocene 264 years, Pleistocene ca. 1,400 years). For species analysis we prepared diatom slides for light microscopy following the common procedure developed by Battarbee et al. (2001). About 0.1 g of freeze-dried sample material was treated with hydrogen peroxide (30%) for several hours and washed before adding 5 × 10⁶ microspheres to estimate diatom valve concentration (DVC). Diatom valves were identified using a Zeiss AXIO Scope.A1 light microscope with a Plan-Apochromat 100x/1.4 Oil Ph3 objective at 1,000× magnification. A minimum of 300 diatom valves (Wolfe, 1997) were counted (mean 316) in each sample. Only 23 valves were found in the deepest sample at 310.5 cm. In parallel to diatom analysis, Chrysophyte cysts and Mallomonas scales were counted but not further identified. Diatom species identification was supported by images from a scanning electron microscope (SEM) at GeoForschungZentrum Potsdam. The diatom species identification was carried to lowest possible taxonomic level. To guarantee correct identification literature was used including Hofmann et al. (2011) and Krammer and Lange-Bertalot (1986–1991) as well as online databases (i.e., http://www.algaebase.org).

Chironomids

A total of 57 samples were studied for chironomid analysis. Sample resolution was 317 years. in the Holocene and ca. 1800 years in the Pleistocene (Supplementary Table). However, the lower (291–366 cm) and the upper (0–40 cm) parts of the core contained less than 50 chironomid head capsule remains per sample, that is necessary for a reliable quantitative paleotemperature reconstruction (Heiri and Lotter, 2001; Quinlan and Smol, 2001). In order to reach the adequate number of head capsules for a reliable estimate of inferred temperature we merged adjacent samples from these two core intervals. Treatment of 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 10 min by adding boiling water and left for 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 and were mounted in Hydromatrix two at a time, ventral side up, under a 6 mm diameter cover slip. 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 was taken from Brooks et al. (2007), Nazarova et al. (2008), Nazarova et al. (2011a), Moller Pillot (2013), Nazarova et al. (2015) and Nazarova L. B. et al. (2017).

Data Processing and Statistical Techniques

For statistical analysis of diatoms and sediment-geochemical variables we used the R environment (R Core Team, 2016).

XRF data was cleaned, i.e., negative values (<1%) resulting from the model were excluded. Therefore, some graphs reveal minor data gaps. For better comparison with lower resolution proxies and to reduce noise, we applied 7-point running mean shown on top of unsmoothed data. We transformed the data to account for compositional data using centered log ratios (CLR) for single elements and additive log ratios (ALR) for ratios of two elements.

Diatom zones were established guided by a constrained incremental sums-of-squares clustering (CONISS) applied on all species counts while using Euclidean dissimilarity after log transformation to downweigh abundant species (Grimm, 1987). Before choosing a suitable ordination technique for diatom data we performed a detrended correspondence analysis (DCA) to calculate gradient length (in standard deviation SD units) using the “decorana” function in the package “vegan” (Oksanen et al., 2020). The DCA gradient length for the original count data was (SD = 3.1) which is above the threshold suggested by Birks (2010). Therefore, to reveal major downcore trends in diatom assemblage data, multidimensional scaling (MDS) technique based on Euclidean distance was chosen (Mardia, 1978) using the “cmdscale” function from the package “stats”. Before MDS analysis we filtered the diatom dataset to include only species present with ≥2.5% in ≥2 samples, excluded the bottom sample at 310 cm because of too scarce valve counts, and applied square-root transformation to downweigh abundant species (but less aggressively than log transformation on unfiltered data applied for CONISS).

We estimated diatom alpha diversity (species richness) based on Hill’s N0, and N2 diversity estimates (Hill, 1973) using all diatom counts excluding the bottom sample that had too low number of counts. We applied rarefaction to correct diversity estimates for sampling differences (Birks et al., 2016) using the minimum base sum of all samples (n = 228). Hill’s diversity number N0 estimates species richness, while N2 estimates the inverse of the Simpson concentration that amplifies abundant species but discounts rare species (Chao et al., 2014). Rarefaction and diversity indices were calculated using the “vegan” R package (Oksanen et al., 2020).

The quantitative chironomid-based reconstruction of mean July air temperatures (T_July) is based on calibration chironomid data sets of the lakes from northern Russia (Nazarova et al., 2015). T_July was inferred by using a North Russian (NR) chironomid-based temperature inference model (WA-PLS, 2 component; r2 boot = 0.81; RMSEP boot = 1.43°C) based on a modern calibration data set of 193 lakes and 162 taxa from East and West Siberia (61–75°N, 50–140°E, T_July range 1.8–18.8°C) (Nazarova et al., 2015). Modern mean July air temperature for the lakes from the calibration data set was derived from New et al. (2002). The T_July NR model was previously applied for paleoclimatic inferences in Europe, Eurasia with emphasis on Arctic Russia (Nazarova et al., 2013; Syrykh et al., 2017; Wetterich et al., 2018; Plikk et al., 2019; Nazarova L. et al., 2021; Nazarova L. B. et al., 2021), and demonstrated high reliability of the reconstructed parameters. The chironomid-inferred T_July 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). Chironomid-based reconstructions were performed and results were displayed in stratigraphic diagrams produced using the software C2 version 1.7.7 (Juggins, 2007). The chironomid data were square-root transformed to stabilize species variance. Zonation of stratigraphy was accomplished using the optimal sum-of-squares partitioning method (Birks and Gordon, 1985) and the program ZONE (Lotter and Juggins, 1991). The number of significant zones was assessed by a broken stick model (Bennett, 1996) using program BSTICK (Birks H.J.B. and Line J.M., unpublished).

We used the pollen dataset from Courtin et al. (2021) to perform July temperature reconstruction. Climate reconstructions are based on harmonized pollen data in accordance with Cao et al. (2013) i.e., woody taxa and common herb/shrub taxa are harmonized on genus level and all other taxa on family level; only terrestrial taxa were included. To minimize the impact of false-positive analogues (Cao et al., 2017), the used modern pollen training dataset was selected from sites within a 2000 km radius around the lake Bolshoe Toko from the Eurasian (Davis et al., 2020), a North American (Whitmore et al., 2005) and East Asian (Cao et al., 2014) modern pollen compilations. T_July and Pann values were extracted for 1701 modern pollen sites from WorldClim 2 data (Fick and Hijmans, 2017; https://www.worldclim.org/) having a T_July range from 3.1 to 24.3°C. The climate reconstructions were performed applying by using the MAT-function in the rioja package version 0.9-21 (Juggins, 2014) for R software (RCore Team, 2016). Cross-validation of modern data set yielded a root mean squared-error of prediction of 1.88°C and a r2 of 0.77 for observed vs. predicted T_July values. In addition, a statistical significance test (Telford and Birks, 2011) was performed for the reconstruction using the randomTF-function in the palaeoSig package (version 2.0-3).

To investigate lead-lag correlations we developed a moving temporal offset approach in the R environment (R Core Team, 2016). To make data sets comparable, new data points were predicted by applying local polynomial regression fitting using the “loess” function in the standard R package “stats”, based on Cleveland and Devlin (1988), and then linear interpolated the newly generated points to a 10-year resolution. This procedure helped us to minimize fluctuations and thus reduced noise in the integrated data set. It further ensured that data of different time series can be compared on a 10-year resolution for correlation.

We calculated the correlation coefficient for each 10-year interval step with a maximum offset of ±1,000 year while keeping one parameter fixed (either Hill’s N0 or N2 diatom diversity). Thus, 200 correlations were calculated for each involved parameter with the selected diatom alpha diversity, using the “cor” function of the package “compositions”, which is based on the Spearman correlation (Van Den Boogaart and Tolosana-Delgado, 2008). Parameters with a correlation coefficient above 0.7 or below −0.7 at any point are included.

Supportive and additional statistical analyses are provided in the Supplementary Material to this paper. To visualize multiproxy data structures, we performed a MDS and a Pearson correlation matrix based on the ensemble of indicator data (interpolated as described above for lead-lag correlation) presented in the study. These include chironomid- and pollen-based temperature reconstructions, sediment-geochemical variables (elements, grain-size, organic carbon estimates), and diatom-inferred variables (Supplementary Figures S1, S2). To compare our lead-lag analysis with other methods we also tested for cross-correlation. We used the same data pretreatment as described above but used 100-year interval steps instead. This change in intervals was necessary to make the data suitable for the correlation calculation in the “ccf” R function, based on Brockwell et al. (2016), for estimating the lead-lag relationships of two selected parameters. We used the outcome of our moving temporal offset approach to test the highest/lowest correlations of interest within the “ccf” function (Supplementary Figures S3–S5). 95% confidence intervals were used from the error margin of the standard error of the “loess” function.

We summarized the basic statistics of parameters used for this study in an overview, which is accessible in the Supplementary Material.

Results

Chronology

Age-depth modelling performed using the R package “Bacon” revealed a continuous age-depth relationship without any reversals (Figure 2). Inferred sedimentation rates from the applied age model show that accumulation during the Holocene (0.02 cm a⁻¹) was higher than in the Pleistocene (0.01 cm a⁻¹). Mean 95% confidence of age estimates over the Holocene were 842 years (±418 years), and 2,595 years (±1,298 years) over the Pleistocene core section. The maximum confidence range was 4,220 years (±2,110 years) at the bottom of the core (367 cm). The full set (100% of the dates) overlap with 95% envelope ranges developed by the Bacon model.

Sedimentology and Biogeochemistry

The Pleistocene part of the sediment core starting at 367 cm (ca. 23.3 cal ka BP), is characterized by medium grey, homogeneous, silty clay. In a lithological transition zone, it passes over to light brown finely laminated silty clay between 288 and 268 cm including the Pleistocene-Holocene boundary. The Holocene part, and especially the upper part starting at 200 cm (7.1 cal ka BP) to the top of the core, consists of homogenous but (partly visible) finely laminated brown organic mud (Figure 3).

FIGURE 3. Variability of grain size and organic geochemistry (TOC, total organic carbon; TC, total carbon; TN, total nitrogen; ∂13C, stable carbon isotopes) in Bolshoe Toko (sediment core PG2208). A stretched core photograph is added below the mean grain-size data at the left side of the figure showing the lithological development and fine lamination.

The generally silty grain-size composition revealed increasing but overall low mean grain size between the bottom and 268 cm (min. 2.4, max. 10.3, mean 4.9 µm). In this section sand is almost absent (mean 2.0 vol%) but clay revealed highest values (mean 31.2 vol%). In samples above 268 cm mean grain size revealed higher, but also strongly fluctuating values (min. 7.9, max. 38.7, mean 17.8 µm) and highest proportions of sand (mean 18.7 vol%), while clay decreased to 9.5 vol%. Mean proportion of silt over the entire core was 70.3 vol%.

The mean difference between TC and TOC is 0.3 wt%, with an exception of two samples at 0.25 cm (4.2 wt%) and 202.5 cm (3.7 wt%), which, regarding the high number of samples measured in both devices (n = 204), is assumed to represent outliers. In the Pleistocene, carbon is almost absent in the core. From the core bottom upward to ca. 288 cm mean carbon values remain below 1%. The transition zone between ca. 288 and 268 cm reveals rising carbon concentrations up to ca. 5 wt%. Carbon concentrations in the above Early Holocene between 268 and 200 cm are slightly but constantly lower (mean 5.5 wt%) than when compared to the mid- and Late Holocene part above 200 cm (mean 7.5 wt%). Fluctuations visible in (Figure 3) are likely an expression of the fine laminations. TOC/TN_atomic ratios generally follow the carbon concentrations increasing from ca. 6 to 16 by 268 cm and subsequently remaining continuously high between 16 and 20 until the core top.

The water content throughout the core increases slowly with lower values between the bottom and 268 cm (min. 30.1 wt%, mean 52.7 wt%) and higher values starting from 268 cm until the top (max. 90.5 wt%, mean 79.9 wt%).

δ¹³C values range from −27.4 to −26.7‰ (mean 27.0‰). While values decline in the Late Holocene, there is an increase to maximum values in the uppermost part of the core which corroborates findings from a nearby short core (PG2203, Figure 1), reported by Biskaborn et al. (2021).

Elements and Minerals

CLR and ALR transformed XRF data suggest a general shift in the sediment core at 280 cm (ca. 14.0 cal ka BP) (Figure 4). At that depth Sr/Rb ratios, Br, Mn/Fe, and Si/Al strongly increase, while Ti decrease and stay approximately uniform with minor but high-frequent fluctuations in the Holocene. The transition is rather abrupt with the exception of Br which starts to increase gradually already at 320 cm. The transition at 280 cm is in accordance with XRD inferred mineral composition which shows a general decline of plagioclase, K-feldspar, and hornblende at that depth, while quartz rises. The maximum (measured) quartz peak was found even earlier at 300 cm.

FIGURE 4. Sediment-geochemical downcore variability in Bolshoe Toko (sediment core PG2208). Elemental data from XRF scanning given as additive log ratios (ALR) or centered log ratios (CLR) with 5-point running averages (black lines). Mineral data from XRD given as semi-quantitative percentage data.

Diatom Species Assemblages

Diatoms occurred between 311 and 0 cm depth (ca. 18.0 and −0.063 cal ka BP) in the sediment core PG2208 (Figure 5). In total, 158 diatom species were identified, dominated by a variety of Cyclotella, Aulacoseira and small achnanthoid taxa. In addition, Pliocaenicus bolshetokoensis appeared in considerable quantities between 5.1 and 34.4% (mean 15.0%) in the entire core record. Diatom valve concentration varied between 1.1 and 9.2 10⁸ g⁻¹ (mean 3.4 10⁸ g⁻¹), started to increase at 200 cm (7.1 cal ka BP) towards the top of the core, revealed a low value at 101.5 cm and peaked at 35.25 cm. The valve preservation index (F-Index) was generally high, ranging between 0.79 and 0.94 in all samples except the lowermost two samples at 300 and 310 cm, that had 0.51 and 0.09, respectively. The lower part of the core between 290 and 311 cm was dominated by Cyclotella sp. C. intermedia and C. ocellata species reaching up to 47.8, 21.7, and 36.3% in some of the samples, respectively. Aulacoseira ambigua started to appear in high frequency between 290 and 250 cm (35.8–10.5%), remained at intermediate values (mean 9.0%) between 250 and 100 cm, and decreased until the top (mean 2.9%). A. subarctica had intermediate values (mean 7.1%) but peaked several times reaching maximum abundance of 23.5% at 250 cm, while A. distans and A. valida had rather minor occurrences with highest values between 290 and 240 cm. C. atomus and C. tripartita started to appear with minor occurrences in the lower part of the core but appeared in higher frequencies above ca. 200 cm, reaching 23.3 and 18.4%, respectively. C. comensis (max. 23.6%) started to appear in high frequencies at 200 cm, and similar to C. cyclopuncta, increased distinctly again between 70 cm and the top. Achnanthidium minutissimum was abundant throughout the record with medium values between 1.9 and 18.5% (mean 10.0%). Other noteworthy taxa occurred throughout the record including Tabellaria (0–10, mean 2.8%), Fragilaria gracilis (0–5.5, mean 2.0%), and Rossithidium pussilum (0–6.2, mean 1.4%), of which the latter decreased in the uppermost part of the core. Mallomonas scales appeared from 270 cm onward and peaked at 170.5, 110.5, 80.5, and 12.5 cm before the Mallomonas index (0–38.3, mean 12.8) reached a constant maximum in the uppermost 3.25 cm of the core. The chrysophyte index peaked between 270 and 210, 140 and 80 cm, and revealed minor peaks in the uppermost core section.

FIGURE 5. Relative abundance of diatom species, diatom valve concentration, Mallomonas and Crysophyte indices, and sums of selected planktonic Aulacoseira and Cyclotella (sensu lato) taxa. Zonation guided by CONISS. Values are plotted next to composite core depth and median age values from the Bacon model. The species shown represent taxa present with ≥2.5% in ≥2 samples including the bottom sample.

MDS analysis depicts grouping of key-planktonic and benthic species along the primary axis and in particular a distinct clustering Cyclotella versus Aulacoseira together with Pliocaenicus and Achnanthidium minutissimum along the secondary axis (Figure 6). Both axes together explain 76.1% of the data variance.

FIGURE 6. Biplot of the first two dimensions (Dim1, Dim2) generated by multi-dimensional scaling (MDS) analysis of diatom species data. Colored circles indicate eco-taxonomical clusters with similar environmental preferences. Explained variance of each dimension is given at the axis label.

4.5 Chironomid Species Assemblages

In total, 65 chironomid taxa were identified. In the lowermost samples (below 290 cm; 15.7 cal ka BP, Figure 7) we found low concentration of chironomid remains represented only by profundal taxa from the genera Heterotrissocladius (H. subpilosus-type, H. grimschawi-type, H marcidus-type, H. maeaeri-type 2). Diversity of chironomid communities remains low: Hill`s N2 varies between 1 and 3.7. Reconstructed T_July range from 7.1 to 7.8°C.

FIGURE 7. Relative proportions of the most abundant chironomid taxa (%), concentration of chironomid head capsules (HC), PCA axes 1 scores for chironomid data, chironomid-inferred mean July air temperature (T_July), Hill’s N2 diversity and relative abundances (%) of littoral, profundal, lotic and semiterrestrial taxa in the sediment core PG2208 of Lake Bolshoe Toko. Grey vertical line at T_July represents modern mean July air T (°C).

From the core depth of 280 cm (14.0 cal ka BP) upwards, the diversity and concentration of head capsules rise. Between 280 and 200 cm (14.0–7.1 cal ka BP) though still at low abundancies, several littoral taxa (Cladotanytarsus mancus-type, Corynoneura-type, Psectrocladius taxa, Orthocladius type S) and semiterrestrial Smittia-Parasmittia appear in the lake. Reconstructed T_July rise from 7.8°C at the core depth 291 cm (15.85 cal ka BP) to 10.0°C at the depth 211 cm (7.6 cal ka BP). Characteristic taxa for flowing waters (lotic) appear after the core depth of 231 cm (8.8 cal ka BP).

The strongest taxonomic shift in the chironomid assemblages, reflected as well by the decline of the PCA 1 below 0, was found at 201 cm depth (7.1 cal ka BP; Figure 7). At this depth we observe rise of species representing littoral, lotic and semiterrestrial taxa. Taxonomic diversity of chironomid communities (Hill’s N2) increases up to 14.6 at the core depth 121 cm (3.9 cal ka BP). Among the profundal taxa we observe a shift towards dominance of acidophobic Micropsectra insignilobus-type. Abundances of semiterrestrial and lotic taxa are the highest in the record and phytophilic taxa from Cricotopus and Orthocladius genera became more abundant. Reconstructed T_July between 201 (7.1 cal ka BP) and 111 cm (3.4 cal ka BP) of the core depth vary between 9 and 13°C.

At a depth of 161 cm (5.6 cal ka BP) we observe a sharp decrease of the head capsules concentration, decrease of the chironomid diversity and a strong increase in the abundance of the semiterrestrial taxon Limnophies-Paralimnophies.

From 151 cm (5.2 cal ka BP) to 111 cm (3.4 cal ka BP) the acidophobic Micropsectra insignilobus-type and lotic and semi-terrestrial fauna declined while profundal acidophilic Heterotrissocladius taxa increase again in line with the rise of Corynoneura arctica-type and several phytophilic taxa from the genera Cricotopus, Orthocladius, Paratanytarsus, Cladotanytarsus. The reconstructed T_July remain slightly above the modern temperature (10.5–11.9°C).

After 111 cm there is a break in the development of chironomid communities. No chironomids have been found till 86 cm, when the chironomids re-appear in the record. The proportion of littoral taxa rises up reaching 66% at 86 cm (2.4 cal ka BP) alongside with the latter increasing abundances of lotic fauna that reach 22% at the depth 71 cm (1.9 cal ka BP) which could indicate lake deepening and new expansion of littoral areas. At the 71 cm Hydrobaenus lugubris-type appears. Reconstructed T_July decreases from 10.5°C at 111 cm to 8.1°C at 71 cm.

Above 71 cm the concentrations of lotic taxa decrease and above 9 cm (ca. 188 cal years BP) they are not found in the chironomid communities anymore. The relative proportion of profundal to littoral taxa remain at ca. 2.2 on average until ca 9 cm. Between 9 and 7.5 cm (ca. 147 cal years BP) the concentration of littoral taxa starts to increase and the concentration of the profundal taxa decreases again. The communities are dominated by profundal acidophobic M. insignilobus-type. Reconstructed T_July has an increasing trend towards 15 cm (11.7°C; 394 cal years BP) then a decrease with a min. of 9.4°C at ca 8 cm and a rise again towards in modern time up to 11.2°C.

Discussion

Our following discussion will start with inferring the paleoclimate history from Lake Bolshoe Toko data and synthesize the findings with adjacent paleoenvironmental records, i.e., along a north-south profile through Yakutia. Subsequently we discuss the ecological developments of the lake and inter-proxy relationships, such as lead-lag correlations between temperature reconstruction, bioindicator indices, and sediment-geochemical trends.

Paleoclimate Variability in Northeastern Siberia

Evaluation of Climate Reconstructions From Bolshoe Toko Sediment Cores

Here we discuss pollen-based and chironomid-based T_July reconstructions performed on samples from two parallel sediment cores from Lake Bolshoe Toko. Chironomid-based reconstruction in core PG2208 provides information derived from aquatic bioindicators; pollen-based reconstruction taken from Courtin et al. (2021, core PG2133) provides terrestrial-based information from the catchment including aeolian pollen transport (Figure 1). Deviations between both reconstructions (Figure 8) can be related to different conceptual proxy and ecosystem settings. Distributions of chironomids are mainly driven by air temperature (Brooks et al., 2007; Self et al., 2011; Nazarova et al., 2015) but can be influenced by changes in limnic habitat preferences (e.g., water depth, substrate, ice cover, oxygen) (Stief et al., 2005; Nazarova et al., 2011b; Engels et al., 2012), while terrestrial vegetation compositions are limited by water supply (e.g., Antonsson et al., 2006). A major palynological constraint is that long-range transport of pollen obscures vegetation-based temperature reconstruction at the single lake site, that is especially true in mountain regions (Rosén et al., 2003) and amplifying uncertainties caused by complex vegetation-landscape interactions and long-term lagged climate response, i.e., in the Holocene (Cao et al., 2019). Therefore, in our study we use T_July (pollen) mainly for inferring climate trends at the Pleistocene/Holocene boundary where pollen shifts are well established (Velichko et al., 2002), and T_July (chironomid) within the Holocene interglacial, when limnoecological conditions were relatively stable. Potential deviations caused by age control or different sampling strategies cannot be ruled out, but are likely neglectable as common features for both reconstructions. One major feature for both data sets is a temperature rise at the Pleistocene-Holocene boundary. Here, pollen data reveal strong short-term fluctuations. Artemisia and Poaceae increased while Betula and Alnus decreased around 12.5 cal ka BP, suggesting evidence for a Younger Dryas cooling (Courtin et al., 2021). Climate cooling and increased continentality during the Younger dryas was reported for other sites in southern (Harding et al., 2020) and eastern Siberia (Velichko et al., 2002; Müller et al., 2009). However, tree-growth (i.e., Larix) persisted during this cold phase at multiple sites in eastern Siberia (Müller et al., 2009; Werner et al., 2010).

FIGURE 8. Graph showing synchronized variables used for lead-lag correlations. Temperature reconstructions (T_July, °C) by pollen are from core PG2133. All other variables are from core PG2208: chironomids, diatom alpha diversity estimates (Hill’s N2, N0), MDS first dimension Dim 1 (multiplied by −1), organic carbon concentration, selected elements (CLR values), and total sand concentrations. Variables are drawn using their ages from PG2208 and PG2133 for pollen (Courtin et al., 2021). Smoothing and confidence intervals were gained using LOESS. The less reliable (lowermost four) samples were excluded from the chironomid data due to the low number of head capsules. HTM, Holocene thermal maximum; P/H, Pleistocene/Holocene transition.

Our pollen-based reconstruction of summer temperatures from Bolshoe Toko suggests that T_July (pollen) fluctuated around 11°C during the late Glacial with slightly higher values during 13 to 12 cal ka BP and a reversal shortly before the early Holocene warming. Holocene temperature ranges around 13°C with a slight decline toward the late Holocene. The T_July reconstruction was significant (p < 0.028) (according to Telford and Birks, 2011). In detail, pollen reconstructions revealed a distinct but ambiguous climatic event between ca. 12.9 and 11.7 cal ka BP (Figure 8) indicating cold temperatures at the beginning and the end of the Younger Dryas period, but a warm interim period (Figure 8). Such a short warming event in the middle Younger Dryas in Siberian records was already considered by Velichko et al. (2002). The lower resolution of our T_July (chironomid) reconstruction in this period can somewhat validate Younger Dryas cooling but not a middle Younger Dryas short warming episode.

The chironomid taxa from the genera Heterotrissocladius that are dominant until 13.0 cal ka BP are profundal, oligotrophic, acidophilic, cold stenotherm (H. subpilosus-type, H. grimschawi-type) to cold tolerant (H. marcidus-type, H. maeaeri-type 2) (Brooks et al., 2007). Such composition of chironomid communities indicates cold and dystrophic lake conditions unfavorable for development of rich chironomid fauna. Head capsule concentration was too low to enable high resolution quantitative temperature reconstructions from samples deeper than 290 cm (15.8 cal ka BP). Nevertheless, the reconstructed T_July (chironomid) was ca 2–3°C below modern.

After ca. 14.0 cal ka BP rising concentration of chironomid remains, diversity of chironomid communities and emergence of littoral and semiterrestrial taxa indicate beginning of development of suitable littoral habitats at the lake shores and slowly rising temperatures (Figure 8). Emergence of the lotic taxa after 8.8 cal ka BP reflects strengthening of the surface water inflow into the lake that can be caused by warming and eventually thawing of cryotic (mountain glacier) catchment features leading to rise of the lake water level.

The strongest change in the chironomid assemblage occurred at 7.1 cal ka BP and is associated with a replacement of acidophilic profundal cold-stenotherm Heterotrissocladius taxa by acidophobic Micropsectra insignilobus-type and a number of littoral, semiterrestrial, phytophilic and lotic taxa. The dominating acidophobic Micropsectra insignilobus-type, though mostly attributed to profundal conditions, is known for thriving as well in the littoral zone of mountain lakes (Lindegaard, 1992; Bitušık and Kubovcık, 1999). We assume that water depth remained high with well-developed littoral zones overgrown by macrophytes. This major shift supports chironomid-reconstructed warm temperatures that reached about 2–3°C above todays values (Figure 8) indicating a mid-Holocene Thermal Maximum (HTM) starting at 7.1 cal ka BP.

However, sharp decrease of the head capsule concentration, decrease of chironomid diversity, and rising abundance of semiterrestrial Limnophies-Paralimnophies possibly indicate unstable lake levels (Massaferro and Brooks, 2002; Brooks et al., 2007) at around 5,600 cal years BP.

Between 5.2 and 2.0 cal ka BP Micropsectra insignilobus-type and lotic and semi-terrestrial fauna declined while cold stenotherm Heterotrissocladius taxa, Corynoneura arctica-type and several phytophilic taxa increased. Corynoneuraarctica-type is littoral, possibly preferring floating leaved macrophytes, especially duckweeds and can be acidophilic because species of this genus were found in lakes with a pH around 4 (Moller Pillot and Buskens, 1990). This taxonomic shift indicates Neoglacial cooling climate conditions at Bolshoe Toko and lowering of the water pH. The pH decrease can be caused by respiration and decomposition processes in shrinking and shallowing littoral zones (Nazarova L. et al., 2017; Nazarova et al., 2020).

Between ca. 3.4 and 2.4 cal ka BP, the temperatures dropped. Chironomids were not found in the sediments during this period, probably as a result of lake level and productivity lowering in response to Neoglacial cooling. This cooling trend was also observed in Kamchatka (Brooks et al., 2015) consistent with the strengthening of the Siberian High (Hammarlund et al., 2015). Furthermore it is prominent in East Eurasian records as it was amplified by thermohaline circulation (Wanner et al., 2011).

At ca. 2.3 cal ka BP the rise of littoral taxa and recovery of Micropsectra insignilobus-type with increasing abundancies of lotic fauna indicate increasing water flow peaking at ca. 1.9 cal ka BP under a relatively cool climatic condition (chironomid T_July drops up to 2°C below modern). Here, emergence of Hydrobaenus lugubris-type confirms strengthening of the water flow. This taxon is specially adopted to the floodplains of large rivers (Moller Pillot, 2013) and might have been brought to the lake with inflowing rivers.

During the last ca. 2000 years the relatively stable ratio of profundal to littoral taxa indicates stable water levels during this interval. Smaller shifts between habitat preferences after ca. 188 cal BP could potentially be associated to water level changes related to centennial/decadal climate disturbance during the little ice age. However, the low number of chironomid head capsules in core PG2208 and the mean Holocene resolution of pollen in core PG2133 of 750 years did not allow sufficient resolution of the reconstructions of ecological conditions during the last millennium to reach centennial scale.

Bolshoe Toko as the Southern End of a North-South Paleoclimate Transect in Siberia

Our study completes the southern end of a N-S transect of paleoclimate records in Yakutia (Figure 9) established by multiple studies along the Lena River over the last decades (Laing et al., 1999; Andreev et al., 2004; Fradkina et al., 2005; Popp, 2006; Müller et al., 2009; Biskaborn et al., 2012; Klemm et al., 2013; Nazarova et al., 2013; Tarasov et al., 2013; Biskaborn et al., 2016; Kostrova et al., 2021). The Holocene warming trend in Siberia is a common and very clear feature reflecting the glacial-interglacial transition in the Northern Hemisphere (Kaufman et al., 2020). However, previous investigations of lake records along the transect showed that the onset of most pronounced warming within the interglacial lagged by up to ∼3 kyrs in southern direction (Biskaborn et al., 2012), contradicting the assumption that warming was primarily related to solar insolation (Bond et al., 2001; Miller et al., 2010). The global compilation of Holocene temperature reconstructions by Marcott et al. (2013) suggested that the Early to mid-Holocene warming was modulated by deglaciation of Northern Hemisphere ice sheets. Biskaborn et al. (2012) found that the spatiotemporal pattern along a latitudinal gradient in eastern Siberia matches with global climate models (Renssen et al., 2009; Renssen et al., 2012) reconstructing deglaciation processes of the Laurentide ice-sheets and associated supply of cold air to Central Siberia until ∼7 cal ka BP. An exception to this latitudinal trend in eastern Siberia is provided by the Lake Emanda record, where diatom opal ∂18O and species assemblage changes suggested slightly earlier warming at 7.9 cal ka BP compared to sites along the Lena transect (Kostrova et al., 2021). An explanation might be the location of Lake Emanda to the east beyond the Verkhoyansk Mountains. As a consequence, the mountain range could have acted as a climate barrier for cold air supply by westerlies.

FIGURE 9. Spatiotemporal variability of the Holocene thermal Maximum (HTM) in Siberia compiled from paleolimnological studies (Laing et al., 1999; Andreev et al., 2004; Fradkina et al., 2005; Popp, 2006; Müller et al., 2009; Biskaborn et al., 2012; Klemm et al., 2013; Nazarova et al., 2013; Tarasov et al., 2013; Biskaborn et al., 2016; Kostrova et al., 2021). Horizontal dashed orange lines indicate the onset of the HTM in each record. This graph was adopted from Biskaborn et al. (2016) and further developed in eastern and southern directions with two new sites.

The new Lake Bolshoe Toko record presented here demonstrates that the cooling mechanism in the Early Holocene could possibly also be applicable for terrestrial aquatic settings close to the Sea of Okhotsk, as Bolshoe Toko is ca. 300 km far from its coast inland. Renssen’s model show that there is a tendency of delayed warming towards latitudes near Bolshoe Toko (Renssen et al., 2009). However, the relative proximity of our study site to the marine coast possibly exceeds the resolution of their model.

On the other hand, pollen and chironomid inferred summer temperatures differ (Figure 8) suggesting, apart from the mechanisms explained above, a bias caused by comparisons of summer versus winter sensitive proxies as highlighted by Baker et al. (2017). Their study demonstrates that winter climate dominated Holocene temperature dynamics is associated with decay of ice sheets in the Northern Hemisphere in continental Eurasia in contrast to regions more proximal to marine coasts.

Several short-term cold events disrupting the overall Holocene climate development were reported by Wanner et al. (2011), which, in most cases, are not found in the Bolshoe Toko record, possibly due to low temporal resolution of the study. A prominent cooling signal that only was found in the chironomid-based T_July reconstruction was dated to 5.6 cal ka BP (5.1–6.2 ka 2σ range) but could not be correlated to common cold events. There is, however, a cooling signal found at around 2.5 cal ka BP in both chironomid- and, to a lesser degree, in pollen-based reconstructions (Figure 8). If taken into account the confidence range of the age model at this core section (2.2–2.7 ka 2σ range) this cold phase could match to cold events between 3.3 and 2.5 ka BP that were found in other Eurasian records associated to Bond events (Wanner et al., 2008).

Paleolimnological Variability and Inter-proxy Relationships

Development of the Sedimentary Depositional Environment

The sediment-geochemical variability in core PG2208 reveals a major change at the Pleistocene/Holocene boundary from siliciclastic erosional regime, represented by high proportions of minerogenic clastic input and low organic matter contents, to higher lake and catchment primary productivity. This is indicated by rising organic matter concentrations and plant remains present in the Holocene (Table 1; Figures 3, 4). This change is accompanied by an increase in sedimentation rates, possibly slightly overestimated because of the low consolidation of the pore-water-rich sediments.

Abiotic proxy data exhibit a decline of detrital elements and minerals (Ti, Feldspar, Hornblende) after ca. 15 cal ka BP (Figure 4). Mn/Fe as oxygenation proxy (Naeher et al., 2013) and Si/Al as biogenic opal and thus diatom proxy (Vyse et al., 2020) increased distinctly and more abrupt at ca. 15 cal ka BP, while production of organic matter (TOC, TOC/TN_atomic) increased more gradually between ca. 15 and 11 cal ka BP together with coarsening grain size. Changes in the mineralogy of detrital components, basically reflected by an increase in quartz at the expense of hornblende, point to modified routes of exogenic sediment supply. Sediments enriched in hornblende indicate sources from the crystalline rocks of the Stanovoy Ranges, while quartz-rich detritus typically points to shore erosion of older Quaternary sediments around the northern fringe of the lake (Biskaborn B. K. et al., 2019). This sedimentary pattern implies a successive delay of bio-related responses such as bioproductivity after physico-chemical changes already occurred. We therefore assume that the sediment-geochemical depositional setting of the lake system first experienced a rather qualitative change driven by provenance changes in a shrinking deglacial mountain catchment. Successive changes, i.e. decreases in winter ice-cover already led to stronger oxygenation and diatom growth, before increased summer warming initiated a fully established Holocene lake ecosystem. As a consequence, the beginning of the interglacial period was accompanied by a decrease of terrigenous sediment supply to the benefit of biogenic sedimentation.

Paleoecological Interpretation of Diatom Taxa

To understand the downcore diatom assemblage trends in a paleolimnological context we attribute all cyclotelloid and Aulacoseira taxa found in the core to principle ecological groups summarized in Figure 5. Lindavia comensis is described as opportunistic and fast-growing (Rühland et al., 2015), oligotrophic reference species that occurs in August in the Great Lakes (Wunsam et al., 1995) and is positively correlated to temperature and thermocline depth (Reavie et al., 2017). Cyclotella cf. atomus occurs in spring in the Great Lakes positively correlated to thermocline depth (Reavie et al., 2017), but not clearly with temperature. Lindavia cyclopuncta is reported as an oligotrophic species (Scussolini et al., 2011). Lindavia intermedia is reported as oligotrophic “bodanicoid” (Novis et al., 2017) and thus not necessarily related to thermal stratification (Rühland et al., 2015), growing in late summer and competing well under low light under P-limitation (Malik, 2016). Lindavia ocellata is described as planktonic and periphytic, ultraoligotrophic to mesotrophic taxon (Wunsam et al., 1995). Lindavia tripartite is correlated with temperature in the Great Lakes (Reavie et al., 2017). Cyclotella cf. iris was reported from fossil records (Khursevich et al., 2001) and thus we believe that a sparce occurrence in the upper part of the core could be misinterpreted. Other Cyclotella-like valves in the record could not be identified to the species level, just a few counts were similar to Cyclotella schumannii, also found in Lake Immandra, Kola Penninsula (Denisov and Genkal, 2018). Given the above assessment of ecological preferences we defined C. comensis, cyclopuncta, tripartite and cf. atomus as planktonic and lightly silicified forms (Figure 5) that dwell under thermally stratified conditions.

The genus Aulacoseira generally builds strongly silicified, heavy, and thus rapidly sinking frustules that are abundant in deep boreal lakes (Laing and Smol, 2003). Euplanktonic forms A. ambigua and A. subarctica, however, differ in their lifestyles to tychoplanktonic taxa A. distans and A. valida. We therefore interpret A. ambigua and A. subarctica as pelagic forms which require water turbulence to remain suspended in the photic zone (Gibson et al., 2003; Rühland et al., 2008) for comparison with abundancies of light cyclotelloid forms preferring stratified water.

A special feature of Lake Bolshoe Toko is the dominance of newly described Pliocaenicus bolshetokoensis (Genkal et al., 2018; Genkal and Gabyshev, 2018), which is, as an endemic planktonic species, ecologically not yet fully understood. Biskaborn B. K. et al. (2019) found that its spatial within-lake variability might relate this taxon to increased nutrient availability at the Utuk river mouth and other small inflow streams. However, this rather is contradicted by the general appearance of Pliocaenicus costatus (sensu lato) forms which are found in cold, oligotrophic, circumneutral to weakly acidic mountain lakes in East Siberia (Cremer and Van De Vijver, 2006).

Benthic species reveal only minor long-term trends, i.e., the most frequent benthic/tychoplanktonic form Achnanthidium minutissimum reveals relative uniform downcore abundances. The main ecological signal is visible in the decrease of Rossithidium pussillum in the last 1.5 cal ka BP parallel to the maximum of DVC, which, because R. pussillum is an oligotrophy indicator can be an indicator that DVC is driven by nutrient supply.

The MDS analysis revealed key-planktonic and benthic species along the first axis and in particular a distinct clustering of light planktonic Cyclotella species versus heavy and pelagic Aulacoseira frustules together with planktonic Pliocaenicus and the tychoplanktonic Achnanthidium minutissimum along the secondary axis (Figure 6). The data structure thus supports that changes in these key species, with ecologies explained above, are i.e., driven by stratification. Sample scores of the first MDS dimension (Dim1, Figures 5, 6, 8) therefore may indicate changes from turbulent to stratified system tendencies.

Late Quaternary Diatom-Assemblage Trends Driven by Paleoclimate Change

The earliest preserved diatoms are represented by only a few (23) valves of Cyclotella and Pliocaenicus species dated to ca. 17.8 cal ka BP. While C. ocellata was reported from a high diversity of cold lakes with differing preferences from ultra-oligotrophic (Cremer and Wagner, 2003) to eutrophic (Schlegel and Scheffler, 1999), C. intermedia is more indicative for low light and nutrient conditions and thus both forms are able to withstand harsh environmental conditions (Malik, 2016). Fragilarioid forms are typical pioneering forms in other boreal shallow lakes (Valiranta et al., 2011) but are missing in this central and deep part of the lake basin.

Valve concentration increased and preservation improved at ca. 15.7 cal ka BP when Aulacoseira taxa and also Achnanthidium joined the assemblage. At the same time July temperatures reconstructed from pollen (and less distinctly from chironomids) increased (Figure 8). The Pleistocene-Holocene boundary is represented by strong increase of euplanktonic Aulacoseira taxa (i.e., A. subarctica) and thus an increase in turbulent lake conditions related to water mixing caused by strong summer winds. Fluctuating and rather cool reconstructed July temperatures support this assumption.

The first strong increase of euplanktonic and light Cyclotella taxa is dated to ca. 7.1 cal ka BP indicating a distinct increase in thermal summer stratification (Rühland et al., 2015; Biskaborn et al., 2021). This shift is accompanied by increases in DVC as a warming indicator (Biskaborn et al., 2012) and in accordance with the onset of the Holocene Thermal Maximum (HTM) indicated by chironomid-inferred July temperature. Increasing temperatures are however not ratified by pollen-based reconstructions (see “T_July (pollen)”, Figure 8).

At ca. 3.4 cal ka BP low values of DVC, TOC, Cyclotella, Crysophytes and Mallomonas and a small clay peak point to a short-term event with low lake bioproductivity related to cooling inferred by chironomid data. The second major increase in light planktonic Cyclotella taxa, together with high DVC values, indicates thermal summer stratification starting at ca. 2.0 cal ka BP. Until recent decades, summers with stratification tend to be frequent. Although this would be contradictory to decreasing insolation over the Late Holocene, this development of the aquatic proxies (partially including chironomid-inferred T_July) can be explained by long-term trends of warming Siberian winters (Meyer et al., 2015) possibly leading to prolonged open water seasons. In this context, increases in the proportion of planktonic taxa after ca. 2.0 cal ka BP (Figure 5) can also indicate water depth increase as a result of increased snow melt in the mountain catchment leading to strengthened river water input. This could partly also have contributed to the formation of the eastern lake terraces (Figure 1) discussed in Biskaborn B. K. et al. (2019) after retreat of the lake level.

For the most recent history, Biskaborn et al. (2021) showed, that there is a distinct increase of summer stratification in Lake Bolshoe Toko during the industrial era as response to recent northern hemispheric warming, accompanied with acidification trends represented by Mallomonas.

Relationships Between Diatom Diversity, Paleoclimate and Sediment-Geochemistry

Reconstructed temperatures and sediment-geochemical parameters were correlated with diatom diversity (Hill’s N0 and N2, Figures 10A,B), revealing strong positive correlations to Hill’s N2 (Figure 10A, Supplementary Figures S1, S2). Reconstructed temperatures from both pollen and chironomids fit well, i.e., chironomid-inferred temperatures show strongest correlation with neglectable time offset. Direct linkage between diatoms and climate due to their fast reproduction rates (Julius and Theriot, 2010), and a general lag of terrestrial vegetation (and thus pollen-inferred temperatures) have been reported by Liao et al. (2020). The good fit between T_July (chironomids) and diatom diversity is supported by near synchronous compositional changes of both chironomids and diatoms over the Holocene reported by Solovieva et al. (2015).

FIGURE 10. Lead-lag correlations of diatom alpha diversity are based on (A) Hill’s N2 and (B) Hill’s N0, and include all involved proxies correlating >0.7 (solid lines indicate positive correlation) or <−0.7 (dashed lines indicate negative correlation). Before analysis the data were LOESS smoothed and re-sampled into an integrated decadal resolution time series. Then, spearman-based correlation was applied in a successive bi-directional 10-year moving temporal offset to maxima of ±1,000 year. Cross-correlation tests based on auto correlation function are provided in the supplementary material.

However, our analyses generally showed that limiting factors for detection of reliable lead-lag relationships is sampling strategy and resolution of the proxies involved. Therefore, to account for deviations in the proxy sampling strategies (Supplementary Material), we only infer directional lead or lag behaviors of proxies and omit illusive absolute temporal durations in this study.

XRF-derived elemental variability revealed positive correlations of redox-related elements (Fe, Mn) and Cl related to water content, and negative correlations of detrital elements (Al, Ti, K), all showing that these sediment-geochemical proxies shifted before diatom alpha diversity responded. Positive correlations to sand and negative correlations to silt and clay confirm this trend. Interestingly, negative correlations ahead of species richness appears slightly stronger for Hill’s N0 index, suggesting that sedimentary changes, i.e., the reduction of very fine (clay, Al) and coarse (gravel) components somewhat allows the initiation of rare taxa growth. The strongest positive correlations are found between measurements of organic proxies (TOC, TC, N) and Hill’s N2, and, to a lesser degree Hill’s N0. The observed centennial-scale delay of organic carbon behind changes in diatom alpha diversity demonstrates a strong causal relationship between species richness and production of organic matter (Marzetz et al., 2017). This finding may motivate further research on diversity-carbon trajectories as understanding of nature-based solutions for carbon sequestration will increasingly important for future climate change management.

Conclusion

The study of radiocarbon dated sediment cores from the remote Siberian Lake Bolshoe Toko enabled new paleoclimate reconstructions that were used to extend a north-south transect in Yakutia. The multiproxy data were further analyzed to gain understanding of the relationships between different paleolimnological system components, i.e., lead-lag relationships between diatom alpha diversity versus sediment-geochemical and paleoenvironmental variability. We summarize our conclusions as follows:

• The PG2208 sediment core revealed 23.3 cal ka BP over a length of 367 cm with higher sedimentation rates in the Holocene (0.02 cm a⁻¹) than in the Pleistocene section (0.01 cm a⁻¹).

• High resolution of pollen-based summer temperature reconstruction T_July (pollen) provides evidence for Younger Dryas cooling between ca. 12.9 and 11.7 cal ka BP and possibly an interim warm period.

• Chironomid-inferred summer temperature reconstruction T_July (chironomid) reveals an onset of Holocene Thermal Maximum (HTM) at 7.1 cal ka BP.

• Extension of a Siberian north-south paleoclimate transect shows that the delay of the HTM reaches far south close to the Sea of Okhotsk.

• Relative shifts in abundances of euplanktonic diatom species Aulacoseira and Cyclotella reveal repeated episodes of thermal stratification in response to mid-Holocene paleoclimate changes.

• Lead-lag correlations between diatom alpha diversity, involved environmental proxy data and temperature reconstructions revealed a lead of geochemical parameters, quasi-synchronous paleoclimate responses, and a lag of sedimentary organic carbon. Centennial-scale delay of organic carbon behind changes in diatom Hill’s N0 and N2 reveals the strongest causal relationship, which implies that species richness can augment the production and accumulation of organic matter in comparable lake systems. This may contribute to assessments of future carbon sequestration during ongoing global warming.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

BB led the sediment-geochemical and diatom analyses and writing of the manuscript. LN and LS led the chironomid analyses including chironomid-based reconstructions. UH led the pollen-based reconstruction. BD and LP led the field work and contributed to the environmental interpretation of the data. TK and GP assisted with data processing and statistical analyses. All authors contributed to analysis of the results and revision of the manuscript.

Funding

This research was supported by the BMBF project PALMOD #01LP1510D (Germany).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling Editor declared a past co-authorship with one of the authors (LN).

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Kim Funck for diatom identification and counting, and Vanessa Warmke for help in the laboratory for high resolution organic carbon and grain-size measurements.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.737353/full#supplementary-material

References