How Does Environmental Inter-annual Variability Shape Aquatic Microbial Communities? A 40-Year Annual Record of Sedimentary DNA From a Boreal Lake (Nylandssjön, Sweden)

To assess the sensitivity of lakes to anthropogenically-driven environmental changes (e.g., nutrient supply, climate change), it is necessary to first isolate the effects of between-year variabilit ...

To assess the sensitivity of lakes to anthropogenically-driven environmental changes (e.g., nutrient supply, climate change), it is necessary to first isolate the effects of between-year variability in weather conditions. This variability can strongly impact a lake's biological community especially in boreal and arctic areas where snow phenology play an important role in controlling the input of terrestrial matter to the lake. Identifying the importance of this inherent variability is difficult without time series that span at least several decades. Here, we applied a molecular approach (metabarcoding on eukaryotic 18S rRNA genes and qPCR on cyanobacterial 16S rRNA genes) to sedimentary DNA (sed-DNA) to unravel the annual variability of microbial community in 40 years' sediment record from the boreal lake Nylandssjön which preserve annually-laminated sediments. Our comparison between seasonal meteorological data, sediment inorganic geochemistry (X-ray fluorescence analyses) and organic biomarkers (pyrolysis-gas chromatography/mass spectrometry analyses), demonstrated that inter-annual variability strongly influence the sediment composition in Nylandssjön. Spring temperature, snow and ice phenology (e.g., the percentage of snow loss in spring, the timing of lake ice-off) were identified as important drivers for the inputs of terrestrial material to the lake, and were therefore also important for shaping the aquatic biological community. Main changes were detected in the late-80s/mid-90s and mid-2000s associated with increases in algal productivity, in total richness of the protistan community and in relative abundances of Chlorophyta, Dinophyceae as well as Cyanobacteria abundance. These changes could be linked to a decline in terrestrial inputs to the lake during the snow melt and run-off period, which in turn was driven by warmer winter temperatures. Even if our data shows that meteorological factors do affect the sediment composition and microbial communities, they only explain part of the variability.

INTRODUCTION
Lake ecosystems are sensitive to anthropogenically-driven environmental changes such as modifications in nutrient supply or climate change (Carpenter et al., 1999;Williamson et al., 2009;Moss, 2012). However, it is often challenging to isolate lake ecosystem responses that are the consequence of local factors (i.e., nutrient supply) and regional/global drivers (i.e., climate) using only short-term monitoring data due to the inherent seasonal and between-year variability in factors such as solar irradiance, temperature and precipitation. In boreal and arctic regions, the inter-annual variability in seasonal conditions influences the timing, duration and strength of snow and lake-ice phenology, which in turn control catchment hydrology and the export of terrestrial matter to the lake. Variations in snow cover, ice-off and run-off may also influence directly and indirectly aquatic communities (e.g., timing, length, and intensity of diatom growing seasons and community composition; Rühland et al., 2015;Warner et al., 2018;Maier et al., 2019). Analyzing changes of complex aquatic communities over a relatively long-term period could facilitate separating such inter-annual variability from longer-term changes in anthropogenic pressures and climate.
Among aquatic communities, the microbial community may be particularly useful in teasing apart the responses of lake ecosystems to inter-annual variations and long-term changes. Microbes act at the base of lake trophic networks, as phototrophs, bacterivores, saprotrophs and parasites, and are thus fundamental for ecosystem functioning (Sherr and Sherr, 1988;Caron et al., 2009;Newton et al., 2011). Among them, the diversity and community composition of aquatic protistan communities (i.e., unicellular eukaryotes) have been well studied using DNA sequencing methods (Simon et al., 2015;Filker et al., 2016;Debroas et al., 2017). However, their temporal dynamics have been studied mostly on short timescales, i.e., from days to seasons Mukherjee et al., 2017;Princiotta and Sanders, 2017), and rarely over timescales relevant to investigate their responses to inter-annual variability (except Cram et al., 2015 for marine bacterial communities) or longerterm pressures. The longest DNA-based monitoring time serieswhich extend up to one decade (e.g., BATS, HOT, SPOT and BBMO; Fuhrman et al., 2015;Giner et al., 2018) -have been useful to identify day length (e.g., luminosity) and temperature as major drivers of community turnover. These studies spanning up to 10 years showed that the aquatic community is strongly influenced by factors with a strong seasonal variability. Such work illustrate the necessity to look at even longer timescales in to effectively separate seasonal and inter-annual variability from longer-term drivers.
Application of DNA-based analysis to sedimentary archives can be a powerful tool to extract long-term data on lake aquatic communities (see Domaizon et al., 2017 for review). These archives can reveal the responses of key components of aquatic food webs to environmental changes over longer timeframes (Capo et al., 2016(Capo et al., , 2017aMonchamp et al., 2018;Tse et al., 2018;Vuillemin et al., 2018). For example, sedimentary DNA (sed-DNA) was used to identify the responses of protistan communities over centennial to millennial timescales to longterm forcing factors such as eutrophication and climate change (Capo et al., 2016(Capo et al., , 2017a. In the present study, we investigate annually resolved changes in the microbial community from the boreal lake Nylandssjön over a 40-year period (1973-2012). In a previous study, Capo et al. (2017b) demonstrated that microbial eukaryotic DNA was well preserved in the recent sediment of Nylandssjön despite significant diagenetic changes such as the rapid loss of carbon (Gälman et al., 2009). In combination with inorganic geochemistry (wavelength dispersive X-ray fluorescence spectroscopy; as used in Rydberg and Martínez-Cortizas (2014), biomarkers for terrestrial and algal organic matter [pyrolysis-gas chromatography/mass spectrometry; as described in Tolu et al. (2015)] and meteorological monitoring data (SMHI) 1 the questions we sought to address were: (1) To what extent is the aquatic ecosystem influenced by the interannual variability in environmental conditions (e.g., snow and ice phenology, modifications in the export of terrestrial matter into the lake)? (2) How has the aquatic microbial community, particularly protistan and cyanobacterial communities, changed over this 40-year period? (3) Can we discern superimposed local or global drivers of ecosystem change over this timeframe, such as land use or climate?

Study Site and Sediment Sampling
Nylandssjön (62 • 57 ′ N, 18 • 17 ′ E) is a circum-neutral, mesotrophic, and dimictic lake in which distinct well preserved annual laminations (varves) covered the past c. 100 year (Renberg, 1986, Figure 1). The lake, which is located at 34 m above the current sea level, was formed as a result of isostatic rebound about 3,000 years ago. Lake and catchment areas are 0.28 and 0.95 km 2 , respectively. A number of sediment freeze cores have been collected over the years (1978-2015cf Supplementary Material). During the period 2004-2015, the freeze cores were then manually sub-sampled into annually resolved samples following the dark winter layer using a scalpel under magnifying lamp in a freezer room. The varves from Nylandssjön thus provide an accurate timescale to date each sediment layer over the 1973-2012 time period.

Sed-DNA Analysis
A total of 40 individual varves covering the years 1973-2012 were sub-sampled from a sediment core collected in 2013 for the DNA-based analysis. DNA extraction, PCR amplification, library preparation and bioinformatics were performed as described in Supplementary Material. Briefly, DNA extraction was performed in duplicate for each of the 40 varved samples (i.e., 80 samples) using approximately 500 mg of wet sediment with the NucleoSpin R Soil kit (Macherey-Nagel, Düren, Germany). Total Cyanobacteria abundance was quantified by performing qPCR analysis using a SYBRGreen assay by targeting the 16S rRNA-ITS region. The values obtained were the mean values of the duplicates from each sediment layer and expressed as the abundance of Cyanobacteria (copies/g sed). For the highthroughput sequencing of eukaryotic DNA, the V7 region of 18S rRNA gene was PCR amplified by targeting a 260 bp long fragment, from ∼25 ng of environmental DNA extract for each sample and using the general eukaryotic primers 960f (5 ′ -GGCTTAATTTGACTCAACRCG-3 ′ ) (Gast et al., 2004) and NSR1438 (5 ′ -GGGCATCACAGACCTGTTAT-3 ′ ) ( Van de Peer et al., 2000). Library preparations and paired-end (2 × 250 bp) sequencing (MiSeq Illumina instrument, San Diego, CA, USA) were performed at Fasteris SA (Geneva, Switzerland). The raw data were deposited in FigShare 2 and processed as described in the Supplementary Material. Briefly, the paired-end reads were merged and cleaned (no undefined bases, minimum sequence length of 200 bp, no sequencing error in primers, removing of chimera) using UPARSE and UCHIME (Edgar et al., 2011;Edgar, 2013) tools allowing to obtain a total of 5,847,791 raw sequences. The remaining DNA sequences were clustered into OTUs (i.e., Operational Taxonomic Units) at a 95% similarity threshold. The taxonomic affiliation was performed by BLAST against an SSURef SILVA database (Pruesse et al., 2007) enriched by lacustrine DNA sequences originating from various studies on lacustrine systems. For each sample, OTUs detected in only one of the two duplicates were removed from the analysis. Finally, the number of DNA sequences obtained for each of the 40 varves was standardized at 75,000 resulting in a total number of 12,311 OTUs. The final OTU abundance table was used to calculate the diversity metrics of protistan community [i.e., Hill's numbers: the species richness (q = 0, OTU richness), the exponential form of Shannon entropy (q =1, transformed Shannon index); and the transformation of Gini-Simpson measures (q = 2, transformed Simpson index)]; R package iNEXT (Hsieh et al., 2016).

Geochemical Composition
The geochemical analyses were made in part to extend the time series previously used  by Rydberg and Martínez-Cortizas (2014). Samples from 1953-2002 and new samples for 2004-2014 were (re-)analyzed using wavelength dispersive X-ray fluorescence spectroscopy (WD-XRF). Prior to analysis, freeze-dried samples were weighed to calculate the bulk density (BD) and dry-mass accumulation rate (DMAR). All samples were analyzed for (i) C and N contents, using a Perkin-Elmer 2400 CHNS/O analyzer operated in CHN-mode (cf. Gälman et al., 2006); and (ii) 23 elements, using a Bruker S8 Tiger wavelength dispersive X-ray fluorescence spectrometer (WD-XRF) instrument. Elemental concentrations were calculated using a calibration method suitable to small sample sizes (50 mg) as developed by Rydberg (2014)

(cf. Supplementary Material).
To facilitate comparison with organic biomarkers and DNA data, the geochemical data were subjected to a principal component analysis (PCA geo ) in order to reduce the 25 geochemical variables to a few components that are related to underlying processes controlling the sediment geochemistry ( Figure S1). The geochemical data was standardized (z-scores) prior to PCA. Because this multivariate technique considers the inorganic sediment matrix as a whole rather than the single variables, the resulting PCs are less sensitive to the analytical uncertainties and spatial heterogeneity that might affect the down-core record of a single variable in a specific core.

Organic Biomarkers
Biomarkers of terrestrial and algal organic matter (OM) were determined using pyrolysis-gas chromatography/mass spectrometry (py-GC/MS) for sediment samples retrieved from a core collected in 2006, which covers the period 1973-2005. The py-GC/MS analysis, peak integration and identification follow the method optimized by Tolu et al. (2015), and is described in detail in Supplementary Material. The pyrolytic organic products used in this study include: (i) 16 lignin oligomers and 4 phenolic compounds described as biomarkers for terrestrial OM (Meyers and Ishiwatari, 1993); and (ii) 6 specific Py-products of chlorophyll (i.e., prist-1-ene, phytenes, and phytadienes; Nguyen et al., 2003), 6 specific Py-products of proteins (i.e., 2,5diketopiperazines; Fabbri et al., 2012) and 13 steroids described as biomarkers for algal OM (Ishiwatari et al., 1991;Nguyen et al., 2003;Ninnes et al., 2017;Tolu et al., 2017). Changes in their fluxes to the sediment were estimated by dividing their signal intensity (i.e., peak area) by the dry mass accumulation rate (DM-AR). The parameters are thus referred to as: DM-AR-normalized signal intensity of lignin oligomers, phenolic compounds, chlorophylls, proteins, and steroids. Detailed information about the identified compounds are given in Table S1.

Meteorological Data
Meteorological data (daily air temperature, precipitation, and snow depth) covering the period 1972-2012 were obtained from the Swedish Meteorological and Hydrological Institute 1 for the weather station in Härnösand, 40 km south of Nylandssjön.
For each year, daily air temperature and precipitation values were converted to seasonal resolution (mean and sum values, respectively). The seasons were defined as follows: winter, November-March; spring, April-May; summer, June-August; and autumn, September-October. Daily snow depth data was used to calculate the daily snow gain and loss and the annual snow accumulation (Snow gain ) during the winter-spring period (i.e., 20 October to 10 May). The percentage of spring snow loss (SSL) was determined by dividing the amount of snow loss during the spring period by the snow accumulation for the entire winter-spring period. We also calculated the annual sum of precipitation during snowfall days (Snow prec ) accounting for the quantity of water released during melting events. The ice off timing (tIOT) of Nylandssjön was calculated by applying a model from Weyhenmeyer et al. (2004) using meteorological data (for the whole period) coupled with monitoring data of ice-off timing in Nylandssjön over the 2002-2012 period (cf.

Data Analysis
A change point analysis was performed on univariate and multivariate data using the function e.divisive from the R package ecp (parameters: R = 999, alpha = 1, min.size = 2; James and Matteson, 2013). A Bray-Curtis dissimilarity matrix was built from the OTU relative abundance table and used to assess (i) the temporal variation in community structure via the use of a principle coordinates analysis (PCoA) and (ii) the clustering of samples by hierarchical clustering tree analysis (iii) the relationships between the community structure and environmental parameters. This analysis was performed using the functions wcmdscale (k = 2, eig = T), hclust and env.fit (permutations = 999, type = series) from the vegan R package (Oksanen et al., 2015) was employed based on the goodness-of-fit statistic (defined as r 2 = 1-ssw/sst, where ssw and sst are withingroup and total sums of squares). Bivariate correlations were calculated between environmental and biological parameters as Spearman rank using a 2-tailed test of significance. Because of relatively high autocorrelation in the data, i.e., the observations are not independent, the effective sample size was reduced according to Dawdy and Matalas (1964) prior to calculating the p-values (i.e., adj-p in the text). Adjusted p-values were not determined for correlations with PC1 geo -scores and transformed Shannon index values due to two low effective sample size for these variables. Finally, an analysis using generalized additive models (GAMs) (Hastie and Tibshirani, 1986) was performed to investigate when and to what extent environmental parameters contribute to the changes in biological parameters (structure, richness and diversity of protistan communities, richness and relative abundance of the 11 most abundant protistan groups and the abundance of Cyanobacteria), as detailed by Capo et al. (2017a). The response variables that do not follow Gaussian distribution were log-transformed prior to modeling. The intervals do not correct for multiple comparisons and do not have 95% coverage. The functions gam and predict.gam from the "mgcv" package in R (Wood, 2015) were used and GAM parameterization was performed following technical recommendations from Simpson and Anderson (2009); Simpson (2018) GAMS were fitted using restricted marginal likelihood (REML), a continuous-time first-order autoregressive [CAR(1)] process in the residuals and the parameters select = TRUE and family = Gaussian(link = "identity"). Certain figures were made with ggplot2 R package (Wickham, 2016).

Temporal Dynamics of the Biological Community
The protistan community from the 40-year Nylandssjön sediment record was characterized by a total richness of 12,311 OTUs. The OTU number in a single varve ranged from 1,234 to 3,371 OTUs with the lowest and highest values recorded in 1976 and 2011, respectively (Figure 2; Datasheet S1). Among the 47 taxonomical groups, 11 accounted for more than 55% of the total richness and 97% of the total number of DNA sequences, with Chlorophyta and Ciliophora as major groups (Table 1; Figure S2). The abundance of Cyanobacteria, measured using qPCR analysis, ranged from 4.5 * 10 5 to 3.7 * 10 6 copies per g wet sediment, with the lowest and highest values in 1989 and 2008, respectively (Figure 2; Datasheet S1).

Modifications in Sediment Geochemistry
The first three principal components (PC1-3 geo ) of the PCA made on the geochemical dataset (spanning 1953-2014) explained 75.5% of the total variance ( Figure S1; Datasheet S1). PC1 geo (30.4% of total variance) reflected variations in the relative proportions of mineral matter and organic matter, because it separates variables associated with the heavy mineral fraction (e.g., DM-AR, Si, Al, Na, K, Ti; positive loadings), from elements of organic matter (i.e., C, N, S; negative loadings). PC2 geo (30.1% of total variance) separated samples rich in the fine-grained detrital fraction, i.e., clay minerals (Y, Ca, Zr, Mg, Sr, Ti, Rb, Al; positive loadings; Peinerud et al., 2001), from samples with high biogenic silica content (represented by the Si:Al-ratio; negative loadings). PC3 geo (15% of total variance) was characterized by positive loadings for P, Fe and Mn, indicating a relationship with nutrient supply and redox conditions in the water column. When examining temporal trends in the three main principal components, we identified three change points in the PC1 geo -scores; first in 1977, when the PC1 geo -scores dropped below the long-term average, then in 1993, when the decreasing trend is temporarily halted, and finally in 2006 when the decreasing trend resumes (Figure 4A). Prior to the 1970s, PC1 geo -scores were relatively stable (around +1; Figure 4A), indicating that a higher proportion of the material reaching the sediment was mineral rather than organic material. In the early to mid-1970s, the PC1 geo -scores decreased, indicating a transition to a more organic rich sediment. This decreasing trend continued toward the top of the record, except for two periods, during the late-1980s and early-2000s, when the PC1 geo -scores temporarily increased. Because the Ti-accumulation rate (Ti-AR) closely followed the PC1 geo -scores, while the C-accumulation rate (C-AR) didn't show any increasing trend (Figure S1), the change in sediment composition toward a more organic rich sediment was driven by a decreased input of mineral material rather than an increase in the input of organic material.
Change point analysis identified two shifts in the PC2 geoscores; 1988 and 1999. Prior to the first change point, no clear trend was detected. The PC2 geo -scores then dropped sharply in the mid-1980s and mid-1990s, indicating a higher proportion of biogenic Si relative to detrital fine grained mineral material. These periods with low PC2 geo -scores were followed by a period with high PC2 geo -scores around the year 2000, and then a return to average PC2 geo -scores.
The change point analysis of the third PC (PC3 geo ) identify two shifts, one in 1978 and the second in 2006. The temporal trend of PC3 geo -scores was marked by a distinct shift to higher scores around 1977-1978. At this time, the P concentration increased from on average 3200 µg g −1 for the period prior to 1977, to about 5,100 µg g −1 for the period from 1978 to 2000. According to local farmers, this shift coincided with a shift in farming practices, from a more traditional style of farming relying mainly on manure to a more intense farming using also commercial fertilizers. The PC3 geo -scores then remain high until the early 2000s when they decrease.
In agreement with PC1 geo and PC2 geo scores, change point analysis on the biomarkers of terrestrial OM, i.e., DM-ARnormalized signal intensity of lignin oligomers (specific of vascular plants) and phenolic compounds (in high abundance in vascular and non-vascular plants), indicates a decrease from the early-80s and late-80s ( Figure 4A).

Inter-annual Variability of Meteorological Parameters
Seasonal air temperatures ranged around a mean summer temperature of 14.6 ± 1.0 • C (range: 10-19 • C) to a mean winter temperature of −3.3 ± 2.2 • C (range: −16 to −5 • C) (Figure 4B; Datasheet S1). Winter mean air temperatures exhibit the highest amplitude from the lowest mean temperature of −7.4 • C in 1986 to the warmest mean temperature of 0.4 • C in 1973. Precipitation was similar for all seasons with mean values of 181 ± 101 mm, but higher values generally occurred in winter and autumn (488 and 407 mm, respectively).
For each year, we used snow accumulation (Snow gain ), sum of precipitation during snowfall days (Snow prec ) and percentage of spring snow loss (SSL) to investigate temporal changes in snow phenology for the lake and its catchment (Figure 4B; Figure S5). Some years or time periods (e.g., 1990 to 1993 and 2005) were characterized by low Snow gain , Snow prec and spring snow loss while others (e.g., 1977 and late 1990s) show high Snow gain , Snow prec and spring snow loss values. A change point was detected for spring snow loss values in 1989. Over the 40-year period, the modeled lake ice-off timing (tIOT index) ranged from the Julian day 112 (25 April) to day 130 (8 May) (Figure 4B). Three years were characterized by early ice-off timing within the observational period of study (i.e., 2002(i.e., -2012(i.e., ): 2002(i.e., , 2007(i.e., , and 2011(i.e., . The years 2002(i.e., and 2007 were marked by late snowfalls (see Snow gain values in mid-December 2002 and mid-January 2007) and relatively weak percentage of snow loss during spring. The year 2011 was instead characterized by a long period of negative air temperature values with a high accumulation of snow and strong loss of snow during spring ( Figure S5).

Yearly Environmental Changes in Nylandssjön and Its Catchment
In boreal systems, a warmer winter and spring is characterized by lower snow accumulation (Snow gain ) and lower percentage of spring snow loss (SSL), which results in less intense spring run-off events as compared to cold winters (Renberg, 1981). The meteorological data for Nylandssjön were consistent with this hypothesis, where negative correlations occur between winter temperature and Snow gain (ρ = −0.70, adj-p < 0.001) as well as between winter temperature and spring snow loss (ρ = −0.59, adj-p < 0.001) ( Table S2). In turn, spring snow loss was dependent on both spring temperatures (ρ = −0.41; adj-p < 0.014) and Snow gain (ρ = 0.63, adj-p < 0.001).
Nylandssjön has only small ephemeral inlets that flow during periods with a very high groundwater table. Such conditions were almost exclusively confined to the snowmelt period, and thus, the transport of terrestrial matter from the catchment to the lake was highly dependent on the intensity of spring run-off during the snowmelt period (Romero-Viana et al., 2008;Petterson et al., 2010). A lower spring snow loss would imply a less intense spring run-off, and thus, a lower transport of mineral material from the catchment to the lake. For Nylandssjön, Rydberg and Martínez-Cortizas (2014) found such a link between the input of coarse grained mineral material to the winter conditions, and such a relationship between PC1 geo and spring snow loss can also be observed in the 40 years record used in this study, i.e., during the mid-1990s when both spring snow loss and PC1 geo scores are low. Similarly, an early ice-off would indicate a longer, and thus, less intense spring run-off period and a lower transport of heavy mineral material to the lake. In our data we see periods with coinciding early ice-off timing and low PC1 geo scores several times in the most recent period (2002, 2007, and 2011). The PC2 geo -scores -where positive scores represent a higher proportion of fine grained detrital material and negative scores a higher proportion of biogenic silicawere not related to spring snow loss, but they are weakly, but non-significantly, negatively correlated to the timing of ice-off (ρ = −0.30, adj-p < 0.085). An early ice-off (low numbers) might imply a prolonged and less intense spring run-off, which will be enough to transport the fine-grained detrital material that give positive PC2 geo -scores but not the heavier mineral particles responsible for higher PC1 geo -scores. The timing of iceoff might also affect the biogenic silica production by diatoms in the lake, and thus, also the PC2 geo -scores (Maier et al., 2019). The PC3 geo -scores were not significantly correlated to any of the tested environmental variables, which supports the interpretation that it is a non-weather related factor. This is consistent with the observation that farming practice was the main driver of PC3 geoscores variation. Altogether, our results indicate that weather conditions do affect the amount and type of terrestrial matter that reached Nylandssjön sediment, however, the relatively weak correlations also imply that there are other factors driving the transport of minerogenic material and nutrients to the lake.
The biomarkers for terrestrial OM inputs were not significantly correlated to any of the tested meteorological variables (Table S2). However, both lignin and phenolic compounds were positively correlated to PC2 geo -scores (Table S2). This shows that concentrations of both lignin and phenolic compounds-which are primarily associated with terrestrial and more degraded organic material  -are favored by a larger input of material from the catchment even if this also implies a higher input of mineral material and reduced total organic matter concentrations in the sediment. In contrast, chlorophyll, protein and steroid concentrations-which were primarily associated with fresh organic material produced within the lake )increased when the mineral material input decrease and the total organic material concentrations in the sediment increase. This pattern is consistent with previous evaluation of C:N-ratios in fresh sediment material from Nylandssjön showing that most of the organic material deposited in the sediment is produced within the lake (Gälman et al., 2006).

Long-Term View of the Aquatic Microbial Community and Its Regulators
The present study brings new insights about the long-term temporal dynamics of the aquatic microbial community (i.e., 40 year) and its regulators at annual resolution. Thus far, most DNA-based monitoring studies have focused on the response to seasonality  and timeframes covering <1 decade (Fuhrman et al., 2015 for a review focused on marine systems).
We identified the inputs of terrestrial material (i.e., PC1 geo scores) as the only regulator of the structure of Nylandssjön's protistan community based on the results from the env.fit analysis ( Table 2). Two biomarkers of algal OM (chlorophyll and proteins) were correlated with changes in cyanobacterial abundance (ρ = 0.49-0.53; adj-p < 0.015; Table S2). Among them, chlorophyll compounds, which are often used as proxy for algal productivity (e.g., Huot et al., 2007), were also found to correlate with fall mean temperatures (ρ = 0.48; adj-p = 0.025, Table S2), which may be explained by the peaks in primary productivity observed during fall (Maier, 2017). PC1 geo and spring snow loss were identified the main parameters reflecting processes occurring in lake catchment that may shape the structure, composition and diversity of microbial communities, the change point analysis capturing the occurrence of their major changes over the 40 years' time series. They were thus selected as predictor's variable for the GAMs analysis, the responses variables selected being the scores of the PCoA axis 1 and 2 (that were reflecting the community structure temporal changes), the diversity metrics of the protistan communities and the richness and relative abundances of the 11 most abundant protistan groups and the Cyanobacteria abundance. The GAMs analysis showed that PC1 geo and spring snow loss explained together 86 and 49% of changes in, respectively, PCoA axis 1 and PCoA axis 2 ( Figure 5; Figure S6). Their contributions to changes in the biological community, were detected mainly during the late-80s/mid-90s (for both predictors) and mid-2000s (for PC1 geo only), concomitantly with an increase of total richness and the occurrence of abundant OTUs within Chlorophyta and Dinophyceae showing that they were the main drivers of community structure changes at these change points, more particularly the PC1 geo while spring snow loss appear to have weak or no effect on shaping the community. Thus, the decrease of terrestrial inputs from the early-80s, and even stronger decline in the mid-90s and mid-2000s, was a key factor shaping Nylandssjön microbial community within the 40-year period studied.
We have considered the occurrence of alternative stable states, as conceptualized in Beisner et al. (2003), in Nylandssjön's microbial community in relation to environmental disturbance. We suggest that the presence of consecutive years with low spring snow loss and thus low inputs of terrestrial material caused a shift in the structure, diversity and composition of the community from late-80s followed by a resilience period to the mid-90s (Figures 2, 3). As defined by Shade et al. (2012), environmental conditions may have acted as a press disturbance on the lake community that reached a break point in late-80s. The short period of change is particularly marked by the decrease of Chlorophyta relative abundance in favor of unclassified taxa from Alveolata and Stramenopiles groups (note that the OTUs from those two last groups were not clearly identified due to incomplete coverage in the reference DNA sequence databases). Interestingly, Apicomplexa becomes highly abundant in this period, but never before or after. Little is known about the ecology of Apicomplexa in lakes except their capability to be fish parasites (e.g., Whipps et al., 2012), and changes in fish populations may also have enhanced such a shift in the trophic network. After the mid-1990s (i.e., end of the resilience period), we may have recorded a new stable state marked by the increase of Chlorophyta relative abundance as well as Dinophyceae, this last taxonomic group being not abundant before the resilience period. The mid-2000s is marked by the increase of photosynthetic groups such as Scenedesmaceae (within Chlorophyta), Dinophyceae (unclassified and may thus also be heterotrophic taxa) and Cyanobacteria. Consistent with this, the varve deposited in 2006 visually appears enriched in algal matter given its strong green color (Figure 1). Altogether, these results highlight the presence of unidirectional changes in the phytoplanktonic community during this time period (as shown by the PCoA plot, Figure 3). While snow phenology appears to be one of the driver of community change-as this change was related to modifications of terrestrial inputs-the lake ice phenology, identified by the presence of several years of early lake ice-off in the 2000s, is a second factor that may have indirectly shaped the temporal dynamics of Nylandssjön's aquatic community. This is consistent with studies that have highlighted the winter/spring transition as a key event for the temporal dynamics of aquatic communities (Lizotte, 2008;Salmi and Salonen, 2016;Grosbois et al., 2017;Hampton et al., 2017;Warner et al., 2018;Maier et al., 2019). The increases of the recurrence of such environmental changes (e.g., early lake ice off, low spring snow loss) and their consequences on Nylandssjön's phytoplanktonic community may be early warnings of the current warming. These factors could result in an increased occurrence and intensity of algal blooms in this lake, a phenomenon already well known to occur in various type of lakes (Paerl and Huisman, 2008;Kosten et al., 2012;Posch et al., 2012).

Effects of Inter-annual Environmental
Variability on Aquatic Microbial Community: The Limit of Annual or Sub-annual Monitoring Data to Isolate the Responses to Longer Term Pressures Due to increasing global warming, a better understanding of the role of climate in past and current changes in lake ecosystems is of key interest. For instance, Capo et al. (2016) showed that protistan communities were driven by the large climatic fluctuations over the last 2000 years, including the Medieval Warm Period, the Little Ice Age and current warming. The longterm effects of an increasing warming trend on phytoplanktonic groups was also shown for several peri-alpines lakes at a centennial timescale (Capo et al., 2017a;Monchamp et al., 2018). However, our findings suggest that the high inter-annual variability of key environmental variables (e.g., weather, snow and ice phenology, terrestrial inputs) may dominate relative to those from global climatic pressures in the Nylandssjön 40-year record. This highlights the limit of drawing firm conclusions concerning drivers of environmental and biological changes at an annual or even sub-annual resolution, and suggests the necessity to study even longer time periods to reveal the effects of climate change.
Anthropogenic pressures, including agriculture and forestry, are identified in the literature as major forcing factors for lake aquatic communities (Jeppesen et al., 2000;Paerl, 2006;Ibelings et al., 2007;Dubois et al., 2018). Although Nylandssjön has not experienced a clear eutrophication, the increase of Chlamydomonaceae in the early-80s may be linked to the transition of local agriculture practices from old-style farming (i.e., natural manure) to industrial farming (with the use of fertilizers) in the late-70s (from 1978 according to information provided by local farmers). This is coherent with the knowledge of the relationships between P concentration and Chlorophyta abundance (Naselli-Flores and Barone, 2000;Soares et al., 2013). Furthermore, using sed-DNA analysis, Capo et al. (2017a) showed significant correlations between high P concentrations (up to 80 µg.L −1 ) and Chlorophyta abundance and more particularly within the Chlamydomonaceae family. Despite the correlations between PC3 geo scores (reflecting changes in P concentrations in sediments) and biological parameters (e.g., diversity metrics of protistan communities), the effects of this local anthropogenic pressures on aquatic communities appear to be overshadowed by the effects of terrestrial inputs in Nylandssjön.

CONCLUSION
The temporal dynamics of Nylandssjön's biological community appear to be driven mainly by the modifications of the inputs of terrestrial matter in the lake that changed between years due to weather variability. Such findings may prevent the use of annual and sub-annual data to predict the changes of biological communities in response to relatively long-term pressures such as climate change or eutrophication. This may be even more important in boreal and arctic lakes in which the snow phenology is a key factor regulating this underlying factor. While calibrations are still required to validate the use of sed-DNA signals as a reflection of the past, as largely discussed in previous papers (Capo et al., 2015(Capo et al., , 2017bDomaizon et al., 2017), this study confirms the potential of sed-DNA to address specific ecological questions. The data here provide evidence about a key time point of changes in boreal lake biological communities: the winter/spring transition, previously shown to affect diatom community (Warner et al., 2018;Maier et al., 2019). The field of microbial ecology would also benefit from DNA-based studies that compare the temporal dynamics of lake communities over contrasting time scales (from days to millennia) to improve our understanding of their adaptive responses to a changing environment.

AUTHOR CONTRIBUTIONS
CB, ID, and EC designed and supervised the study. EC, JR, and JT performed the lab work for the sed-DNA based analysis, inorganic geochemistry, and organic geochemistry, respectively and provided input to the interpretation of results and contributed to the writing of the manuscript. DD provided inputs for bioinformatics. All authors contributed to the revision of the manuscript, read, and approved the submitted version.

FUNDING
Funding for this project was provided by a grant from the Région Rhône-Alpes, the EC2CO INSU program (France) supporting the program REPLAY and the Swedish Research Council (VR, contract nr. 90432301).

ACKNOWLEDGMENTS
We thank Cécile Chardon for her technical help in molecular biology and Pr Gesa Weyhenmeyer for her guidance to model lake ice-off timing from meteorological data. We also thank the Région Rhône-Alpes for their financial support of the Ph.D. thesis of EC through the ARC Environnement as well as the Umeå Plant Science Center for making the Py-GC-MS available to us. We are grateful for the useful comments and suggestions from the two reviewers and the Associate Editor, Gavin L. Simpson.