Changing nutrient cycling in Lake Baikal, the world’s oldest lake

Significance Lake Baikal (Siberia) is the world’s oldest and deepest lake and a UNESCO World Heritage Site. Containing an exceptionally high level of biodiversity and endemism, in addition to a fifth of global freshwater not stored in ice sheets, the lake has been cited by UNESCO as the “most outstanding example of a freshwater ecosystem.” Using geochemical and climate data, we demonstrate that rates of nutrient supply to the lake’s photic zone have risen to unprecedented levels in the last 2,000 y through the 20th and 21st centuries. Linked to increases in wind speed enhancing deep ventilation, we show that these changes are capable of altering lake primary production and community dynamics, including the balance between endemic and cosmopolitan species.

constitute 50-90% of phytoplankton biomass in the lake and consume up to 90% of available silicic acid in the photic zone 25,37 . Following atmospheric warming of the surface waters a summer stratification develops, which leads to low surface water nutrient conditions and to diatoms being replaced by pico-cyanobacteria as the dominant phytoplankton 23,35,38 . As summer stratification breaks down in autumn, the return of turbulent mixing triggers a second diatom bloom that is notably smaller than the spring bloom in terms of magnitude and silicic acid utilisation 35,37 . Today, productivity in the photic zone of Lake Baikal is co-limited by nitrogen and phosphorus 39,40 . Therefore, increased upwelling of deep water rich in silicon, nitrogen and phosphate to the photic zone (630 mmol SiO2 m −2 year −1 ; 93 mmol NO3m −2 year −1 ; 5 mmol P m −2 year −1 ) 34 will increase both biological productivity and rates of silicic acid utilisation.
Silicon exported by diatoms from the photic zone will, depending on rates of dissolution and other sedimentation processes, either be remineralised in deep waters or transferred into the sediment record. Combined with a long silicon residence time in Lake Baikal of 100-170 years, this results in deep water nutrient concentrations for nitrate, phosphate and silicate that are significantly higher than those at the surface 34,[41][42][43][44] . The mixing of these waters back to the photic zone, therefore, provides an essential source of nutrients to fuel further primary productivity. This internal cycling of regenerated silicon back to the photic zone significantly exceeds riverine inputs to the lake of 312 mmol m -2 yr -1 33,34 , with riverine inputs quickly diluted by pre-existing lake waters 37 . Consequently, changes in the rate at which these nutrient-rich deep waters are mixed into the photic zone have the potential to significantly alter biogeochemical cycling, aquatic productivity and food-web interactions in Lake Baikal.
Previous studies have argued that silicic acid concentrations have decreased throughout Lake Baikal's water column from 1995-2001 CE, in response to an increased export of silicon by diatoms into the sediment record 44,45 . These measurements, however, are potentially biased by sampling, which varies in any given year from May to September and so is influenced by aforementioned seasonal/monthly variations in biological productivity and nutrient cycling to the photic zone. The observation of a decline in Lake Baikal silicic acid concentrations is also not in agreement with nutrient mass-balance calculations for the lake 34 . Within the δ 30 Sidiatom record presented in this study, there are insufficient samples which allow further investigation of the reported decline in silicic acid concentrations from 1995-2001 CE 44,45 , although our reconstructed record indicates a moderate increase in rates of silicic acid supply to the photic zone between two samples dated to 1990 CE and 2000 CE respectively.

Silicon isotopes
Silicon has three stable isotopes: 28 Si (92.2% relative abundance on Earth), 29 Si (4.7%), and 30 Si (3.1%) 46 . The silicon isotope ratios of diatoms (δ 30 Sidiatom) are expressed as δ and reported in part per mille (‰): where R is 30 Si/ 28 Si and "reference" refers to NBS-28: the standard reference material for silicon. There is a well constrained mass dependent fractionation between δ 30 Si ( 30 Si/ 28 Si) and δ 29 Si ( 29 Si/ 28 Si) such that: During biomineralisation, diatoms preferentially uptake the lighter 28 Si isotope over the heavier 29 Si and 30 Si from silicic acid in the water column. This results in a progressive enrichment of 30 Si in both the dissolved and particulate phases within the photic zone, so that δ !" Si #$%&'( can be used to infer changes in the rate of silicic acid utilization. The magnitude of this isotopic fractionation between diatoms and dissolved silicic acid is defined by (α): Most commonly, this isotopic fractionation is defined in terms of the enrichment factor ε: Whilst marine diatoms grown under laboratory conditions have shown evidence of speciesdependent silicon isotope fractionation 47 , such a process has not been demonstrated in lacustrine diatoms/systems. Instead, previous work analysing diatoms in sediment traps through Lake Baikal's water column has constrained ε at −1.61‰ 48 . Reassessment of the raw data, assuming a normal distribution for proxy data uncertainty and using 10,000 replicate Monte Carlo simulations with the Monte Carlo package in R 49,50 indicates an uncertainty for ε of 0.11‰ (1σ). Given this, and using the Rayleigh fractionation model to describe changes in the isotopic composition of δ 30 Sidiatom and δ 30 Si(OH)4 51 , measurements of δ 30 Sidiatom can be interpreted to reflect changes in the relative rate of nutrient (silicic acid) utilisation and supply in the photic zone (Fig. S2).
An increase (decrease) in δ 30 Sidiatom could be caused by: 1) an increase (decrease) in biogenic silicic acid utilisation due to the isotope fractionation associated with this process; 2) a decrease (increase) in nutrient supply to the photic zone, which replenishes the pool of nutrients and their isotope composition; or 3) a combination of these two processes. In calculating rates of photic zone silicic acid utilisation and supply (see Equations in Methods Section which are repeated below for clarity) a series of assumptions are required, each of which are examined in the sections below:

Dissolution
The use of (δ 30 Sidiatom) in palaeoenvironmental research requires that the isotopic composition of living diatoms in the photic zone is faithfully transported into diatoms buried in the sediment record following their deposition at the surface-sediment interface 52 . This issue is particularly important in Lake Baikal, with only c. 1% of diatoms in the photic zone ultimately preserved in the sediment record 53,54 and dissolution indices indicating that up to 60% of diatom frustules preserved over the last 1,000 years have undergone some form of dissolution 18,55 . However, measurements from sediment traps distributed down the entire depth of Lake Baikal's water column and at the surface sediment show no alteration in δ 30 Sidiatom, indicating that any dissolution during sinking or sedimentation in Lake Baikal does not modify δ 30 Sidiatom 48 .
The high rate of diatom dissolution does potentially impact sediment biogenic silica concentrations and so calculation of Si(OH)4 supply (equation 6). The Si(OH)4 supply record presented in this study assumes that the degree of biogenic silica dissolution through the water column and within the sediment record has remained unchanged over the analysed interval. This assumption can be tested through the use of a diatom dissolution index 56 : in which n is the number of diatom frustules (n ≥ 300), x1i is the number of pristine diatom frustules and x2i is the number of diatoms frustules showing any sign of dissolution. From this the raw BSi mass accumulation rate (MAR) concentration record can be corrected for dissolution and used to re-calculate changes in silicic acid supply to the photic zone (equation 6): The results from this show that whilst reconstructed rates of silicic acid supply are higher when the BSi data is corrected for dissolution, the overlying trend of a significant increase in nutrient (silicic acid) supply from the 20 th Century onwards remain valid and is in agreement with concordant increases in Ulnaria acus, a diatom often associated with high dissolved silica concentrations 57,58 (Fig. S3). In this study we elect not to use the dissolution corrected nutrient supply record due to the uncertainties associated with individual indices 59,60 . We also avoid using diatom biovolume accumulation rates in place of BSi MAR since, although relatively well established for Lake Baikal 61 , there remains considerable uncertainty in the calculation of species biovolume coefficients and over their potential to vary with time and changing environmental conditions, including changes in nutrient availability.
Although nutrients are primarily supplied to the lake through riverine inputs 33,34,62 , intra-and interannual nutrient cycling in Lake Baikal is regulated by vertical mixing in the water column, which transports regenerated nutrients from deep waters in Lake Baikal into the photic zone 33,34 . Due to this cycling and the volume of water in Lake Baikal, the isotopic and DSi signature of riverine waters entering the lake is minimal beyond the immediate area around a river mouth 37 . As such, the δ 30 Si composition of silicic acid (δ 30 Si(OH)4) supplied to the photic zone and utilised by diatoms is best described by the δ 30 Si composition of deep water (δ 30 Silake), constrained at 1.71‰ (1σ = 0.10‰) 37 . Whilst it is not possible to account for how δ 30 Silake may have varied over time, the residence time of Si(OH)4 in the south basin of Lake Baikal is 100-170 years 43,44 , minimising the risk that large changes have occurred over the studied time interval.
Core chronology All 14 C dating was completed at the Gliwice Radiocarbon Laboratory, Silesian University of Technology. The total organic carbon was prepared from the sediment by treatment with 0.5M HCl, subsequent rinse with DI water and drying. The samples were combusted in VarioMicroCube elemental analyser (Elementar™) coupled to AGE-3 graphitisation equipment 63 . 14 C concentration was determined in DirectAMS laboratory, Bothell, WA, USA 64 (Table S1).
New 14 C and existing 210 Pb dates 65 were combined using the Bacon software 66 to establish an age-depth model as well as an age interval for each analysed sample (Fig. S4). The curve IntCal13 was used for calibration of 14 C results 67 . For calculations, the cores were divided into 0.5-cm-thick sections. The priors used were accumulation rate 200 yr/cm of gamma distribution with shape 1.5. The default memory prior was used, which is a beta distribution with the parameters mem.strength = 4 and mem.mean = 0.7. The models were extrapolated to match the range of proxy data and the ages and uncertainties for required depths generated. Throughout the paper, the weighted mean age was used as an approximate age, along with 1σ uncertainty.

Ekman transport
Due to Lake Baikal's depth (maximum depth of 1,642 m) and climatic conditions, a key contribution to vertical mixing is given by periodic deep mixing events controlled by thermobaricity, i.e. the combined dependence of water density on temperature and pressure. According to this physical property of water, wind-driven inshore Ekman transport can trigger thermobarically unstable conditions, generating deep coastal downwellings that can reach down to the bottom of Lake Baikal. Deep mixing events typically occur in December/January, before ice formation on the lake, and in May-June after ice-out in late spring when the lake is weakly inversely stratified 41,42,[74][75][76][77][78][79] . Further deep water renewal events have been detected under-ice, although it remains unclear to what extent these were driven by Ekman transport before the onset of ice-cover on the lake 80 (Fig. S5). This trend is consistent with a progressively increased probability of more intensive deep ventilation events and so deep water nutrient supply to Lake Baikal's photic zone during the 20 th and 21 st Century and follows previous analyses that showed deep ventilation in Lake Baikal is particularly sensitive to climate change through changes in wind forcing 79 . Whilst the horizontal resolution of the CERA-20C data is relatively coarse (125 km), ERA5 reanalysis data 84 from 1979 CE, with a finer horizontal resolution (30 km), averaged over the same area covered by the CERA-20C data, accurately captures wind seasonality over Lake Baikal (Fig. S6), confirms the wind patterns and the inter-annual variability and trends described by CERA-20C over the same interval (Fig. S5), thereby providing validation of the coarser reanalysis product.
Equations for calculating Ekman transport in the lake are reported in the methods section of the main manuscript, with anomalies for a given downwelling season (s) of year (i) calculated relative to the mean Ekman transport from 1990-2000 CE (Fig. S7). By reasonably assuming that all variables in Equations 4 and 5 in the Methods Section of the main text are constant, except for wind speed parallel to Lake Baikal's coast (W), Ekman transport (M) anomalies can be calculated as the anomalies of the square of directional wind (W) over Lake Baikal around the barycentre of the lake:

\\\\\ (1990 − 2000, )
where the subscript c indicates the cumulative calculated across periods when the wind both blows from the north-east and is interrupted by winds from other directions for less than 1 day. The overbar indicates averaging over each downwelling season (s) of a given year (i), to obtain a quantitative measure of the potential intensity of deep ventilation in that season and year.

Diatom assemblage data
Changes in Ekman transport nutrient supply to the photic zone, from the late 19th Century onwards, coincide with changes in the diatom record (Fig. S8 a-e). In addition to reductions in the ratio of autumn/spring blooming taxa, significant declines also occur through the 20 th and 21 st Century in the wider diatom community (PCA axis 1) and composition turnover -the change in species composition and relative abundances -(DCCA axis 1) 85 .
Changes in composition turnover, beginning in the late 19 th Century and pre-dating other anthropogenic impacts on the lake, are strongly associated with reduced relative abundances for the endemic Crateriportula inconspicua (r = +0.75; p < 0.005), (Fig. S8f). Whilst the ecology of Cyclotella inconspicua is poorly understood, phytoplankton monitoring has found C. inconspicua cells to be most abundant during autumn overturn 53 . However, Cyclotella minuta is the only pelagic diatom in Lake Baikal with a notable autumnal bloom, attributed to its ability to survive warm summer surface water temperatures (SWT) that result in a competitive advantage at the end of summer stratification 86 . Although C. minuta is only weakly associated with declining DCCA score (0.141; p = 0.511), together C. inconspicua and C. minuta suggest a significant decline in autumnal productivity from the 19 th Century (Fig. S8 f-g).
In addition to the autumnal decline in C. inconspicua and, C. minuta, the increase in diatom composition turnover is strongly associated with increases in the endemic Aulacoseira skvortzowii (r = −0.68, p < 0.005) and cosmopolitan Ulnaria acus (r = −0.83; p < 0.005) (Fig. S8 h-i). A. skvortzowii is adapted to cold water conditions 87 whilst U. acus is a thermally tolerant taxon associated with higher dissolved silica concentrations in the water column 58,88 . Whilst previous work has attributed recent community changes in Lake Baikal's diatom community to warmer SWT and enhanced summer stratification (Fig. S8j) 55,89, the concordant change between nutrient supply and diatom community changes from the late 19 th Century, including increases in U. acus, suggests that community shifts over the last 150 years are a function of both nutrient supply and SWT. A. skvortzowii and U. acus also both require strong winds/currents in Lake Baikal for their cells to be transported from littoral to pelagic waters. The correlated increases in wind speed, Ekman transport driven nutrient supply, A. skvortzowii and U. acus therefore lends further support to our suggestion that increased photic zone nutrient supply has significantly impacted recent diatom changes in Lake Baikal (Fig. S8h-j), by accounting for 23.4% (p = 0.002) of community variation. All ordinations were calculated using CANOCO5 90 .
The other dominant taxon in Lake Baikal, the endemic Aulacoseira baicalensis, exhibits a weak curvilinear trend through the 20 th -21 st Century (Fig. 8k). Whilst A. baicalensis needs relatively high concentrations of both silicic acid and turbulence in the photic zone to grow, the taxa is also adapted to cold water, low light conditions 91 . Although speculative, the increase in A. baicalensis up to c. 1950 CE may be linked to increasing photic zone nutrient supply, with the more recent (post-1950 CE) decline linked to increasing SWT, causing its heavily silicified cells to sink deeper into the water column ( Fig. S8j-k).

Generalized additive model
Generalized additive models (GAM) were calculated with restricted maximum likelihood (REML) smoothness selection using the mgcv package in R 50,92,93 . To account for temporal autocorrelation, all models included a continuous-time first-order autoregressive (CAR(1)) process 94 . Intervals of significant temporal change in a GAM were detected using the first order derivative of the fitted trend and a 95% simultaneous confidence interval using the gratia package in R 50,94,95 . Additional breakpoint analyses were carried out using the segmented package in R 50,96 .

Fig. S1
. Dominant seasonal controls on diatom productivity in Lake Baikal a, Ice/snow cover over the lake from January to May inhibits diatom blooms, although important under-ice growth occurs, with an increase during ice break-up. b, Ice break-up from May onwards allows turbulent mixing within the uppermost sections of the water column. Combined with coastal downwelling (typical months = May-June and December-January), and the associated upwelling and build-up of nutrient-rich deep waters in the photic zone, large diatom blooms occur. c, Thermal stratification of the lake in summer reduces nutrient availability for diatoms, with productivity instead dominated by autotrophic pico-cyanobacteria. d, Breakdown of the stratification in autumn months restarts the turbulent mixing of nutrients into the photic zone. Due to the absence of coastal downwelling and deep water supply of nutrients to the photic zone, the autumnal bloom is typically smaller than the spring bloom. Fig. S2. δ 30 Sidiatom in cores BAIK13-1C and BAIK13-4F in Lake Baikal (see Fig. 1 in the main text for core location). Given the proximity of the two sites and their well constrained age-models, these records are combined to create a composite record which are then used to reconstruct changes in photic zone silicic acid utilisation/supply (see main text). Shaded polygons reflect the 1σ uncertainty derived from Monte Carlo simulations (10,000 replicates).    Fig. S8. Diatom community changes in Lake Baikal at core site BAIK13-7 18,55 (Fig. 1). a, CERA-20C 83 and ERA5 84 wind speed reanalysis anomaly data, relative to a baseline period of 1990-2000 CE, from the barycentre of Lake Baikal during the typical periods of downwelling (May-June and December-January). b, Changes in photic zone silicic acid supply (relative to a value of 100% at 2005 CE). c, Detrended canonical correspondence analysis (DCCA) axis 1 scores reflecting diatom composition turnover. d, Principal components analysis (PCA) axis 1 scores reflect changing diatom taxa composition. e, Ratio of autumn/spring taxa. f, Relative abundance of C. inconspicua. g, Relative abundance of C. minuta. h, Relative abundance of A. skvortzowii. i, relative abundance of U. acus. j, Mean lake surface water temperatures (SWT) (May-October) from shoreline locations across Lake Baikal 15 . k, Relative abundance of A. baicalensis. For all panels, trend lines (solid) and confidence intervals (dashed lines) show a generalized additive model fitted to each time series.

± 105CE
Dataset S1 (separate file). δ 30 Sidiatom and absolute analytical uncertainty (2σ) together with reconstructed rates of silicic acid utilisation and supply to the photic zone with the 1σ uncertainty derived from Monte Carlo simulations (10,000 replicates).
Dataset S3 (separate file). CERA-20C wind speed reanalysis data from the barycentre of Lake Baikal during the typical periods of coastal downwelling (May-June and December-January). All data are anomalies relative to a baseline period of 1990-2000 CE.
Dataset S4 (separate file). Modelled Ekman transport anomalies during the typical periods of coastal downwelling (May-June and December-January). All anomalies are relative to a baseline period of 1990-2000 CE.