Decoupling of silica, nitrogen and phosphorus cycling in a meromictic subalpine lake (Lake Iseo, Italy)

Silica (Si), nitrogen (N) and phosphorus (P) loads and stoichiometry are key factors controlling the trophic status of lakes and coastal seas. In the hydrographic network, lakes also act as biogeochemical reactors, controlling both nutrient retention and fluxes. This work aimed to examine the coupling of Si, N and P cycling, together with their stoichiometry in a deep meromictic subalpine lake (Lake Iseo, Northern Italy). Si, N and P mass budgets were calculated by quantifying loads in the inlets and in the outlet over a period of 30 months (May 2016−October 2018), in-lake sedimentation rates and net nutrients accumulation in the water body. Lake Iseo acts as a biogeochemical filter, which differentially retains the external Si, N and P loads. Retention of Si and P was similar (75–79%), but considerably higher than N (45%), evidencing a decoupling of their fate due to in-lake processes. This differential retention is likely to be exacerbated by meromixis which enhances Si and P accumulation in the monimolimnion, while impairing denitrification, thus limiting N removal. Such decoupling resulted in an increase of the N:Si and N:P ratios in both the epilimnion and in the outlet compared to the inlets, whereas the ratios decreased in the monimolimnion. As a result, there may be a stronger Si and P limitation of the photic zone, leading to a shift towards more oligotrophic conditions. This transient equilibrium could be impaired in the case of water overturn produced by extreme climate events—a highly relevant issue, considering that a growing number of deep lakes are turning from holo-oligomictic to meromictic as a result of combined eutrophication and climate change.

the minimum nutritional requirement of a given phytoplankton or macrophyte species (Paerl et al. 2016). Since the first report of the Organization for Economic Co-operation and Development, P has been recognised as the main limiting factor influencing lake primary productivity (Vollenweider and Kerekes 1982), after which N was also identified as a limiting nutrient, especially in marine ecosystems (Howarth et al. 2011). Less attention has been given to Si, although its role as an essential element for diatoms has been recognised (Struyf and Conley 2012;Carey et al. 2019). Furthermore, Si is acknowledged as playing a key role for several aquatic macrophytes species by increasing structural resistance against water flow, preventing the assimilation of toxic elements (as Cd + ) and protecting against herbivory (Humborg et al. 1997;Serediak et al. 2014;Schoelynck and Struyf 2016). More recently, studies have highlighted the multifactorial character of eutrophication, which depends on both stoichiometry and speciation of Si, N and P (Billen and Garnier 2007;Glibert 2017). A balanced diatoms growth in marine and freshwaters is considered to occur at the molar ratio C:N:Si:P = 106:16:16:1 and at C:N:Si:P = 106:16:40:1, respectively (Redfield et al. 1963;Billen and Garnier 2007;Dupas et al. 2015). When an imbalance in the reactive forms of Si, N and P occurs, harmful algal blooms can take place, causing detrimental effects on ecosystem suitability for fishing, aquaculture and bathing (Ittekkot et al. 2000;Billen and Garnier 2007;Seitzinger et al. 2010;Glibert 2017). Despite the coupling of Si, N and P, how these elements are processed together have not been fully addressed at ecosystem scale in inland waters, resulting in a lack of perception of eutrophication complexity, especially with regard to surges and development (Duarte et al. 2009;Howarth et al. 2011;Glibert 2017;Maranger et al. 2018).
Nutrient abundance and stoichiometry in aquatic ecosystems are affected by hydrology and biogeochemical processes in the aquatic system itself, and by river runoff, which undergoes wide variations and imbalances due to the soil and its use in the watershed (Elser and Hamilton 2007;Hillebrand et al. 2014). In recent decades, human activity has significantly modified the transport and storage of Si, N and P across heavily exploited watersheds by changing processes which control terrestrial inputs from natural to human-managed, including untreated sewage, partially treated wastewater, agricultural soil runoff and atmospheric deposition (Harrison et al. 2009;Baker et al. 2014;Serediak et al. 2014;Carey and Fulweiler 2016;Maranger et al. 2018;Royer 2020). Furthermore, climatic changes together with altered watershed hydro-geomorphology could increase the quantity of nutrients delivered to aquatic environments, thus changing their relative proportions (Strayer and Findlay 2010;Baron et al. 2013; Barone et al. 2019;Pilotti et al. 2021).
The hydrographic network is a complex system of diverse, interconnected aquatic subsystems which act as reactors and filters, not only transporting, but also transforming, nutrient stocks (Meybeck and Vörösmarty 2005;Maranger et al. 2018). In this framework, lakes and reservoirs function as biogeochemical reactors able to retain and regulate the dissolved nutrient input, converting it into biomass or transforming it into non-reactive species (Ittekkot et al. 2000;Harrison et al. 2009;Verburg et al. 2013;Frings et al. 2014). These processes can greatly differ among elements as a result of their different physical-chemical properties, metabolic pathways, and fate. P and Si have a predominantly sedimentary cycle with the particulate fraction which tends to settle. In open waters, dissolved silica (DSi) is continuously assimilated by phytoplanktonic diatoms and a few other species to form biogenic silica (BSi), which settles after organisms die, accumulating in the bottom water layers and sediments, where it is partially retained (Hobbs et al. 2010;Frings et al. 2014). Similarly, P retention is due to the assimilation by phytoplankton, along with physical-chemical processes, such as adsorption and co-precipitation, especially with metals (i.e. iron, aluminium and calcium), followed by sedimentation and sediment burial (Reynolds and Davies 2001;Serediak et al. 2014). In contrast, in heavily human-impacted systems, reactive N mainly occurs as dissolved nitrate  ) and, to a lesser extent, as dissolved ammonium (NH 4 + ), with both ions able to be assimilated by primary producers, leading to subsequent N sedimentation. In addition, under oxic conditions, NH 4 + can be oxidized by nitrifying bacteria to NO 3 − , which is highly soluble; on the other hand, under anoxic conditions, NO 3 − can be reduced via bacterial denitrification to di-nitrogen gas (N 2 ), through a multistep reduction which also includes the formation of the greenhouse N 2 O gas (Saunders and Kalff 2001;Baron et al. 2013).
Denitrification takes place mainly in the sediment both in littoral areas and in deep waters (Bruesewitz et al. 2012;Nizzoli et al. 2014Nizzoli et al. , 2018. Under anoxic conditions, anaerobic NH 4 + oxidation (ANAMMOX process) and dissimilative NO 3 − reduction to NH 4 + (DNRA) could also occur despite their role in lakes seems quantitatively less important compared to denitrification (Nizzoli et al. 2018).
Despite the recognition of lakes as biogeochemical reactors and nutrient traps, nutrient retention in lakes has mainly been explored with the focus on a single element (Duarte et al. 2009;Howarth et al. 2011;Glibert 2017;Maranger et al. 2018). A small number of studies on mass balance have investigated N:P stoichiometry in freshwater ecosystems (Howarth et al. 2011;Verburg et al. 2013;Maranger et al. 2018), with an even smaller amount of research focussing on the simultaneous fate of Si, N and P (Cook et al. 2010;Viaroli et al. 2013;Carey et al. 2019;Royer 2020). In addition, there is variable nutrient retention within lakes which depends on factors such as water residence time, temperature, redox potential and trophic status (Seitzinger et al. 2006;Finlay et al. 2013;Verburg et al. 2013;Nizzoli et al. 2018). Mixing rates and the occurrence and persistence of thermal stratification can drive in-lake biogeochemical processes, hence influencing nutrient retention (Rempfer et al. 2010;Baron et al. 2013;Nizzoli et al. 2018;Rogora et al. 2018). In meromictic lakes, for example, persistent stratification allows the formation of a deep-water layer, called the monimolimnion, which is usually anoxic, sometimes rich in H 2 S and sporadically mixes with water from the upper layers (Zadereev et al. 2017). Reduced water mixing prevents the upward flux of accumulated nutrients to the mixed layer, promoting the accumulation of the released nutrients in the monimolimnion (Reynolds and Davies 2001;Siipola et al. 2016;Nizzoli et al. 2018;Viaroli et al. 2018;Leoni et al. 2019;Lau et al. 2020). In the context of climate change, such a process represents a key issue, as temperate deep lakes may evolve towards a reduced frequency and depth of mixing events and even meromixis (Fenocchi et al. 2018), which would greatly alter nutrient pathways and fate.
The aim of this work was to analyse the coupling of Si, N and P biogeochemical cycles throughout their travel in a meromictic, deep subalpine lake. The specific objectives were (a) to evaluate how the lake controls and regulates the downstream fluxes of nutrients and their stoichiometry, thus acting as a filter in the hydrographic network; and (b) to explore the main processes within the lake that influence Si, N and P retention under meromictic conditions. It was hypothesized that retention differs for Si, N and P, with retention for P and Si being the highest and for N the lowest, as a result of differences in chemical binding properties and processes for each element in the sediment. A change in their stoichiometry, therefore, was expected to be detected between inlets and outlet. It was also hypothesized that silica metabolized by diatoms would undergo sedimentation and retention in the monimolimnion due to incomplete mixing. To verify these hypotheses, Si, N and P loads, their composition and stoichiometry were analysed in both the inlet and outlet waters. Internal lake mass budgets of each element were also calculated to assess to what extent meromixis affected their fate.

Study area
The study was carried out in Lake Iseo (Italy), a deep subalpine lake located on the southern slopes of the Alps (Fig. 1). Among the deep subalpine Italian lakes, Lake Iseo is the fourth largest in terms of volume and surface area (Salmaso and Mosello 2010).
The watershed of Lake Iseo is elongated, with an NE-SW orientation and covers a total area of 1807 km 2 (Pilotti et al. 2021). The average elevation of the lake watershed is 1401 m a.s.l. The catchment area is mainly characterized by sedimentary rocks, with only the northern part presenting crystalline formations (Bini et al. 2007). Population density around the lake is 109 inhabitants km −2 ; urban wastewater is treated by 50 plants, only five of which serve more than 10,000 inhabitants. Agricultural use accounts for 22% of the watershed surface (Table 1), mainly as pasture or rough grazing (386 km 2 , 21.4%). Vineyards (11 km 2 , 0.61%), chestnuts (6 km 2 , 0.33%) and corn (3 km 2 , 0.16%) represent the main crops. In term of livestock units (1 LSU = 1 adult dairy cow), livestock bulk accounts for 17 LSU km −2 ; cattle (18241 LSU), poultry (8145 LSU) and sheep (2307 LSU) are the main farmed groups (ISTAT 2010, http:// dati. istat. it/).
In the period 2000-2020, watershed rainfall measured in Sarnico ranged from 638 to 1660 mm year −1 , air temperature was between − 4.7 and 28 °C (Regional Agency for Environmental Protection of Lombardy Region) and annual average water inflow to the lake was between 33 and 84 m 3 s −1 (Supplemental Info Fig. S1). The main natural tributary is the river Oglio, contributing more than 38% of the lake's inflow (Fig. 1). An artificial channel, the Industrial Channel, which serves an electric power plant, has a similar discharge and accounts for 50% of the total water inflow. These two tributaries together drain a total area of 1444 km 2 . Other minor natural tributaries, including the Bagnadore, Borlezza and Rino torrents, drain an area of 307 km 2 (Pilotti et al. 2021). The water level in Lake Iseo is regulated by a dam located at the lake outlet in Sarnico, which was built in 1933 and is currently managed by the Oglio Consortium; the surface water outflow is regulated to supply water for irrigation and hydroelectric power production and to prevent the risk of flooding.
Like the other subalpine lakes, Lake Iseo has an elongated shape running from North to South (surface 60.8 km 2 ), with steep banks and a flat bottom in a cryptodepression-characteristics reflecting its glacial origin (Bini et al. 2007). The western shore is steep and rocky, while the eastern shore is steep in the northern part but becomes progressively flat towards the south. Since the early 1980s, Lake Iseo has been meromictic, although a complete overturn occurred in 2005 and 2006, and it is now classified as meso-eutrophic (according to the average phosphorus concentration in the mixolimnion); catchment nutrient load and meromixis are the main causes of this condition (Salmaso et al. 2018). The decrease of complete overturn events in Lake Iseo has been attributed to the effect of salinity on vertical density distribution (Ambrosetti and Barbanti 2005). Meromixis has been strengthened by the warming of the water column, which amplifies thermal differences between surface and bottom waters (Valerio et al. 2015). Since the last complete overturn, oxygen  (Northern Italy, left) and detailed map of Lake Iseo (right). The three sampling sites for loadings estimation are shown (•, Oglio inlet = Og. In; Industrial Channel inlet = Ch. In; Oglio outlet = Og. out), together with the sampling sites of littoral and deep benthic fluxes (♦), sedimentation fluxes (▲) and water depth profiles (■) depletion has been expanding upwards: currently, the water mass is anoxic below a depth of 100 m . Temperature, conductivity and oxygen profiles indicate that the water column can be divided into 2-3 compartments with a variable vertical range across years (Lau et al. 2020). As for the other subalpine deep lakes, a pronounced thermal stratification develops beginning in April, with a maximum stability in August-September . The epilimnion, comprised on average between 0 and 20 m (1.1 km 3 ), accounts for 14% of lake water volume; the surface area down to a depth of 20 m (52.4 km 2 ) accounts for 86% of the total surface of the lake. The hypolimnion, comprised between 20 and 90 m (3.2 km 3 ), accounts for 40% of lake water volume. Finally, the monimolimnion extends on average from a depth of 90 m downwards (3.6 km 3 ) and accounts for 46% of the lake water volume; the area at a depth of 90 m (36.7 km 2 ) accounts for 60% of the total surface of the lake.

Inlet and outlet nutrient concentration
The survey of the two main inlets and the outlet took place over a period of 30 months (May 2016−October 2018). A minimum of one sampling per month was carried out, totalling 37 dates for each of the three sites: the Oglio River and the Industrial Channel at the inlet, and near the Sarnico dam at the outlet ( Fig. 1). Samplings were carried out during both base flow and high flow events to account for the impact of flow variability on water chemistry. For each sampling date and station, three replicate water samples were collected with a custom made bottle, consisting in a PVC tube equipped with an one-way valve at the bottom that allowed to collect an integrated water sample of 1 m depth. Following collection, an aliquot for each replicate was immediately filtered (Whatman GF/F) and stored in polyethylene vials and in glass vials, the former ready for dissolved silica (DSi), ammonium (N-NH 4 + ) and nitrate + nitrite (referred as N-NO x − in the following) analyses, and the latter for soluble reactive phosphorus (SRP), total dissolved phosphorus (TDP) and total dissolved nitrogen (TDN) analyses. A second aliquot for each replicate was filtered on pre-weighed filters (Whatman GF/F), with the material collected on filters being subsequently analysed for the measurement of total particulate phosphorus (PP) and particulate nitrogen (PN). Total phosphorus (TP) was estimated as TDP + PP, while total nitrogen (TN) was estimated as TDN + PN. A third aliquot of water for each replicate was filtered using polycarbonate filters (Whatman Nucleopore) and then frozen to determine biogenic silica (BSi). The total silica (TSi) was estimated as DSi + BSi.
Si, N and P content and stoichiometry of settling particulate matter in the water column In order to quantify Si, N and P content and stoichiometry of settling particulate matter, two multitrap systems (Hydro-Bios Ltd., Kiel, Germany) were moored to collect total suspended solids at two depths: below the epilimnion (20 m) and near the hypolimnion-monimolimnion interface (90 m) ( Fig. 1, Lau et al. 2020). Each trap was equipped with 24 collecting bottles and two cylinders (154 cm 2 ) for duplicate sampling. The bottle content was preserved with pre-added formalin to prevent the decomposition of settled material. The traps were put in position from 8 April 2016 to 5 April 2017, with a 31-day sampling interval. Samples were immediately taken to the laboratory, stored in the dark at 4 °C and processed within two weeks. The contents of the recovered bottles were concentrated, freeze-dried prior recording of the dry mass by weighing. BSi, total N and P were estimated on dry material.

Water column Si, N and P profiles
For the period 2009-2017, chemical data for the water column Si, N and P concentrations were provided by the Regional Agency for Environmental Protection of the Lombardy Region (ARPA Lombardy) ( Fig. 1). The following variables were selected: N-NO 3 − , N-NH 4 + , TN, SRP, TP and DSi all reported as µg L −1 ; data on total silica were not available. The total inorganic nitrogen (DIN) was calculated as DIN = N-NH 4 + + N-NO 3 − . The sampling frequency for the entire water column (0-256 m) ranged from four per year (one per season, six depths) from 2009 to 2011 to twelve per year (one per month, seven depths) from 2012 to 2017.
Dissolved inorganic nutrient fluxes across sediment− water interface at the deepest lake site Six sediment cores were sampled on 26 October 2018 at 256 m depth to analyse nutrient fluxes across the water sediment interface. The cores (Plexiglas tubes with a 5 cm inner diameter and a height of 30 cm) were collected using a gravity corer (Uwitec core sampler) connected to an electronic winch. Once collected, the cores were immediately submerged in anoxic water retrieved from the same depth and closed with a rubber stopper in order to minimize the exposure to atmospheric oxygen. They were then put in a refrigerated box and transported to the laboratory within 4 h of collection. In the laboratory, sediment cores were housed in a thermostatic room, submerged in an incubation tank (75 L) containing water from the monimolimnion. Incubation was carried at the same ambient temperature (± 1 °C). Anoxia in the tank was guaranteed by covering the incubation tank with a plastic bag and by bubbling N 2 into the water (Nizzoli et al. 2018). Cores were subsequently left to stabilize overnight (~ 12 h). The following day, cores were incubated in dark conditions to simulate ambient light in the deep monimolimnion. Incubation took place for 64 h. During maintenance and incubation, the water inside the cores was gently stirred using magnetic stirrers suspended within each core and driven with a magnet rotated at 30 rpm by an external motor. At the beginning of incubation, the water level in the tank was lowered below the core edges and the upper opening of the cores was closed with a floating Plexiglas lid. Water samples were collected through a valve in the lid using plastic syringes and filtered through Whatman GF/F glass fiber filters at the beginning of the incubation and after 18, 41 and 64 h. Samples were collected in polyethylene vials to determine DSi, N-NO x − and N-NH 4 + , while samples for SRP analysis were preserved in glass tubes. Samples were stored at 4 °C and analysed within 24 h.
Analytical methods N-NH 4 + (Koroleff 1970), N-NO x − (APHA 1998), SRP (Valderrama 1981) and DSi (Golterman et al. 1978) in river water and core incubations were determined using standard spectrophotometric methods (Perkin Elmer, Lambda 35). An analytical blank undergoing the same procedure as the samples, including filtration, was always analysed to correct for sample contamination.
Filters with particulate matter were oven-dried at 70 °C for 24 h. Subsequently, they were analysed to quantify particulate N and P following an alkaline digestion (Maher et al. 2002). Filtered samples were processed with the same digestion to determine total dissolved N and P.
The BSi content of suspended particulate matter and from sedimentation trap samples was analysed following DeMaster (1981): 30 mg of dry sediment was digested with 30 ml of 0.1 M Na 2 CO 3 in polypropylene bottles for 5 h at 85 °C. Subsamples were collected after 3, 4 and 5 h and analysed for DSi using the molybdate blue spectrophotometric method (Golterman et al. 1978). Before analysis, each sample (1 ml of extraction solution) was neutralized with 9 ml of HCl 0.021 M. To correct for Si resulting from mineral dissolution, the Si content of the subsamples was plotted against dissolution time. The y-axis intercept of the linear regression line represents the estimated BSi content. All BSi dissolved within the first 3 h.
Total N and P content in samples of sedimentation traps was determined at the Leibniz Institute of Freshwater Ecology and Inland Fisheries (IGB): N was determined using a CHNS elemental analyzer (Fa. Elementar, Hanau), while total P was determined as SRP after digestion of approximately 5 mg dried material in 2 ml 10 M H 2 SO 4 and 30% H 2 O 2 at 150 °C (see Kleeberg et al. 2010).

Nutrient load, sedimentation rates and sediment water flux estimation
The annual nutrient loadings into and out of the lake were calculated as the product of the discharge weighted mean concentration by the mean annual discharge calculated over the three years (Quilbé et al. 2006) as follows (1): where L is the annual loading (t year −1 ), Ci is the concentration at day i (mg m −3 ), Qi is the mean daily discharge at day i (m 3 s −1 ), Qm is the mean annual discharge (m 3 s −1 ) of the sampling period and k is the conversion factor (31.54 × 10 -3 ) from mg s −1 to t year −1 . Given the absence of data on nutrient concentration and discharge in the minor tributaries, we assumed the same concentrations for the river Oglio in order to provide an estimation of the overall load flowing into the lake.
Daily average total water inflows and outflows were provided by the Oglio Consortium, while data on the Industrial Channel discharge were provided by the Paravisio hydro-electrical company. The total daily inflow to the lake is quantified by the Oglio Consortium using a mass balance approach calculated as the difference between the net daily change of water volume in the lake (quantified from water level changes) and the discharge at the outlet.
Lake nutrient retention was estimated as the difference between the nutrient loading flowing into the lake and the nutrient load exported through the lake outlet. The average retention efficiency over the three years relative to the input was calculated as (2) where L in and L out are the average Si, N or P loads during the three years (t year −1 ) in inlet and outlet waters, respectively, and estimated according to (1).
Atmospheric N input (N Dep ) was estimated from wet and dry deposition data for reactive N (N r ) from the Cooperative Programme for Monitoring and Evaluation of the Long-range Transmission of Air Pollutants in Europe (EMEP 2010, http:// www. emep. int/). In the Lake Iseo area, such N Dep is 950 kg N km −2 year −1 . Unfortunately, no P and Si deposition has been estimated, as data for atmospheric deposition of these nutrients are rarely gathered. However, the few estimates of TP (20 kg km −2 year −1 ) and TSi (111 kg Si km −2 year −1 ) deposition suggest that they are negligible when compared with other inputs to lake Iseo (Hofmann et al. 2002;Kopáček et al. 2011).
Sedimentation rates (mg m −2 day −1 ) were calculated by multiplying the areal dry mass flux rate with the nutrient content in settled material for each time interval. Hourly fluxes of nutrients (mg m −2 h −1 ) were calculated with linear regression as the slope of concentration versus incubation time and then multiplied by 24 h to obtain daily fluxes (mg m −2 day −1 ).

Statistical analysis
Nutrient concentrations in rivers and in-lake sedimentation rates were both divided into two periods: a spring-summer period from April to September (summer) and an autumn-winter one from October to March (winter) which encompass the periods of winter mixing and summer stratification of the epilimnion.
River water nutrient concentrations and stoichiometry were analysed with two-way generalized least square (GLS) models to assess differences between sites (inlets and outlet) and seasons (summer and winter). Nutrient concentrations in particulate matter collected from sedimentation traps, together with their stoichiometry, were also analysed with two-way generalized least square (GLS) models to assess differences between depths (20 and 90 m) and seasons. To deal with observed heteroscedasticity, the argument "weights" was used within the function gls(), and the function varIdent() was implemented to specify variance models. A posteriori comparison of the means was performed using a post hoc Tukey test. Statistical analysis was conducted using the "nlme" (Pinheiro et al. 2016) and "emmeans" (Lenth et al. 2018) packages in the R software v. 3.6.0 (The R Core Team 2019). Nutrient ratios were log transformed prior to analysis as suggested by Isles (2020). The relation of nutrient concentrations to the discharges in the sampling sites (inlets and outlet) was further investigated using the Pearson correlation coefficient. Finally, individual regressions were run between the DSi, SRP and DIN concentrations in monimolimnion and the time (year). The regression slopes and the correlation coefficients were used to estimate annual accumulation (µg L −1 year −1 ) of dissolved nutrients in the deeper layer of the lake. If not otherwise specified, all results in the text, figures and tables are presented as mean ± standard error.

Nutrient concentrations and stoichiometry in inlet and outlet waters
Over the three years, average daily discharge rates underwent natural variability in the river Oglio tributary (annual mean = 24.6 ± 2.9 m 3 s −1 ; min = 10.8 m 3 s −1 ; max = 103.8 m 3 s −1 ), while they were regulated in the Industrial Channel feeding the hydroelectric power plant (mean = 27.2 ± 2.3 m 3 s −1 ; min = 5.4 m 3 s −1 ; max = 50.4 m 3 s −1 ). Water discharge was also regulated at the lake outflow (annual mean = 47.3 ± 4.0 m 3 s −1 ; min = 23.4 m 3 s −1 ; max = 110.8 m 3 s −1 ). Flow rates increased in summer to supply the demand for irrigation and power generation, while they were low in late winter to allow for water storage to be used in the following months (Supplemental Info Fig. S2; Table S1).
TSi, TN and TP concentrations resulted significantly different among the three rivers with TSi also significantly lower in summer compared to winter (Supplemental Info Table S2; Fig. 2). On average, TSi and TP concentrations were not statistically different between the two inlets (Supplemental Info  Table S2; Fig. 2) and were mainly affected by river discharge, particularly in the river Oglio (Supplemental Info Figs. S2, S3 and S4). In the two inlets, particulate pools of all nutrients were positively correlated with discharge (BSi r = 0.69; PN r = 0.67; PP r = 0.61; p < 0.001 for all), whereas DIN (r = 0.50, p < 0.01) and DSi (r = 0.59, p < 0.01) were negatively correlated (Supplemental Info Fig. S5). During flood events, BSi, PN and PP peaks were detected, reaching up to 2181 μg Si L −1 , 350 μg P L −1 and 1172 μg N L −1 .
TSi, TN and TP concentrations were significantly lower in the outlet compared to the inlets (Supplemental Info Tables S1, S3; Fig. 2). The observed changes of TSi and TN coincided with the decrease of the dissolved forms, which accounted for more than 70% of TSi and TN pools, while the differences between the two inlets and the outlet were not statistically different for BSi and PN concentrations (Supplemental Info Tables S1, S3; Fig. 2). On the other hand, TP concentrations were much lower in the outlet than in the inlets (Supplemental Info Tables S1, S3) due to the significant decrease of both SRP and PP, the latter accounting for more than 50% of TP (Fig. 2).
In the outlet waters, TSi concentrations showed a clear seasonality, being higher in autumn-winter than in summer-spring ( Fig. 2; Supplemental Info Table S1), with most of the observed differences coinciding with the significant decrease of DSi concentrations in summer. In contrast, TP and TN concentrations were similar in both seasons ( Fig. 2; Supplemental Info Table S1).
Total nutrient ratios differed between inlets and outlet, but the differences depended on the season and the specific nutrient ( Fig. 3; Supplemental Info Table S2). Overall, changes in nutrient ratios highlighted an in-lake retention of Si relative to N and P, in particular during summer (Fig. 3). During the winter period TP:TSi did not change between inlet and outlet waters, while a significant increase (from 0.015 ± 0.003 to 0.045 ± 0.011) was observed in summer (Supplemental Info Table S3). The TN:TSi ratio significantly increased (Supplemental Info Table S3) from 1.38 ± 0.06 in the inlet to 4.57 ± 0.77 in the outlet. In the outlet, the TN:TSi ratio was higher in summer compared to winter. The TN:TP ratio increased from 84 ± 5 in the inlet to 116 ± 10 in the outlet (Fig. 3), but no differences were observed between seasons.
Both DIN:DSi and DIN:SRP were higher in the outlet compared to the inlets, in addition to PN:PP, with no significant differences detected between seasons (Fig. 3). The mean PP:BSi ratio (0.17 ± 0.02) was one order of magnitude greater than SRP:DSi (0.013 ± 0.004) in both inlet and outlet waters. In inlet waters, the PP:BSi ratio increased up to 0.20 ± 0.03, especially in autumn and winter during high flow events (Fig. 3).
Nutrient loading to and nutrient export from Lake Iseo The mean annual loadings to and exports from the lake are reported in Table 2. Lake Iseo received a mean annual TSi loading of 3044 ± 303 metric tons Si year −1 and exported 747 ± 164 metric tons Si year −1 , equivalent to a 75% net Si retention. The decrease was mainly due to the DSi decrease, which accounted for 86% of TSi and, to a lesser extent, to the BSi retention. Total P exhibited a similar pattern to Si, passing from 126 ± 22 metric tons P year −1 in the inlet to 26 ± 3 metric tons P year −1 in the outlet, equivalent to a 79% in-lake retention, almost equally accounted for SRP (80%) and PP (84%) uptake. In contrast, although decreasing from 2174 ± 139 metric tons N year −1 (inlet) to 1197 ± 71 metric tons N year −1 (outlet), much less of the TN load was retained (45%) compared to TSi and TP. The dissolved inorganic pool was the main component of the TN pool, with DIN composed of more than 90% nitrates, while the N-NH 4 + contribution was almost negligible in both inlets and outlet waters.

In-lake nutrient concentrations and budgets
Since the winter of 2006, there has been no complete water overturn in Lake Iseo and a stable and persistent stratification has been established. Under these conditions, nutrient concentrations differed between water layers (Figs. 4, 5).
In the period 2009-2017, the mean TN concentration in the whole lake ranged from 386 to 1390 µg N L −1 and did not exhibit a clear vertical profile. On average, TN was only slightly higher in the epilimnion (750 ± 24 µg N L −1 ) compared to the monimolimnion (521 ± 28 µg N L −1 ) (Fig. 5). Average DIN concentration was also slightly higher (579 ± 15 µg N L −1 ) in the epilimnion compared to the monimolimnion (388 ± 19 µg N L −1 ); DIN concentration varied over time, mainly due to DIN depletion, especially in the monimolimnion and, to a lesser extent, in the hypolimnion. In the monimolimnion, the average DIN concentration decreased from 2009 (769 ± 30 µg N L −1 ) to 2014 (219 ± 54 µg N L −1 ) (− 117 ± 10 µg N L −1 year −1 , Fig. 4, Supplemental Info Table S4), after which it remained nearly constant. A progressive, but slighter, decrease was also detected in the hypolimnion (Fig. 4, Supplemental Info Table S4). The composition of DIN also varied with depth with N-NO 3 − contributing up to 97% of DIN in the epilimnion and no more than 50% in the monimolimnion, where N-NH 4 + contribution progressively increased from 2009 to 2017.
The dissolved N:Si:P stoichiometry showed a depletion of both DSi and SRP compared to DIN in the epilimnion, whereas both DSi and SRP were stored in the monimolimnion in a larger amount compared to DIN (Fig. 6). The average DIN:DSi molar ratio was higher in the epilimnion (1.2 ± 0.1) compared to the monimolimnion (0.2 ± 0.0). Similarly, the DIN:SRP ratio decreased from 531 ± 21 in the epilimnion to 8 ± 1 in the monimolimnion. In contrast, the SRP:DSi molar ratio increased from 2.7 × 10 -3 ± 0.2 × 10 -3 in the epilimnion to 35 × 10 -3 ± 2 × 10 -3 in the monimolimnion. From 2009 to 2017, the DIN:DSi ratio progressively decreased in the monimolimnion from 0.58 to 0.10 (Supplemental Info Table 4S; Fig. S7), with the same trend found for the DIN:SRP ratio, which decreased from 16 to 4 (Supplemental Info Table S4; Fig. S7). The SRP:DSi ratio, however, remained stable in the monitored period (0.034 ± 0.004), confirming the similar fate that connects P and Si cycling (Supplemental Info Table S4; Fig. S7).

Si, N and P sedimentation rates and fluxes across the water sediment interface
In open waters (90 m deep), the average sedimentation rates were 113 ± 13 mg Si m −2 day −1 for BSi; 17 ± 2 mg N m −2 day −1 for PN; and 2.0 ± 0.2 mg P m −2 day −1 for PP. In the shallow coastal zone (< 20 m deep), sedimentation rates of BSi and PP fluxes were comparable to those in the deeper zone, whereas PN fluxes were only slightly higher (GLS, F 13,1 = 1.997, p = 0.181) ( Table 3).
In the deep sediments, fluxes of DSi and SRP were from the sediment to the water column. A net DIN release was also revealed for N-NH 4 + efflux, which was partially counterbalanced by  Table 3). The N, P and Si stoichiometry of benthic fluxes highlighted an excess of DIN and SRP compared with DSi (SRP:DSi = 0.04 ± 0.02 and DIN:DSi = 1.38 ± 0.11), with both P:Si and N:Si flux ratios greater than those of settled matter (Fig. 7).
Benthic fluxes and inherent stoichiometry of Si, N and P in the deep zone were compared with those of the shallow littoral zone reported in Scibona et al. (2019) (Table 3). In littoral sediments, SRP fluxes were negligible (0.04 ± 0.33 mg P m −2 day −1 ), whereas a N-NO x − uptake greater than a N-NH 4 + release induced a net DIN uptake. Overall, DIN was taken up in the littoral zone and recycled into the water column in the deep sediment, whereas the benthic system was predominantly a net source of DSi to the water column in both littoral and deep sediment.

Discussion
The results of this study demonstrate that Lake Iseo, a deep meromictic lake, acts as a biogeochemical filter by trapping a fraction of the external silica, phosphorus and nitrogen loads (Fig. 8). This is not unexpected, as many other studies have shown that biogeochemical processes occurring in lakes can regulate nutrient transport in the hydrographic network (see among others Seitzinger et al. 2006;Harrison et al. 2009Harrison et al. , 2012Finlay et al. 2013;Verburg et al. 2013;Frings et al. 2014;Maranger et al. 2018;Nizzoli et al. 2018 and references therein). The study described here is unique in analysing Si, N and P simultaneously, revealing that the three nutrients are differently processed and retained as they move through the lake ecosystem. Fractions of the external TSi (75%) and TP (79%) loads retained within the lake were similar and were almost twice as high as TN load (45%). Their stoichiometry was impacted by this difference in element retention, with TN:TSi and TN:TP molar ratios being greater in the outflowing water relative to the tributaries, an imbalance allowing Si to be considered as potentially limiting for siliceous algae (Billen and Garnier 2007). These results are in line with the few studies that had previously addressed large-scale effect of lakes on nutrient stoichiometry (Maranger et al. 2018;Goyette et al. 2019) or had directly quantified the relative retention of N and P in lakes (Verburg et al. 2013). Findings are also consistent with Cook et al. (2010) who reported a significant increase of N:P and N:Si compared to external loads in a shallow tropical lake. In the present study, the observed higher retention of Si and P relative to N could be a consequence of the different composition of the external loads and the biogeochemical pathways of the three elements (see below).

Main pathways governing Si, P and N retention
Nutrient retention is the combined result of four main processes: (a) sedimentation and burial of particulate nutrients delivered by inflowing waters; (b) assimilation of dissolved nutrients into living biomass and subsequent export of particulate forms from the photic zone through sedimentation and burial; (c) net accumulation of dissolved nutrients in the water column and (d) nitrogen export following conversion of reactive forms to unreactive gaseous N 2 (i.e.   (Dupas et al. 2015). Mean values and standard error are reported Fig. 8 Main fluxes and pathways of Si, N and P as dissolved inorganic, dissolved organic and particulate pool in Lake Iseo. In addition to inlet and outlet nutrient loads, the following data are reported: particulate nutrient deposition on sedimentation trap at 20 m and 90 m depth; nutrient regeneration from littoral zone (Reg. 20 m total surface 52.4 km 2 ); and deep sediment (Reg. 256 m total surface 36.7 km 2 ). All fluxes are reported as mg m −2 day −1 with standard error in brackets. Arrow size is proportional to areal flux rates denitrification). In meromictic lakes, in particular, the monimolimnion can play a major role as storage of dissolved Si and P (Lehmann et al. 2015;Viaroli et al. 2018). Overall, differences between input and output loadings of Si were relatively well matched with Si sedimentation rates. The difference between the total Si retained in the lake, calculated as mass budget (103 mg m −2 day −1 ), and the total amount of Si exported from the upper mixed layer, measured with sedimentation traps (115 mg m −2 day −1 ), is about 12% (12 mg m −2 day −1 ), which could be due to some unaccounted deposition and DSi recycling in the littoral zone which seemed to act as a DSi source for the epilimnion (Table 3 and Fig. 8; Scibona et al. 2019). Conversely, P and N mass budgets exhibited relevant differences between sedimentation and net retention: the P sedimentation lacked approximately 2.5 mg m −2 day −1 (~ 56% of the retained P), while the N sedimentation lacked approximately 57% of the retained load (25 mg m −2 day −1 ).
With regard to mass balance studies, there are various uncertainties, depending on the method used to calculate loads (Quilbé et al. 2006), the sampling frequency of the monitoring data given the variability of concentrations and their dependency on hydrology, and the unknown nutrient sink and sources. The missing P load can be explained by the prompt sedimentation of particulate P at the lake tributary inlet. The hypothesis that there is a higher proportional sedimentation of TP with respect to TSi and TN is supported by the relatively high particulate to total P ratio (72%) compared to Si (14%) and N (10%) in external loads. The amount of PP retained by the lake (3.4 mg m −2 day −1 ) was higher compared to the average P sedimentation recovered by traps (1.9 mg m −2 day −1 ) and only slightly lower compared to the PP load flowing into the lake (4.1 mg m −2 day −1 ). In contrast, BSi (14 mg m −2 day −1 ) and PN (6 mg m −2 day −1 ) retained by the lake were lower compared to the average BSi (115 mg m −2 day −1 ) and PN (19 mg m −2 day −1 ) sedimentation rates and only slightly lower compared to the BSi (19 mg m −2 day −1 ) and PN (11 mg m −2 day −1 ) load flowing into the lake. In all likelihood, such differences depend on early deposition of a consistent fraction of the particulate phosphorus delivered to the lake in the proximity of the mouths of the tributaries. This deposition, close to the lake inlets, was probably not collected by the traps which were deployed in open waters, several km from the mouths of the tributaries, so probably only represented the sedimentation occurring in the pelagic portion of the lake. Moreover, the particulate matter recovered from the traps was qualitatively different from that delivered by the tributaries: the settled matter was relatively richer in N and Si (molar Si:P in the range 41-83; and N:P in the range 15-28) and closer to a balanced ratio (Dupas et al. 2015), whereas the external particulate load was richer in P (molar Si:P in the range 0.8-12.7; and molar N:P in the range 3.0-7.9). It can, therefore, be inferred that most of the P retention depended on the sedimentation of allochthonous particulate phosphorus, while a major fraction of the settled particulate Si and N was autochthonous. Consequently, the amount and stoichiometry of dissolved nutrients exported from the lake was influenced by the seasonal evolution of phytoplankton production in the water column. In detail, DSi and DIN concentration decreased from winter to summer and the DIN:DSi ratio increased, suggesting DSi was potentially limiting primary production, especially in early spring. This could be the result of dissolved nutrient assimilation and conversion into phytoplankton biomass and of a shift of primary producer assemblages from diatom-dominated communities in early spring to ones dominated by cyanobacteria in summer (Havens 2008;Salmaso et al. 2018;Leoni et al. 2019).
On the other hand, an explanation of the missing N load could be that up to 70% of the N load is accounted for by nitrates, which can be converted into non-reactive N 2 gas by microbial denitrification under anoxic conditions (Bruesewitz et al. 2012;Nizzoli et al. 2018Nizzoli et al. , 2020. For example, in the nearby Lake Idro, denitrification was found to account for ~ 30% of the retained N load (Nizzoli et al. 2018). Assuming that all the missing N in Lake Iseo can be accounted for by denitrification, areal rates resulted in 22 mg m −2 day −1 (8049 kg km −2 year −1 ), a figure which falls within the range of values detected in other lakes (1.9 -84 mg m −2 day −1 ; Nizzoli et al. 2018 and references therein). This could, however, be an overestimation of denitrification since other processes are likely to contribute to N retention-e.g. unaccounted deposition, or assimilation by macrophytes and microphytobenthos in the littoral zone. Nevertheless, even if overestimated, the contribution of in-lake denitrification to N retention appears lower compared to values found for other lakes (62-100%; David et al. 2006), which, in all probability, depends on the meromictic conditions of Lake Iseo (see below).
Overall, the mass balance approach followed in this study suggests that the main pathways governing the retention of Si, N and P differ between the three elements. When the lake is considered as a whole, early sedimentation of allochthonous particulate phosphorus represents the major pathway for P retention, while assimilation by phytoplankton followed by sedimentation represents the major removal pathway of DSi. Nitrogen retention, on the other hand, partly depends on assimilation followed by sedimentation and burial, whereas up to 60% of retained N is released into the atmosphere through denitrification as inert N 2 .

Effect of meromixis on nutrient retention
As a result of climate change, many lakes are experiencing altered mixing patterns, resulting in lower water-overturn frequency and a shift towards holooligomixis and even meromixis (Danis et al. 2004;Valerio et al. 2015;Ficker et al. 2017;Rogora et al. 2018;Salmaso et al. 2018). This represents an important issue because circulation patterns play a significant role in regulating chemical conditions and primary production in lakes (Wilhelm and Adrian 2008; Salmaso et al. 2014Salmaso et al. , 2018. As a consequence of the permanent and stable stratification of the water column, in meromictic lakes, the monimolimnion can evolve into a permanent nutrient storage, preventing water overturn and nutrient return to the photic zone, ultimately affecting nutrient retention and, therefore, the trophic status of the lake (Nizzoli et al. 2018;Salmaso et al. 2018;Viaroli et al. 2018;Leoni et al. 2019;Lau et al. 2020).
In the specific case of Lake Iseo, the effect of the meromictic condition on nutrient retention capacity can be evaluated by comparing the Si, N and P retention efficiencies measured in this study with theoretical retention efficiencies calculated using widely used published equations (see caption to Fig. 9). With regard to the lake surface (60.8 km 2 ), areal Si, N and P retention was 37.7, 16.1 and 1.6 g m −2 year −1 , respectively; such values are in the upper range of those reported for other lakes worldwide and derived either from model computation or direct measurements (Reynolds and Davies 2001;Cook et al. 2010;Harrison et al. 2012;Verburg et al. 2013;Frings et al. 2014). It should be noted, however, that the fraction of incoming nutrient loads retained in lakes is variable and becomes greater with increasing theoretical water residence times (Seitzinger et al. 2006;Finlay et al. 2013;Verburg et al. 2013;Frings et al. 2014). In Lake Iseo, with a theoretical water residence time of 4.5 years, the predicted retention of TSi (61%) and TP (68%) is actually lower than the retention calculated in this study (75% and 79%, respectively). In contrast, TN retention (45%) is lower compared to that estimated with the use of reference predictive models (53-71%; Seitzinger et al. 2006;Finlay et al. 2013;Rissanen et al. 2013), as already found in the meromictic Lake Idro (Nizzoli et al. 2018) but higher compared to model prediction by Harrison et al. (2009). This comparison, although based on several assumptions, clearly shows that the Si and P retention in Lake Iseo is higher than predicted values, suggesting that meromixis impairs nutrient recycling when compared to mixed lakes. In the light of this analysis, meromixis can be said to increase the capacity of lakes to filter Si and P while, at the same time, likely diminishing the retention capacity of N. Concerning the unbalanced retention of Si and P compared to N, this result can be attributed to the different effect of meromixis on the biogeochemistry of these elements: P and Si have sedimentary cycles (Ruttenberg 2003;Teodoru et al. 2006), whereas the N cycle is dominated by dissolved reactive forms, especially nitrates, that are easily soluble and act as a substrate for denitrification (Seitzinger et al. 2006). Under meromictic conditions, the complete overturn is either much less frequent or does not occur at all, thus causing deep-water mass isolation, oxygen depletion and anoxia (Lehmann et al. 2015;Zadereev et al. 2017;Viaroli et al. 2018). Under these conditions, dissolved nutrients recycled by microbial mineralization of the settled organic matter in deep-water sediments are prevented from being periodically redistributed in the whole water mass, thus continually accumulating in the monimolimnion (Lau et al. 2020). Results from core incubations evidence that sediment in the monimolimnion is a net source of dissolved inorganic silica and phosphorus. Accordingly, from 2009 to 2017, the average net accumulation of DSi (1089 ± 17 metric tons Si year −1 ) and SRP (17 ± 1 metric tons P year −1 ) in the monimolimnion was estimated equivalent to an areal regeneration of 81 ± 15 mg m −2 day −1 DSi and 1.2 ± 0.4 mg m −2 day −1 SRP, which was in the same order of magnitude of benthic fluxes. An approximate evaluation of the fate of reactive Si and P was carried out by comparing internal pathways (sedimentation and regeneration) with external loads, aware that outcomes are affected by both spatial and temporal heterogeneity. Such a comparison suggested that 50-70% of settled Si and 60-100% of settled P were recycled, working their way back to the water column, where they accumulated as dissolved reactive forms. Other studies, such as that by Lau et al. (2020), had already noted that sedimentation is the major P vector in Lake Iseo and that P is quickly released from the sediment to the monimolimnion. Results from the present study highlight the fact that Si sedimentation followed by efficient recycling and accumulation in the monimolimnion also represents an important retention mechanism, accounting for an internal loading which is equivalent to 40-60% of the external annual Si load.
In contrast, despite evidence from core incubations demonstrating that the deep benthic system is a net source of ammonium, the dissolved inorganic nitrogen bulk in the monimolimnion was found to decrease at a rate of 420 ± 37 metric tons year −1 from 2009 to 2014, with an areal consumption of 31 mg m −2 day −1 , after which it remained relatively stable. The onset of meromictic conditions and anoxia, which extended up to 45% of the water volume in Lake Iseo in 2016-2018 (Lau et al. 2020), may have favoured denitrification and the consumption of the nitrate reserve in the monimolimnionthis would explain the observed net DIN decrease. Prolonged water stratification, however, represents a threat to N removal: although reducing conditions stimulate NO 3 − depletion, acting concurrently they can impair nitrification, thus only allowing N recycling through ammonification (Nizzoli et al. 2018).
Dissolved nutrient stoichiometry in the monimolimnion can also be altered by the difference in Si and P recycling and pathways compared to N. Stoichiometry along the depth profile changes from Si and P shortage relative to N in the epilimnion to a Si and P excess in the monimolimnion. On the whole, Si and P accumulated over time in the monimolimnion, due to the lack of complete water mixing; this nutrient bulk, however, is not available in the photic zone on account of the steady stratification. Nevertheless, it is possible for complete overturn to occur as a result of strong perturbations, e.g. prolonged cold and windy weather conditions in late winter like those which occurred in 2005 and 2006, making large quantities of nutrients available in the photic zone (up to 75 µg P L −1 ) and moving anoxic waters towards the surface layers (Holzner et al. 2009;Salmaso et al. 2014;Rogora et al. 2018).

Conclusions
The current research is one of the few studies simultaneously exploring Si, N and P pathways and fate in a deep meromictic lake, thus allowing for an estimation of how and to what extent meromictic lakes act as nutrient filters in the hydrographic network. Results from this study demonstrate that Si and P retention is similar, with both being greater than N retention. In the light of this outcome, lakes similar to Lake Iseo could potentially alter the relative amount of nutrients along the downstream hydrographic network, with nitrogen (nitrate) enrichment compared to dissolved silica and phosphorus.
The Si and P retention capacity is likely to be exacerbated by meromixis: Si and P dynamics are mainly driven by the sedimentation of particulate forms that are then efficiently recycled through the water column as dissolved forms, but accumulate in the monimolimnion and cannot return to the photic zone until stratification persists. On the other hand, prolonged water stratification caused by meromixis could represent a threat for the N removal capacity of lakes, limiting NO 3 − removal. Extreme reducing conditions impair nitrification and, in turn, denitrification due to a shortage of NO 3 − ; as a result of this, ammonium recycled from sediments accumulates in the monimolimnion.
In conclusion, the meromixis enhances Si and P retention in the deep-water mass, which is isolated from the epilimnion. The photic zone, therefore, can become Si and P depleted, undergoing apparent mesotrophic, or even oligotrophic, conditions. This could represent a disruptive scenario if extreme climate events occur favouring water circulation. At the same time, a nutrient disequilibrium in the river and lentic waters downstream can be caused by the Si and P retention vs delivery of reactive N, potentially producing unfavourable conditions for siliceous phytoplankton species, while stimulating flagellates and cyanobacteria. This is a relevant issue, since with climate change, deep lakes are turning from holooligomictic to meromictic, with these conditions likely to become more severe in the future if they are exacerbated by concurrent eutrophication. Further research is, therefore, needed in order to provide new insights with which to understand the role of lakes in regulating nutrient availability and fluxes-these studies would serve as a guide for the maintenance of water quality and the ecological conservation of water systems.