Effects of climate change and industrialization on Lake Bolshoe Toko, eastern Siberia

Industrialization in the Northern Hemisphere has led to warming and pollution of natural ecosystems. We used paleolimnological methods to explore whether recent climate change and/or pollution had affected a very remote lake ecosystem, i.e. one without nearby direct human influence. We compared sediment samples that date from before and after the onset of industrialization in the mid-nineteenth century, from four short cores taken at water depths between 12.1 and 68.3 m in Lake Bolshoe Toko, eastern Siberia. We analyzed diatom assemblage changes, including diversity estimates, in all four cores and geochemical changes (mercury, nitrogen, organic carbon) from one core taken at an intermediate water depth. Chronologies for two cores were established using 210Pb and 137Cs. Sedimentation rates were 0.018 and 0.033 cm year−1 at the shallow- and deep-water sites, respectively. We discovered an increase in light planktonic diatoms (Cyclotella) and a decrease in heavily silicified euplanktonic Aulacoseira through time at deep-water sites, related to more recent warmer air temperatures and shorter periods of lake-ice cover, which led to pronounced thermal stratification. Diatom beta diversity in shallow-water communities changed significantly because of the development of new habitats associated with macrophyte growth. Mercury concentrations increased by a factor of 1.6 since the mid-nineteenth century as a result of atmospheric fallout. Recent increases in the chrysophyte Mallomonas in all cores suggested an acidification trend. We conclude that even remote boreal lakes are susceptible to the effects of climate change and human-induced pollution.

intermediate water depth. Chronologies for two cores were established using 210 Pb and 137 Cs. Sedimentation rates were 0.018 and 0.033 cm year -1 at the shallowand deep-water sites, respectively. We discovered an increase in light planktonic diatoms (Cyclotella) and a decrease in heavily silicified euplanktonic Aulacoseira through time at deep-water sites, related to more recent warmer air temperatures and shorter periods of lake-ice cover, which led to pronounced thermal stratification. Diatom beta diversity in shallow-water communities changed significantly because of the development of new habitats associated with macrophyte growth. Mercury concentrations increased by a factor of 1.6 since the mid-nineteenth century as a result of atmospheric fallout. Recent increases in the chrysophyte Mallomonas in all cores suggested an acidification trend. We conclude that even remote

Introduction
In the debate about the role of humans in changing our environment, it is necessary to study the role of natural processes versus the impact of anthropogenic activities. The overwhelming role of human activity basically comprises two aspects: (1) climate change from enhancement of atmospheric greenhouse gases by combustion of fossil fuels, and (2) pollution of the natural environment by industrialization and other activities associated with human population growth (IPCC 2018). Whereas atmospheric change acts mainly at a global scale, the impact of pollution occurs mainly at regional and local scales. Global influences on the environment can be recognized in pristine ecological systems that are not directly affected by local anthropogenic activities, such as in the vast boreal areas of the Northern Hemisphere (Jones et al. 2015).
Since the onset of industrialization in the midnineteenth century, accelerated increases in air temperature (Abram et al. 2016) and human land-use have induced ecological consequences in northern freshwater ecosystems (AMAP 2017). Boreal and highelevation regions are particularly sensitive to climate change because of strong snow albedo and surfacebased feedbacks (Pepin et al. 2015). Over the past decades, Russian landscapes have experienced faster warming than sites elsewhere (Biskaborn et al. 2019b). In these regions, lakes are effective sentinels of environmental change that rapidly integrate information about catchment changes.
Climate change is altering physicochemical and biological lake properties (Heinecke et al. 2017;Vyse et al. 2020;Wang et al. 2015). The cumulative effect of warmer air temperature on Russian lake ecosystems and the role of local factors that influence lake catchment responses remain unclear. Understanding the amplitude of lake responses to global and local environmental changes is crucial for the sustainability of these fragile ecosystems. A powerful approach for deciphering the impacts of climate change and industrialization on these remote northern lakes is paleolimnology. To infer changes, lake-sediment cores are dated and analyzed for biological remains and physicochemical properties (Cohen 2003).
In northern lakes, diatoms (unicellular siliceous algae) often represent the main aquatic primary producers. Their opaline frustules are well preserved in the sediment record ) enabling identification at the species level. The presence and variety of diatom taxa in aquatic environments are linked to the physicochemical variables of water (e.g. light, temperature, nutrients, pH, salinity, and the content of oxygen and organic matter) (Smol and Stoermer 2010;Van Dam et al. 1994). Collectively, these characteristics make diatoms effective indicators of environmental change and are especially useful bioindicators in cold-region water bodies such as the Russian boreal lakes (Biskaborn et al. , 2016Hoff et al. 2015;Huang et al. 2020;Palagushkina et al. 2017). Moreover, widespread responses of planktonic taxa to anthropogenic warming indicate that northern lakes have already crossed sensitive ecological thresholds (Rühland et al. 2008;Smol and Douglas 2007).
One of the well-known effects of industrialization in the Northern Hemisphere is acidification of lakes as a consequence of the emission, transport, and deposition of sulfuric and nitric acid aerosols by river input or directly from the atmosphere (Farmer 1990). Atmospheric fallout of anthropogenic pollutants, including heavy metals, leads to contamination of soils and subaquatic substrates. Historically emitted mercury (Hg) is a well-established tracer of heavymetal pollution in lakes (Jones et al. 2015;Wang et al. 2004). Although it is mainly adsorbed to sediment, it is still a dangerous threat to aquatic ecosystems and can affect the composition of, for example, diatom assemblages (Pérès et al. 1997). Modern atmospheric Hg accumulation in lake beds of remote and formerly pristine aquatic systems is a strong sign of contamination, and can be up to 24 times higher than preindustrial levels (Lamborg et al. 2002;Wang et al. 2004).
Lake Bolshoe Toko is possibly the deepest large lake in Yakutia, eastern Russia (Zhirkov et al. 2016). The lake is considered a pristine system to study the impact of environmental changes without direct anthropogenic activity (Biskaborn et al. 2019a). As shown in this and other studies, within-lake distributions of environmental indicators such as microfossils differ across space, and mixing processes do not homogenize the microfossil assemblages before deposition in such large lakes (Stewart and Lamoureux 2012). It is therefore possible to study the different responses to external environmental changes in deepand shallow-water habitats (Biskaborn et al. 2019a).
In this study, our primary goal was to assess the effects of global climate change and industrialization on remote Lake Bolshoe Toko. We also sought to examine if climate change after the mid-nineteenth century affected diatom assemblages, and if there were different responses in deeper and shallower water communities related to species habitat preferences. Therefore, the objectives of this study were to: (1) reconstruct diatom communities from sediment cores taken at different water depths, to evaluate differences before and after the onset of industrialization and (2) discover whether there was a change in sediment geochemical composition related to post-industrial human activity such as air temperature warming, or heavy metal pollution.

Study site
Lake Bolshoe Toko (Russia, Sakha Republic, 56°15 0 N, 130°30 0 E, 903 masl) is a deep oligotrophic freshwater lake, 15.4 km long and 7.5 km wide, with a surface area of about 84.3 km 2 and a maximum water depth of 72.5 m (Fig. 1). The lake is fed by the Utuk River, which runs through Lake Maloe Toko and brings water from the southern igneous catchment of the Stanovoy Mountains. Historical documents for the region suggest that no permanent settlement was ever established in the study area, although a temporary mining settlement located 17 km northwest of the lake ( Fig. 1) was built in 2011 CE and was still active at the time of our sampling campaign in 2013. The El'ginsky coal deposits will likely cause future impact to the lake ecosystem and its catchment (Konstantinov 2000;Sobakina and Solomonov 2013).
The study area is situated in the East Siberian boreal forest, within the continental temperate climate zone and a predominant westerly wind system (Shahgedanova 2002). About 8 km NNE of the lake, the Toko RS meteorological station (131°01 0 E, 56°10 0 N, WMO code 31137) has collected data on temperature and precipitation since 1960 CE, which were used in our study to validate the most recent climate trends.  Konstantinov (2000) and slightly modified from Biskaborn et al. (2019a) to fit the measured water depths at sites, and topographic features based on Kornilov (1962) The lake basin is of tectonic and glacial origin, adjoined to the northern slope of the eastern Stanovoy Mountain Range (Imaeva et al. 2009). The catchment is characterized by Precambrian igneous rocks in the south and Mesozoic sandstones and mudstones in the north. The northern part of the lake is bordered by moraines of three different glacial sub-periods (Kornilov 1962) (Fig. 1).

Fieldwork
In a collaboration between the Alfred Wegener Institute Helmholtz Centre for Polar and Marine Research (AWI) and the North Eastern Federal State University in Yakutsk (NEFU), we conducted fieldwork between 19 March and 14 April 2013, as part of the German-Russian expedition ''Yakutia 2013.'' To retrieve short sediment cores, we drilled 250-mmdiameter holes in the ice cover using a Jiffy ice auger. Water depth was measured using a portable Echo sounder (HONDEX PS-7 LCD) and a calibrated rope for verification. We retrieved short cores from different sites in the lake at different water depths using a UWITEC gravity corer (Fig. 1, ESM Table 4). We selected four cores for study, taken at shallow, intermediate and deep sites, to explore the responses of diatom communities to climate change, as related to lake bathymetry. We subsampled them onsite to avoid sediment mixing during transport. Sediment samples were transported in sterile Whirl-Pak TM bags to the AWI laboratories in Germany and stored at 4°C in a dark room prior to freeze-drying and analyses. Hydrochemical measurements were performed during fieldwork and are given in detail in Biskaborn et al. (2019a).

Dating
Chronologies were developed for short cores PG2139 and PG2203 by measuring 210 Pb, 226 Ra, 137 Cs and 241 Am activities in freeze-dried sub-samples 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). Counting times varied from 1 to 3 days depending on the activity of the sample. The criterion used was the number of 210 Pb counts, generally at least 500 for sufficient accuracy. 210 Pb was determined via its gamma emissions at 46.5 keV, and 226 Ra by the 295 keV and 352 keV gamma rays emitted by its daughter radionuclide 214 Pb, following 3 weeks storage in sealed containers to allow radioactive equilibration. 137 Cs and 241 Am were measured by their emissions at 662 keV and 59.5 keV, respectively. The absolute efficiencies of the detectors were determined using calibrated sources and standard sources of known activity. Corrections were made for the effect of self-absorption of low-energy gamma rays within the sample (Appleby et al. 1992).

Diatom analysis
We analyzed diatoms in 0.5-cm sediment slices taken from short cores, between the sediment surface (0 cm) and down to as much as 13 cm core depth. Samples were collected from the following four cores, each taken at a different water depth: PG2208 (68.3 m), PG2203 (36.1 m), PG2133 (26.0 m), and PG2139 (12.1 m). In total, 26 samples, each 0.1 g dry sediment, were retrieved from the sediment cores for diatom analysis. Each sample spanned 0.5 cm in the core, except for core PG2133, which had 1-cm slices removed above 4 and below 6.5 cm. Diatom slides were prepared and analyzed for species assemblages following standard procedures .
We treated all samples with hydrogen peroxide for removal of organic matter and added 5 9 10 6 microspheres to estimate diatom valve concentration (DVC). For species identification, we used several diatom floras, including Hofmann et al. (2011) and Krammer andLange-Bertalot (1986-1991). For rare and endemic taxa, we conducted a literature search, which included  for Pliocaenicus. A minimum of 300 (but up to 361) valves were counted in each sample, except one, using a Zeiss AXIO Scope.A1 light microscope with a Plan-Apochromat 100 9 /1.4 Oil Ph3 objective at 1000 9 magnification. Identifications of small species were validated with a scanning electron microscope (SEM) at the GeoForschungsZentrum Potsdam. In the sample from 8.5 to 8.0 cm in core PG2203, diatoms were scarce and only 204 valves could be counted. Diatom valve preservation was estimated while counting and calculated as the ''F-index'' (Ryves et al. 2001). In all samples, scales of the chrysophyte Mallomonas were counted along with diatom valves, but not identified to species level. The Mallomonas index, i.e. Mallomonas relative to diatom cells, was estimated to assess lake acidification, following Smol (1984).
The change in relative abundance of Cyclotella (sensu lato), and two groups of Aulacoseria taxa (group 1: A. distans & A. valida and group 2: A. ambigua, A. subarctica) at the onset of industrialization was calculated by subtracting mean percentages in bottom samples (before * 1850 CE) from the mean in top samples (after * 1850 CE) for each core. According to the number of samples present and inferred accumulation rates, samples used to calculate means varied (top versus bottom samples as indicated red and blue in Fig. 3: PG2208: 4 vs. 4, PG2203: 3 vs. 3, PG2133: 2 vs. 2, PG2139: 2 vs. 2). To exclude the influence of earlier climate fluctuations, before the pre-industrial conditions, from the calculations, we did not extrapolate dates from 210 Pb age model, but estimated sedimentation rates instead to infer a higher bottom-boundary for the samples from shallow-water sites.
There is no 210 Pb age model available for cores PG2208 and PG2133, therefore we used the dated cores (PG2139 and PG2203) to estimate approximate ages in the undated cores, using lithological similarities between the two deepest cores (PG2208, PG2203) and the two shallowest cores (PG2133, PG2139).

Organic geochemistry
To assess the accumulation of organic matter in the lake, we analyzed total organic carbon (TOC) and nitrogen (N) in 61 freeze-dried and milled samples from the top 30.5 cm of short core PG2203 (Fig. 1) using a Vario EL III CNS and a Vario MAX C analyser. The measurement accuracy was 0.1 wt% for both TOC and TN. To gain additional information about the productivity, we analyzed stable organic carbon isotopes (d 13 C) in 61 samples using a Finnigan Delta-S mass spectrometer. Results are presented as d 13 C values relative to the PDB standard (%) with an error of ± 0.15%. The information on carbon chemistry is provided at a resolution of 0.5 cm in the short core.

Inorganic geochemistry
To characterize the chemical variability of the lake sediment in historical times, we analyzed samples from dated short core PG2203, taken at intermediate water depth. This site is located in the central part of the lake, at a depth close to the average water depth for the whole lake, and thus is believed to best represent the sediment chemical properties (Fig. 1). For analysis of heavy metal pollution, we used mercury (Hg) concentrations. We sent 14 samples from sediment core PG2203 to the chemical-analytical laboratory of Vinogradov Institute of Geochemistry, Siberian Branch of the Russian Academy of Science (IGC SB RAS). The laboratory measured mass fraction of total mercury using an RP-91S/RA-915 ? mercury analyzer. Total sedimentary mercury was detected by atomization of mercury in a two-section pyrolyzer and subsequent determination by the flameless atomic absorption method with the mercury analyzer. Measurement error was 0.03-0.06 ppm and is shown next to the data.

Statistical analyses
For multivariate statistical processing of diatom data we used the R environment (R Core Team 2016). A constrained incremental sum-of-squares clustering (CONISS) on all species counts was performed after chi-squared standardization and use of the Bray-Curtis dissimilarity index. A principal component analysis (PCA) was performed on species data after filtering out rare species. Only species with C 10 valve counts in C 2 samples were kept in the dataset for PCA and included in the species distribution diagrams and biplots of the first two principal components (ESM Fig. 1-4). We chose PCA as a suitable ordination method because the gradient length was below 2.5 standard deviation (SD) units (Birks 2010), as determined by detrended correspondence analysis (DCA). PCA was performed using the 'prcomp' function in the basic R package 'stats'. DCA was performed using the function 'decorana' in the package 'vegan' (Oksanen et al. 2012). DCA axis lengths from the core datasets are as follows (in SD units): PG2208 = 0.93; PG2203 = 1.18; PG2133 = 1.25; PG2139 = 1.57. The total variances of the core datasets explained by the first two principal components derived from the PCA are: PG2208 = 61.05%; PG2203 = 81.67%; PG2133 = 69.43%; PG2139 = 69.22%.
Diatom species richness was estimated after sample-size normalization, using a rarefaction analysis of Hill numbers in the iNEXT package (Hsieh et al. 2016). To estimate the rate of change of species composition between samples, we estimated beta diversity using the function beta.div in the R package 'adespatial' (Dray et al. 2016). This function estimates the total beta diversity (BDTotal) from the total variance in the community data and decomposes the BDTotal into local contributions to beta diversity (LCBD). Adjusted p values were estimated using the 'Holm' correction method in the same R package (ESM Table 3). Prior to the estimation of beta diversity, count data were transformed using Hellinger dissimilarity. Relative changes for other sediment variables, as a consequence of industrialization were calculated using an approach similar to that used for diatom taxa, that is, by subtracting the mean of samples before * 1850 CE from the mean of samples after * 1850 CE for each core.

Dating
Sedimentation rates in PG2139 and PG2203 are very low (Fig. 2, Table 1). 210 Pb/ 226 Ra equilibrium, corresponding to a period of around six 210 Pb half-lives (* 130 years), occurs at sediment depths of \ 5 cm. In PG2139 the unsupported 210 Pb activity declines exponentially with depth (Fig. 2e), showing that sedimentation rates have remained relatively constant since at least the late nineteenth century. The mean accumulation rate determined from the exponential gradient of the 210 Pb decline is 0.018 ± 0.002 cm year -1 . A higher accumulation rate of 0.033 ± 0.006 cm year -1 in PG2203 (collected from 36.1 m water depth) is probably a consequence of sediment focusing, as supported by its higher mean 210 Pb supply rate (74 Bq m -2 year -1 ) compared to PG2139 (45 Bq m -2 year -1 ).
The 210 Pb dates are in reasonable agreement with the 137 Cs results. High concentrations of this artificial radionuclide, generated by the atmospheric testing of thermo-nuclear weapons, were detected in the uppermost 2 cm of PG2203. Presence at similar depths of low but significant concentrations of 241 Am, also a product of nuclear weapons fallout, shows that this section of the core includes material from the early 1960s. The 210 Pb dates place 1963, the year of maximum 137 Cs fallout, within the 1.5-2 cm slice. Highest 137 Cs concentrations in PG2139 occur in the uppermost 1 cm, reflecting its slower accumulation rate. The 210 Pb dates for this core place 1963 CE within the 0.5-1 cm slice.
The low resolution of our age models only enables us to estimate the timing of the industrial signal in the lake. Therefore, we performed a top-bottom comparison between pre-industrial and industrial conditions, as reported in the literature for the Northern Hemisphere (Abram et al. 2016;Gabrielli et al. 2020;Holtgrieve et al. 2011). According to this literature, we chose 1850 CE to approximate the onset of industrialization. We are aware of the fact that industrial activity in Europe started as much as 100 years earlier (IPCC 2018). Noticeable effects, however, were delayed (Holtgrieve et al. 2011) and are best expressed since the mid-nineteenth century (Abram et al. 2016). Accordingly, we used the 210 Pb and 137 Cs chronologies for cores PG2203 and PG2139 to identify the 1850 CE boundary between top and bottom samples and projected the boundaries onto cores PG2208 and PG2233.

Diatoms
Preservation of diatom valves was generally excellent, with a mean F-index of 0.9, which means 90% of the valves counted were pristine (Flower and Ryves 2009). In sediment samples from core PG2208 (Fig. 3a), taken at 68.3 m water depth, 70% of diatom valves are planktonic, 25% are epiphytic, and 5% are epipelic species (means of all counted samples). Cyclotella taxa are represented by 45.2% and Aulacoseira by 7.4%, on average. The difference in mean Cyclotella abundance between the upper four samples, from 0 to 2.5 cm, and the lower four samples, from 3 to 13 cm, is 1.4%. The differences from the same sample intervals for Aulacoseira group 1 (A. distans & A. valida) and Aulacoseira group 2 (A. ambigua, A. subarctica, A. sp.) are -0.5% and -4.7%, respectively (Fig. 5). Species richness increases from 38.6 to 40.7, DVC (mean 3.9) increases to 6.0 10 8 g -1 at 0.75 cm depth, and the Mallomonas index increases from 25.6 to 35.7. Pliocaenicus bolshetokoensis is present at 14.9% on average in all samples, with no clear trend, whereas Achnanthidium minutissimum averages 8% and exhibits an increasing trend in the uppermost samples. Total beta diversity (BDTotal) is 0.17; LCBD varies between 0.09 and 0.16, but none of the values is statistically different. CONISS cluster analysis shows a clear distinction between surface samples and samples from the lower part of the sediment core.
Sediment core PG2203 (Fig. 3b), from 36.1 m water depth, has 68% planktonic, 13% epiphytic, and 19% epipelic species (means of all counted samples). Cyclotella taxa are represented by 46.6% and Aulacoseira by 6.2%, on average. The change in the mean abundance of Cyclotella species between three upper samples, from 0 to 4.5 cm, and three lower samples, from 6 to 13.5 cm, is 9.3%. The differences from the same sample intervals for Aulacoseira groups 1 and 2 are 0.2% and -2.4%, respectively (Fig. 5). Species richness is 49.5 on average in all samples, with strong fluctuations (35.6-82.9), but no clear trend. DVC (mean 5.1) increases to 10.4 9 10 8 g -1 at 0.25 cm depth. The Mallomonas index only applies to the two uppermost samples (mean 17.5). Pliocaenicus bolshetokoensis is present at 13% on average, with a slight decrease in the uppermost sample (9.5%). Achnanthidium minutissimum reaches 6.5% on average. BDTotal in the core is 0.20 and LCBD varies between 0.1 and 0.22, but none of the values is statistically different. CONISS cluster analysis shows a clear distinction between the surface sample and the lower part of the sediment core.
Sediment samples from core PG2133 (Fig. 3c), collected at 26 m water depth, have 70% planktonic, 24% epiphytic, and 6% epipelic species. Cyclotella taxa are represented by 37.7% and Aulacoseira by 11.0%, on average. The change in the mean abundance of Cyclotella species between two upper samples, from 0 to 3 cm, and two lower samples, from 4 to 6.5 cm, is 5.5%. The differences from the same sample intervals for Aulacoseira groups 1 and 2 are 2.2% and -4.9%, respectively (Fig. 5). Species richness has a mean of 32.5 in all samples without any clear trends. DVC (mean 4.3) peaks at 8.5 cm, with 9.6 9 10 8 g -1 , then decreases, but peaks again at 0.5 cm, with 6.4 9 10 8 g -1 . Mallomonas is only present in the surface sample, giving an index value of 16.3. Pliocaenicus bolshetokoensis is present at 20.8% on average, with no clear trend, but a notably low value in the surface sample (15.3%). Achnanthidium minutissimum accounts for 7.4% on average; Rossithidium pusillum is present in some samples (max. 4.6%), showing a steady decrease towards the top of the core, before disappearing in the surface sample. BDTotal in the core is 0.16 and LCBD varies between 0.12 and 0.20, but none of the values is statistically different. CONISS cluster analysis indicates a separation of the two uppermost samples.
In sediment samples from core PG2139 (Fig. 3d), taken from 12.1 m water depth, 37% of valves are planktonic, 47% epiphytic, and 16% epipelic species. Epiphytic species exhibit a slight increase in the record; planktonic species show a peak at 2.25 cm. Cyclotella taxa are represented by only 1.7% and Aulacoseira by 18.9%, on average. The  Fig. 3 Diatom assemblage changes in short sediment cores from Lake Bolshoe Toko. The red shaded area indicates the period of global industrialization starting in the mid-nineteenth century, as determined by 210 Pb and 137 Cs dating. The blue area indicates samples used to represent pre-industrial values. Species occurrence is given as relative abundance (%, scales vary to show rare taxa) and species richness as rarefied number of taxa. Planktonic/benthic ratio, Mallomonas index, and local contribution to beta diversity (LCBD) are also shown change in the mean abundance of Cyclotella species between two upper samples, from 0 to 2.5 cm, and two lower samples, from 4 to 6.5 cm, is negligible (-0.03%). The differences from the same sample intervals for Aulacoseira groups 1 and 2 are -8.2% and 0.4%, respectively. Shallow site PG2139 therefore exhibits the strongest increase in Aulacoseira species (Fig. 5). Species richness has a mean of 47 in all counted samples, and an increasing trend, reaching its maximum value in the surface sample (52.6). DVC (mean 4.9) generally decreases from 6.7 9 10 8 g -1 (13.5-8 cm), but slowly rises again and peaks at 5.5 9 10 8 g -1 at 0.25 cm.

Geochemistry
To gain reliable information on the chemical changes related to industrialization, the organic and inorganic geochemical composition was characterized for the top 30.5 cm of dated short core PG2203 (Fig. 4). The concentration of total organic carbon (TOC) is 5 wt% on average (3.7-6.5 wt%), with only minor decreases of about 1% at 15 cm, 6 cm, and in the uppermost 2 cm. Mean nitrogen (N) is 0.45 wt% (0.35-0.54 wt%). Changes in N are subtle and generally parallel changes in TOC between 30 and 4 cm. In the uppermost 4 cm of the core, N, in contrast to TOC, increases to maximum values. Accordingly, the TOC/ N atomic ratio (9.58-15.6, mean 12.9) exhibits high values in the lower part of the core and distinctly lower values in the uppermost 3 cm. Stable carbon isotope (d 13 C) values, relative to the PDB standard, have a mean of -26.85%, with a minimum value of -27.18% and a maximum value of -26.53%. Values tend to decrease slightly and steadily from the bottom of the core towards 6 cm, after which there is a remarkable increase in d 13 C, with high values in the uppermost 4 cm (mean -26.73%).

Chronological constraints
Because of very low sediment accumulation rates in lakes of northern cold regions, the sediment record in Lake Bolshoe Toko provides limited temporal resolution. Each 0.5 cm slice represents * 15 years of accumulation in PG2203 and * 30 years in PG2139, and the mean diatom sample resolution (1.7 cm) is * 52 years in PG2203 and * 94 years in PG2139 (Table 1). Hence, we used the short core data only for a top-bottom comparison, by averaging samples above and below the sediment depth dated to 1850 CE. This approach is supported by 210 Pb results that suggest bioturbation is not a major problem. Nevertheless, another possible consequence of the low sedimentation rates is that even a small amount of mixing can blur environmental signals in the short cores. Therefore, we used the data only to detect qualitative changes and trends since the onset of industrialization.
Sediment accumulation rates differ between deeper and shallower sites (Table 1). We attribute this to winnowing (Crusius and Anderson 1995). Lake Bolshoe Toko is a flow-through lake with only one outflow and we suspect stronger currents in the northern shallower part. The central basin is deep and calm, so that more sediment particles settle, and are possibly augmented by focusing effects. Particle focusing may be intermittent, as suggested by the irregular nature, i.e. fluctuating activity, of the 210 Pb record from PG2203 (Fig. 2a, b).

Diatom responses to post-industrial warming inferred at different water depths
The cores at water depths between 68.3 and 26 m (PG2208, PG2203, PG2133) reveal similar trends, with an increase in Cyclotella and decrease in Aulacoseira species during the mid-nineteenth century (Fig. 5), the approximate onset of the industrial period (Abram et al. 2016). This general change in planktonic taxa is a response of deep-water pelagic groups to environmental forcing. The shifts in single species, as tested in the PCA biplots (Electronic Supplementary Material [ESM] Figs. 1-4), provide only limited potential for more precise inferences of ecological forcing. A careful interpretation of the PCA results reveals a slight tendency towards clustering of Cyclotella species (i.e. PG2203), and a slightly similar response for Pliocaenicus and some Aulacoseira species in PG2203 and PG2139 cores.
The genus Aulacoseira is abundant in deep, cool, and unstratified boreal lakes (Laing and Smol 2003). Species vary from euplanktonic (A. ambigua and A. subarctica) to tychoplanktonic (A. distans and A. valida) forms with strongly silicified, heavy, and thus rapidly sinking frustules, which require water turbulence to remain suspended in the photic zone (Gibson et al. 2003;Rühland et al. 2008). Cyclotella on the other hand represents small planktonic diatoms with less silicified and thus lighter frustules, with high surface area to volume ratios (Rühland et al. 2003). These taxa therefore prefer strong thermal stratification regimes and thrive under a warmer climate (Rühland et al. 2015;Smol et al. 2005). Moreover, Cyclotella sensu lato as a group can be used to detect recent warming in lakes that stratify (Reavie et al. 2017). Accordingly, diatom assemblage shifts at Fig. 4 Geochemical properties of sediment core PG2203, showing total organic carbon (TOC), nitrogen (N), atomic TOC/ N ratios, stable carbon isotope values (d 13 C), and mercury concentrations (Hg) in parts per million. The red shaded area indicates the industrialization of the Northern Hemisphere, starting in the mid-nineteenth century as determined by 210 Pb and 137 Cs dating. Ages are given in common era next to core depths. Dashed lines indicate mean values for the whole sequence deeper-water sites, from euplanktonic Aulacoseira (A. group 2) to Cyclotella (Fig. 5), indicate prolonged thermal stratification (Rühland et al. 2015). Among other common driving factors, such as light availability, wind speed, turbidity, and nutrient availability, or a combination of these factors, diatom shifts were mainly related to the response of Lake Bolshoe Toko to increasing air temperature, associated with shorter ice cover and longer summer thermal stratification. Influence of recent atmospheric warming on planktonic diatoms, caused by stronger thermal stratification, has also been observed in Lake Baikal ) and in lakes of western Russia (Solovieva et al. 2008).
The slight increase in A. minutissimum observed in PG2208 (Fig. 3a), which is regarded as a pioneer species, can indicate an increase in available habitats, including for diatoms that have tychoplanktonic and benthic lifestyles, during longer growing seasons. In contrast to the deep-water sites, core PG2139 reveals how shallower-water communities respond to warming. A distinct increase in species richness suggests that the development of new habitats, possibly resulting from macrophyte growth associated with changes Fig. 5 Diatom assemblage changes at the onset of industrialization in the mid-nineteenth century in Lake Bolshoe Toko. Changes in Cyclotella (sensu lato) and two groups of Aulacoseira (Aul.), calculated as difference between mean abundance before and after 1850 CE, as estimated by 210 Pb and 137 Cs dating of cores PG2203 and PG2139 and mapped onto cores PG2208 and PG2133. Genus Aulacoseira mainly consists of euplanktonic species (A. ambigua, A. subarctica and A. sp.) at deep sites (PG2203, PG2208). Shallow site PG2139 (and to a lesser degree intermediate site PG2133) reveals a higher proportion of tychoplanktonic A. distans and A. valida in light and nutrient availability during longer ice-free periods, leads to increasing diatom assemblage complexity (Rühland et al. 2015). At shallow site PG2139, euplanktonic A. subarctica still decreases (Fig. 3), while A. distans and A. valida increase (Fig. 5), as the latter species are considered littoral tychoplanktonic forms, often found in shallow waters with low transparency (Krammer andLange-Bertalot 1986-1991). The observed lower total beta diversity of diatoms in deep-water sites, compared to higher diversity in the shallow-water site (PG2139), corresponds well to the findings of Stoof-Leichsenring et al. (2020). Diatom metabarcoding of sedimentary DNA from Bolshoe Toko revealed a negative correlation between diatom beta diversity and water-depth-related differences in habitat and light availability.
We attribute the DVC increase in surface samples of all cores, i.e. in planktonic communities at deeper sites, to enhanced productivity as a consequence of warming Douglas and Smol 2010). However, this trend was initially interrupted at shallower sites, likely because habitat changes and associated species turnover caused stress to the diatom community and thus led to lower reproduction rates.
Our interpretation of the observed shifts in diatom assemblages, following recent air temperature increases and associated reductions in the duration of the annual ice-cover, is supported by instrumental measurements of air temperature since 1960 at the TOKO RS weather station, * 8 km NNE of Lake Bolshoe Toko (ESM Fig. 5). According to the data from the National Oceanic and Atmospheric Administration (NOAA) Database (www.noaa.gov; station ID RSM00031137, TOKO RS) mean summer and winter air temperatures (MSAT, MWAT) increased steadily until the date of sampling in 2013. During that period, MSAT increased by 0.4°C per decade and MWAT by 0.6°C per decade, suggesting a strong influence on lake ice. Our sediment records suggest that the onset of atmospheric warming at Lake Bolshoe Toko occurred about 100 years before acquisition of meteorological data from the ground-based instrumental weather station began.
Geochemical impacts to the lake ecosystem at the onset of industrialization Concentrations of mercury (Hg) in sediment core PG2203 increased since the late-nineteenth century (Fig. 4). Hg contamination in aquatic systems can come from atmospheric deposition, erosion sources, and urban, agricultural, mining, and combustion/ industrial discharges (Wang et al. 2004). We can dismiss significant impacts from the mining settlement 17 km NW of the lake because it was established only two years before core retrieval. Analyses of Hg in sediments of Lake Baikal (Gelety et al. 2007) and other lakes (Jones et al. 2015;Lamborg et al. 2002) showed that Hg concentrations in sediments of remote natural lakes reflect variations in the Hg concentration in the atmosphere. In industrial areas, Hg adheres to aerosol particles and is transported over thousands of kilometers from its source for up to one year before being deposited (Lucotte et al. 1995). It has been reported that mercury in lake sediments is sufficiently immobile to generate historically reliable downcore records (Lockhart et al. 2000) and that high metal concentrations in Russian lake sediments are amplified by acidification trends (Moiseenko and Gashkina 2011). Therefore, we interpret the increasing upcore Hg concentrations in core PG2203 to reflect greater late-nineteenth century atmospheric Hg fallout in the catchment of Lake Bolshoe Toko as a consequence of industrialization in the Northern Hemisphere. The increase in Hg concentration by a factor of 1.6 in the lake deposits is low compared to other core studies in the Northern Hemisphere, which show twofold to sevenfold increases (Schuster et al. 2002). Nevertheless, total Hg concentrations are about double those from short sediment cores collected in Kamchatka (Jones et al. 2015). The Hg peak occurs at 2.75 cm (1931 CE), with slightly decreasing values afterwards, which could reflect a decrease in heavy metal pollution since the 1980s (Thevenon et al. 2011).
Diatom community responses to heavy metal contamination can induce changes in abundance of certain indicator species, for example, increases in Nitzschia (Pérès et al. 1997;Tlili et al. 2011). Samples from Lake Bolshoe Toko, however, show only minor changes of Nitzschia, around 0.5-1.6% in the cores from deep-water sites (PG2208, PG2203, PG2133), and slightly larger changes at the shallow site (PG2139), around 1.6-4.0%, without clear increasing or decreasing trends. Thus, no impact of heavy metal contamination on the lake diatom assemblages can be deduced from the current data.
Parallel increases in values of d 13 C and Hg (Figs. 4,6) can indicate that less C3 plant material (Biskaborn et al. 2019a) and/or less terrestrial plant material is delivered to the sediments (Erlenkeuser et al. 1974). This raises the question of whether Hg is positively associated with algal-derived organic matter accumulation. An inverse relationship between d 13 C and planktonic/benthic diatom ratios was reported by Wang et al. (2013). In our records, however, only core PG2208 exhibits a decreasing planktonic/benthic ratio, and it remains unclear if increasing d 13 C values reflect autochthonous or allochthonous processes.
Increasing nitrogen concentrations coincide with decreasing ratios of organic carbon to nitrogen (TOC/ N atomic ), indicating an increase in nitrogen loading independent from carbon input to the surface sediments (Figs. 4, 6). Higher ratios of TOC/N atomic in the southern end of the lake, compared to the northern end, however, suggest differences in catchment characteristics and a high supply of terrestrial organic matter from the Utuk River (Biskaborn et al. 2019a). High N favours rapid reproduction rate of small diatom species such as Cyclotella (Jordan and Bender 1973). Aulacoseira subarctica, on the other hand, is a better competitor when there is high silica loading (Gibson et al. 2003). Higher nitrogen availability enhances growth of some less silicified Cyclotella taxa, even though they compete with other phytoplankton such as Mallomonas (Fig. 3). Furthermore, nitrogen enhances macrophyte growth (Farmer 1990), which explains the slight increasing trend for epiphytic species in PG2139, also demonstrated by increases in alpha and beta diatom diversity.
The delayed timing of nitrogen increase, starting ca. 1950 CE (Fig. 4), indicates that climate change was the first cause of observed diatom shifts (Korosi et al. 2020). Increased N likely originated from a combination of greater catchment weathering inputs and direct atmospheric fallout of nitric acid aerosols (Gruber and Galloway 2008). The reliability of Lake Bolshoe Toko for capturing and preserving records of atmospheric fallout is demonstrated by the quality and consistency of the 210 Pb, 137 Cs and 241 Am results.
Bioindicators in all short cores indicate secondary effects of catchment warming, i.e. a strong increase in the Mallomonas index (Figs. 3a-d, 6). Mallomonas is a group of flagellated golden-brown algae used as indicators of lake acidification ) and nutrient enrichment (Munch 1980). Lake acidification is supported by a decrease in R. pusillum, an alkaliphilous species (Bigler and Hall 2002), in sediment cores PG2233 and PG2139, and it even disappears in the surface sample at the shallow-water site.
According to Kremleva and Moiseenko (2017), susceptibility of lakes to acidification is largely determined by the prevalence of silicic rocks in the catchment, which enable the inflow of acid-forming compounds. They also report that lakes in the northern and central taiga are the most susceptible to acidification among lakes in western Siberia. Given the geological setting of Bolshoe Toko, we believe that in spite of low precipitation in the region (ESM Fig. 5), the observed lowering of alkalinity occurred in the course of warming and atmospheric acid supply during industrialization in the Northern Hemisphere. It must be noted that parallel increases in air temperature and acidification make it difficult to attribute changes in diatom relative abundances to a single driver, and it is more likely that the community responded to the combined effects of both. Vörösmarty et al. (2010) demonstrated that 80% of the world's population is exposed to threats of water security and associated declines in freshwater biodiversity, mostly related to densely settled areas with intense agriculture. In their study, remote areas such as Yakutia are subject to less threats because pollution comes mainly from atmospheric fallout. Nevertheless, our study demonstrates that even extremely remote freshwater ecosystems experienced human impact following the onset of nineteenth-century industrialization, with negligible delay.

Conclusions
This investigation of oligotrophic freshwater Lake Bolshoe Toko revealed biogeochemical responses to climate change and industrial pollution in the Northern Hemisphere. Although there is almost no direct human impact on the lake, it has experienced changes forced by climate and long-range transport of human-generated pollutants. Small-scale changes in the lake started with the onset of industrialization in the mid-nineteenth century, the timing of which is supported by well-preserved anthropogenic contaminants such as mercury, which increased in the lake sediments by a factor of 1.6.
Effects of limno-ecological change can be attributed to climate change since 1850 CE, which drove shifts in the diatom community, followed by acidification, as indicated by increases in Mallomonas. Changes in planktonic diatom abundance suggest increasing summer stratification, likely a consequence of rising winter and summer air temperatures and associated shortening of the duration of lake-ice cover. Species indicative of macrophyte habitat availability and nutrient loading increased at some sites. Deepwater communities underwent changes in key planktonic taxa relative abundances, i.e. from euplanktonic Aulacoseira to Cyclotella species, within otherwise relatively stable species assemblages, and significant species turnover (beta diversity) was found only at the shallowest site, which generally had a higher diversity of available habitats to support different life forms.
Regional warming is supported by local meteorological data that cover only the last few decades. The sedimentary diatom records enabled us to place this regional warming into longer-term perspective, extending back to the second half of the midnineteenth century. Our findings underline the high sensitivity of remote boreal lakes to human-induced climate change and pollution.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open Access funding enabled and organized by Projekt DEAL.
Data availability The full dataset is available at PANGAEA (https://www.pangaea.de; Data Publisher for Earth & Environmental Science).