Vertical distribution of picocyanobacteria in deep lakes: the in�uence of inorganic turbidity

Picocyanobacteria (Pcy) represent the dominant photosynthetic fraction in aquatic systems, contributing signi�cantly to global primary production and playing a key role in global biogeochemical cycles. Based on a 20-years dataset of in situ observations in four deep Andean North-Patagonian lakes, we analyzed and presented a simple model to understand how the input of inorganic particles affects light penetration and in�uences the vertical distribution of freshwater Pcy during summer strati�cation. The analyzed temporal series includes two important events (volcanic eruption and glacial recession) that substantially affected lake turbidity. Thus, our mechanistic model was constructed as a function of changes in light extinction coe�cient (Kd PAR ) and mean irradiance of the mixing layer (I m ). Our modeling approach using Bayesian inference and a continuous non-monotonic function successfully predicted changes in Pcy vertical distribution. The obtained model was successful in �tting data of different minerogenic particles (volcanic ashes and glacial clay) and in predicting changes under sharp increases in turbidity (volcanic eruptions) as well as in more steady changes (glacial recession). Pcy maximum abundance increased with transparency (lower Kd PAR values) and the amplitude of the vertical pro�le increased with higher I m values. Using our model, we achieved a full prediction of Pcy vertical distribution under different scenarios of lake transparency and lake thermal structures.


Introduction
In aquatic ecosystems, the vertical distribution of phytoplankton is driven by light (Kirk 1994;Falkowski and Raven 2007; Leach et al. 2018).In the upper layers of transparent environments, irradiance may be too high, and photosynthesis will be reduced due to photoinhibition (Platt et al. 1980).Consequently, phototrophs are constrained to intermediate light levels where they frequently develop deep peaks of biomass called Deep Chlorophyll Maxima (DCM) (Coon et al. 1987;Estrada et al. 1993;Gervais et al. 1997).In addition, lake thermal structure is also important since strati cation enhances light supply by decreasing mixing depth (Diehl 2002).In this sense, changes in vertical mixing affect the composition of phytoplankton communities (Diehl et al. 2002).
Picocyanobacteria (Pcy) are key organisms in most pelagic freshwater and marine systems (Callieri 2010;Flombaum et al. 2013;Visintini et al. 2021) and they are often important components of DCMs (Gervais et al. 1997;Camacho González 2006;Callieri et al. 2007).Due to a higher surface-to-volume ratio, Pcy has a competitive advantage over larger-sized phytoplankton in the uptake of nutrients when they are scarce (Danger et al. 2007).Therefore, Pcy contribution to total phytoplankton is more important in oligotrophic habitats, such as open ocean and clear mountain lakes (Callieri et al. 2012).Water warming and reduced nutrients will bene t Pcy over bigger cell size groups, thereby altering carbon export (Flombaum et al. 2013; Martiny et al. 2013).
Water transparency is a fundamental parameter that responds to changes in climate and land use (Rose et al. 2014).Differences in dissolved organic carbon (DOC) and colored dissolved organic matter (CDOM) have been identi ed, in early studies that included a wide range of temperate lakes, as the primary regulators of variation in diffuse attenuation coe cients (Kd) (Morris et al. 1995;Rose et al. 2009).
Particularly, the differential light absorbance determines the dominance or coexistence of green and red Pcy in lakes and seas (Stomp et al. 2007a;Stomp et al. 2007b).However, a high concentration of inorganic suspended particles in the water column increases light scattering reducing transparency (Kirk 1994).The inputs of inorganic particles into aquatic ecosystems are common, i.e., meltwater due to changes in glacier recession rates and volcanic eruptions, that, in turn, alter the input of minerogenic particles into lakes and oceans (Hamme et al. 2010;Slemmons et al. 2013).Glacial particles can play an important role in regulating the attenuation of both ultraviolet radiation (UVR) and photosynthetically active radiation (PAR) in glacially-fed lakes (Modenutti et al. 2000;Gallegos et al. 2008;Rose et al. 2014).
These lakes often exhibit low DOC concentration values (Morris et al. 1995), but they are not highly transparent due to higher light attenuation near glacial in ows (Modenutti et al. 2000; Hylander et al.

2011).
The resulting light attenuation by turbidity due to minerogenic inorganic particles strongly affects the composition and distribution of planktonic communities (Slemmons et al. 2013).In this sense, the vertical location of the DCM can be very sensitive to changes in minerogenic particle inputs (Hamilton et al. 2010;Modenutti et al. 2013a;Bastidas Navarro et al. 2018).In a study including 100 temperate lakes, light and thermal strati cation resulted in important drivers for DCM structure (Leach et al. 2018).The study was based on in situ chlorophyll uorescence and covers a broad range of lake size, maximum depth, transparency, dissolved organic carbon concentrations, and trophic state.However, predictive and mechanistic models considering changes in minerogenic particles are still scarce.These quantitative approaches would be valuable to assess the sensitivity and recovery ability of aquatic organisms facing environmental changes.
Glacier recession and volcanic eruptions affected, in the last decades, the glacial Andean North-Patagonian Lake district (41°S) causing substantial changes in turbidity at both temporal and spatial scales (Hylander et al. 2011;Modenutti et al. 2013b;Bastidas Navarro et al. 2018).While volcanic eruptions have short-term impacts on the ecosystems (Lindenmayer et al. 2010), glacier recession is expected to gradually increase, and so is the transparency of proglacial lakes (Elser et al. 2020).Here, we modeled the vertical distribution of Pcy during lake summer strati cation as a function of changes in water turbidity due to minerogenic particle inputs.We calibrated our model with a temporal series (20 years) of temperature, light and Pcy vertical pro les in four deep Andean Patagonian lakes, and we assessed its predictive accuracy with cross validation methods.We choose Pcy because of its important contribution to total lake production in the DCM of Andean lakes (Callieri et al. 2007).We took advantage of two important events that affected transparency: a-the Puyehue-Cordón Caulle volcanic eruption that heavily affected Andean transparent lakes (Modenutti et al. 2013b) and b-the Mount Tronador glacier recession that includes a glacial lake outburst ood (GLOF) event (Bastidas Navarro et al. 2018).We will evaluate whether the vertical distribution of Pcy displays a non-monotonic function with depth, where the vertical position of the maximum Pcy abundance depends on turbidity (i.e., light attenuation coe cient) while the magnitude of the maximum depends on the mean irradiance in the mixing layer.

Study area and dataset construction
We constructed a dataset using temporal series of Pcy abundance, temperature, and light vertical pro les obtained from eld summer samplings (2001-2022, n = 534) in four deep lakes (more than 100 m depth) of the Andean North-Patagonian glacial lake district: Lakes Espejo, Correntoso, Frías and Mascardi (Fig. 1).These lakes have a monomictic thermal regime and are oligotrophic (Total Phosphorus < 5 µg L − 1 , epilimnetic Chlorophyll a < 1 µg L − 1 , Dissolved Organic Carbon < 1 mg L − 1 ), and very clear with wide euphotic zones (Morris et al. 1995;Modenutti et al. 2013a;Balseiro et al. 2022).Lakes Mascardi and Frías are glacially-fed lakes (Hylander et al. 2011;Schenone et al. 2020).Lake Mascardi has two arms and the Tronador arm receives glacial water from Manso River producing a longitudinal light gradient due to glacial clay (Modenutti et al 2000).The suspended solid input into the lake declined in the last 10 years (Bastidas Navarro et al 2018), due to the rebuilding of a headwater proglacial lake caused by a glacial outburst ood (GLOF) event that occurred in May 2009 (Worni et al. 2012).On the other hand, Lakes Espejo and Correntoso were affected by the eruption of the Puyehue-Cordón Caulle complex in June 2011 (Modenutti et al. 2013b).Lakes Espejo and Correntoso were sampled at a single central point, while lake Mascardi and Frías were sampled in different sampling points (6 and 4, respectively) due to spatial turbidity heterogeneity (Fig. 1).

Light, temperature and Pcy vertical pro les
Light and temperature vertical pro les were obtained from the upper 60 m of the water column using a PUV500B submersible radiometer (Biospherical Instruments™).Diffuse attenuation coe cients of PAR (400-700 nm) (Kd PAR , m − 1 ) were estimated as regression coe cients from light pro les obtained with the radiometer in the eld.
Eq. 1 where I 0 is the irradiance at the surface and Z is the depth (m).
The depth of the mixing layer (Z m ) was determined using 'rLakeAnalyzer' R package on the temperature pro les (Winslow et al. 2019).The PAR mean irradiance of the mixing layer (I m , dimensionless) was calculated according to Sterner et al. (1997).Kd P AR *Z m using an Olympus BH 50 epi uorescence microscope tted with blue excitation (U-MWB) and green excitation (U-MWG) light lters.Phycocyanin-rich (PC-rich) cells were not registered in the samples or observed in very low abundances (less than 1% of total Pcy abundance only in Lake Frías).Cells were counted using an image analysis system (Image ProPlus; Media Cybernetics).
Modeling absolute abundance or biomass would need to include multiple variables such as nutrient availability (Stockner 1991

Model description
We built a non-monotonic model for vertical Pcy abundance based on the probability density function of the Beta distribution: Eq. 3 this equation has the potential to display non-monotonic functions with only two parameters, and (Gupta and Nadarajah 2004).The shape of the function depends on both parameters: the relation between and de nes the position of the maximum, while the magnitude of the parameters de nes the amplitude of the curve.A normalization constant is also present to ensure that the total probability is 1.However, this function alone has a narrow domain (D: [0, 1]), thus, we standardized our response and predictor variable and replaced them in the former equation: Eq. 4 where the response variable 'Pcy' corresponded to the normalized Pcy abundance obtained by dividing Pcy abundance at each sampling depth by the average Pcy abundance of all depths at that sampling point.On the other hand, the predictor variable ' ' corresponds to the depth levels divided by 120 m, which is the maximum depth that light can reach at maximum transparency (Kd PAR ~ 0.03, Sargasso Sea, Kirk (1994)), to standardize this variable to values in a range from 0 to 1.
Finally, we added Kd PAR as the predictor variable affecting parameters α and β and I m affecting the whole formula: Eq. 5 where parameter is the constant from the previous equation and ω is a conversion factor for transforming Kd PAR as a dimensionless variable.The value of ω is 1, and its units are meters (m), which is the inverse of those of Kd PAR .We set Kd PAR dividing and multiplying since we expect that increasing turbidity will increase over and thus the function maximum will be located at lower depths.In addition, we set multiplied by γ as an exponent of the whole equation, as we expect that a higher I m will increase photoinhibition at surface layers, and enhance the maximum at deeper layers, hence increasing the amplitude of the normalized Pcy abundances.Finally, we divided by Kd PAR to ensure that the average of each Pcy vertical pro le is 1.

Model t and predictive accuracy
We t the four parameters of our model ( , , and from Eq. 5) with our dataset with a Bayesian statistical approach performed using STAN code interfaced with R (R Core 2019) through 'brms' package (Bürkner 2017).We ran a sensitivity analysis of our model by changing the prior distributions of model parameters.In addition, we compared the effect of the different types of inorganic turbidity (volcanic ashes vs glacial clay) by tting our model with the dataset divided by the predominant source of inorganic turbidity: lakes Espejo and Correntoso (volcanic ashes) vs Mascardi and Frías (glacial clay).We assessed the predictive accuracy of our model with a K-fold cross validation approach that iteratively split our datasets into 90% training data and 10% test data withheld from model tting (Boyce et al. 2002).We repeated this procedure ten times (K = 10), with no test data repetition between folds (Wenger and Olden 2012).We then validated our model using mean absolute error (MAE) calculated for each of the ten test datasets.

Results
Our dataset showed a wide range of Kd PAR values, ranging from 0.089 m − 1 in the transparent lake Espejo to 2.172 m − 1 in the glacially-fed lake Frías (Fig. 2).Throughout the 20-year sampling period we observed that the two transparent lakes (Espejo and Correntoso) were drastically affected by the Puyehue Cordón Caulle eruption increasing K dPAR values up to ~ 0.3 m − 1 and then values decreased to the pre-eruption ones.Lake Mascardi was analyzed considering turbid sites (sampling stations 1, 2 and 3) and clear sites (sampling stations 4, 5 and 6).Both sites differ in K dPAR and showed a global decreasing trend in the last 15 years following the GLOF event in 2009 (Fig. 2).In addition, the mean irradiance at the mixing layer (I m ) also showed a wide range from 0.112 (lake Mascardi) to 0.709 (lake Correntoso) (data not shown but see Fig. 3).Pcy vertical distribution changed with Kd PAR (Fig. 3a-c) and with I m (Fig. 3d-e).The position of the observed maximum Pcy abundance in the vertical pro les was deeper with higher transparency (lower Kd PAR values, Fig. 3a-c).In addition, the amplitude of the Pcy vertical pro le changed with the mean irradiance at the mixing layer (I m ).When I m was low (high turbidity or deeper mixing layer) the Pcy abundances showed little vertical variation (Fig. 3d, all boxes around the dotted line).However, as the I m increased (low turbidity or shallower mixing depths) the amplitude of the vertical distributions increased (Fig. 3e-f).
Our model showed an excellent t to the dataset as most points were located around the 1:1 observed vs. predicted line (Fig. 4 MAE = 0.42).The sensitivity analysis of the model parameters showed a high overlap between their posterior distributions, regardless of the selected priors (Fig. S1).In addition, when comparing inorganic turbidity sources, we found a high overlap between the posterior distribution of all C α β C γ model parameters tested ( = 37.3%, = 37.9%, = 80.8% and = 41.4%),indicating that tted parameter values of our model were independent of the type of inorganic particles (volcanic ashes and glacial clay) (Fig. 5).
We found an excellent predictive accuracy of our model for the 10-fold partitions of the dataset as we observed small differences between predicted normalized Pcy abundances and the observed ones (Fig. S2, MAE = 0.35 ± 0.12).Finally, based on the obtained model we achieved a full prediction of Pcy vertical distribution under different scenarios of lake transparency (Kd PAR ) and lake thermal structures (mixing layer depth) (Fig. 6).As we expected, the model predicts deeper maximum Pcy abundances as Kd PAR decreases (comparative Fig. 6a-d).In addition, the amplitude of the Pcy vertical distribution was well predicted by the depth of the mixing layer through the I m (Fig. 6, line colors within each graph).
Consequently, lower transparency (higher Kd PAR values) and deeper depth of the mixing layer (lower I m values) dampen Pcy vertical distribution.

Discussion
Here, we analyzed and presented a simple model to understand how the input of inorganic particles affects light penetration and in uences the vertical distribution of freshwater Pcy.In this study, we combined 20-year of eld sampling with modeling to understand the vertical distribution of Pcy during summer lake strati cation.The analyzed temporal series includes four deep lakes and two important events that substantially altered the input of inorganic particles that affected lake turbidity (volcanic eruption and glacial recession).In addition, the location of each lake and their wind exposure along with the temporal variation in wind intensity allowed us to have different mixing layer depths and thus different I m .The obtained model for vertical Pcy distribution was successful in tting data from lakes that differed in their main source of inorganic turbidity (volcanic ashes and glacial clay).Thus, our modeling approach using Bayesian inference and a continuous non-monotonic function successfully predicted changes in Pcy vertical distribution under different turbidity and mixing layer depths.
As in many other approaches analyzing phototroph vertical pro les (Platt et al. 1980 Organic turbidity by CDOM signi cantly reduces the penetration of green and blue wavelengths relative to red ones (Kirk 1994).A previous study in Andean Patagonian lakes indicated a prevalence of blue-green light at DCM in deep lakes, while shallow lakes with high CDOM were dominated by red light (Pérez et al. 2002).These optical differences generate distinct spectral niches for Pcy (Stomp et al. 2004) being the PE-rich cells dominant in clearer environments with blue-green wavelength prevalence (Stomp et al. ). Accordingly, in Andean lakes phytoplankton photosynthetic accessory pigment concentrations and speci c absorption spectra showed more e cient harvesting of the available blue-green light at deep levels of the water column (Pérez et al. 2002).However, inorganic turbidity due to volcanic ashes and glacial clay showed a similar scattering of blue, green and red wavelengths (Modenutti et al. 2013b), explaining why we still found the dominance of PE-rich cells despite increasing inorganic turbidity.
Lake thermal structure and vertical stability also affect Pcy abundance, community composition, and vertical distribution (Crosbie et al. 2003;Callieri 2010).Vertical mixing is important for light because mixing can lead to a shortage of light if planktonic organisms are frequently dragged down to the bottom, whereas strati cation enhances light supply by decreasing mixing depth (Diehl 2002).Andean mountain lakes are exposed to high and variable wind affecting the depth of the mixing layers, hence the Im (Modenutti et al. 2008).Here, we considered these different lake conditions by including I m as a predictor variable in our model.Higher I m will increase the probability of photoinhibition decreasing the Pcy abundance in the upper layers.In our dataset, we observed that the amplitude of Pcy vertical distribution was related to I m .In the model, I m was included as an exponent affecting the amplitude of the nonmonotonic curve.The resulting equation successfully represents the observed pattern and has good predictive accuracy.This was particularly evident when we simulated the Pcy vertical pro les for increasing Kd PAR values and different mixing layer depths (Fig. 6).In this simulation, we were able to highlight how the two variables (Kd PAR and I m ) interact in driving the Pcy vertical distribution.
Due to climate change, the frequency and magnitude of events with high inputs of suspended solids (i.e.glacial recession) are expected to increase (IPCC 2021).In addition, volcanic eruptions have shaped life over geological times through abrupt changes with conspicuous short-term impacts on the ecosystems (Lindenmayer et al. 2010).Here, we showed that turbidity induced by volcanic ashes and glacial clay have similar effects on the vertical distribution of Pcy.In this sense, our model was able to predict changes under sharp increases in turbidity (volcanic eruptions) as well as in more steady changes (glacial recession).In the former, with abrupt changes, Pcy distribution moves upwards and then returns to deeper depths as the lakes recover their transparency.In the second scenario, glaciers are expected to disappear, and in the meantime, GLOFs are expected to increase their frequency (Zemp et al. 2015) and the transparency of proglacial lakes will gradually change (Elser et al. 2020), driving the maximum Pcy abundances more smoothly as predicted by our model.In addition, changes in wind patterns will affect the depth of the mixing layer (Le et al. 2023) and thus light availability in the epilimnion for phototrophs (Diehl 2002).Finally, other climatological events, such as changes in precipitation patterns and subsequent runoff, are important turbidity drivers of aquatic systems (Zhang et al. 2007), and our model has the potential to predict the consequences of these changes.

Declarations
Competing interests   Values of Kd PAR (m -1 ) of sampled lake during summer strati cation from our temporal series.In Lake Mascardi, clear sites correspond to sampling stations 1 to 3 while turbid sites correspond to sampling stations 4 to 6. References: L. = Lake.Posterior probability density functions of our model parameters calibrated with the lakes divided by the predominant source of inorganic turbidity (Glacial Clay: lakes Mascardi and Frías; Volcanic Ashes: lakes Espejo and Correntoso).

;
Eilers and Peeters 1988;Leach et al. 2018, among others) light (Kd PAR ) was one of the main factors triggering the nonmonotonic vertical distribution of Pcy.Different processes of photosynthetic cell metabolism such as photosynthesis and growth are light-dependent(Platt et al. 1980).Pcy appears to be favored by low light conditions(Wehr 1993;Callieri et al. 2007), conforming deep maxima as seen in our dataset and successfully predicted by our model.In addition, underwater light climate is an important factor in niche differentiation for Pcy with different pigment compositions(Stomp et al. 2007a;Holtrop et al. 2021).

Figure 1 Map
Figure 1

Figure 3 Observed
Figure 3

Figure 4 Model
Figure 4

Figure 5
Figure 5 ;Weisse 1993;Callieri et al. 2012;Leach et al. 2018).Therefore, we normalized the Pcy vertical distribution to avoid the noise of absolute abundances.In this way, all our vertical pro les average 1, which allowed us to generalize the vertical Pcy distribution.