Anthropogenic impact on the historical phytoplankton community of Lake Constance reconstructed by multimarker analysis of sediment‐core environmental DNA

During the 20th century, many lakes in the Northern Hemisphere were affected by increasing human population and urbanization along their shorelines and catchment, resulting in aquatic eutrophication. Ecosystem monitoring commenced only after the changes became apparent, precluding any examination of timing and dynamics of initial community change in the past and comparison of pre‐ and postimpact communities. Peri‐Alpine Lake Constance (Germany) underwent a mid‐century period of eutrophication followed by re‐oligotrophication since the 1980s and is now experiencing warm temperatures. We extended the period for which monitoring data of indicator organisms exist by analysing historical environmental DNA (eDNA) from a sediment core dating back some 110 years. Using three metabarcoding markers—for microbial eukaryotes, diatoms and cyanobacteria—we revealed two major breakpoints of community change, in the 1930s and the mid‐1990s. In our core, the latest response was exhibited by diatoms, which are classically used as palaeo‐bioindicators for the trophic state of lakes. Following re‐oligotrophication, overall diversity values reverted to similar ones of the early 20th century, but multivariate analysis indicated that the present community is substantially dissimilar. Community changes of all three groups were strongly correlated to phosphorus concentration changes, whereas significant relationships to temperature were only observed when we did not account for temporal autocorrelation. Our results indicate that each microbial group analysed exhibited a unique response, highlighting the particular strength of multimarker analysis of eDNA, which is not limited to organisms with visible remains and can therefore discover yet unknown responses and abiotic–biotic relationships.


| INTRODUC TI ON
Since the beginning of the 20th century, many freshwater lakes were subjected to anthropogenic perturbations resulting commonly in aquatic eutrophication, due to the faster pace of human population growth and urbanization.The increased use of phosphate-based water softener (pentasodium triphosphate) in laundry detergents and of agricultural mineral fertilizers at that time has been recognized as two of the main contributors of eutrophication of freshwater bodies especially in North America and Europe (Ashley et al., 2011;Jenny et al., 2016;Smil, 2000).In Europe, oligotrophic lakes in the peri-Alpine area experienced eutrophication, which was first noticed in, and monitored since the 1950s for Lake Constance (Figure 1a) (Wessels et al., 1999) (for more details, see Material and Methods).In response, from the 1960s onwards, measures were taken to reduce phosphorus loads that led to improvement of water quality in these aquatic systems and re-oligotrophication ensued, as was the case for Lake Constance since the 1980s (Anneville et al., 2005;Güde et al., 1999;Stich et al., 2018).While the phosphorus level of today indicates a return to levels of the early 1950s (Figure 1a), the resilience, response and reversibility of the lake ecosystems are not sufficiently understood.Even for lakes with a history of monitoring reaching back for decades, there is a lack of knowledge about (a) the lake's biota diversity and composition in the pre-eutrophication state, (b) the exact timing and dynamics of community change in response to rising phosphorus load and (c) the trajectories of the communities during re-oligotrophication when pre-eutrophication communities are considered as a reference.An ecosystem's ability to respond to external environmental impacts relies on its capacity to absorb changes to its physiological state and keep persisting F I G U R E 1 (a) Annual average concentration of phosphorus in Lake Constance water as monitored since 1950 until 2017 (dark red) and diatom-inferred phosphorus concentration (light red) as reconstructed from a Lake Constance sediment core.Box I depicts the rate of change of phosphorus levels as calculated by posterior simulation of continuous intervals; bold black indicates significant increase and bold grey significant decrease.(b) Air temperature anomalies illustrated as the homogenized annual fluctuations of air temperatures in the Upper Lake Constance region (of annual averages expressed in 10-1°C), as recorded by the Historical Instrumental Climatological Surface Time Series Of The Greater Alpine Region (Auer et al., 2007).Box II depicts a significant increase in the rate of rise in air temperatures since the 1950s.(c) Image of the reference sediment core (used for dating) as obtained from Upper Lake Constance showing the laminated layers.The five main trophic phases across the past century are identified according to dating of the individual strata using known flooding events.(d, e) Gam fitted regression curves depicting changes in the effective species numbers in these three data sets as calculated from the logarithmic function of the effective number of species of Shannon-Wiener index (ENS-Shannon) and multiplicative inverse of the Berger-Parker values as reconstructed by selective PCR amplification and barcode sequencing of eDNA archived in the sediment core from Upper Lake Constance, when targeting microbial eukaryotes (dotted dark green), cyanobacteria (dashed green) and diatoms (solid brown line)  (Hodgson et al., 2015).In this study, we define an ecosystem as being "ecologically resilient" for its ability to resist change through self-organization to retain its essential structure and as "ecologically reversible" for its ability to recover from the state of disturbance back to its original state (Allison & Martiny, 2008;Holling, 1973;Quinlan et al., 2016).
While anthropogenic change in trophic states has been a key factor for determining community composition of lakes during the past century (Figure 1a), anthropogenically induced increase in annual average air temperatures (Figure 1b) is thought to be of high relevance for current processes of peri-Alpine lakes (Anneville et al., 2005;Capo, Debroas et al., 2017;Monchamp et al., 2019;Stich & Brinker, 2010).During the last decades, water temperatures increased globally in lakes (O'Reilly et al., 2015) and in Lake Constance (Fink et al., 2014;Straile et al., 2003).A number of recent studies suggested consequences of rising temperatures on the aquatic communities of peri-Alpine lakes (Anneville et al., 2005;Berthon et al., 2014;Capo et al., 2016;Capo, Debroas et al., 2017).
In aquatic ecosystems, bioindicators act as natural "gauges" that are used to study the responses of ecosystems to disturbances.The development and use of several bioindicators prominently started in the 1960s following the recognition of widespread effects of global industrialization and urbanization (Carignan & Villard, 2002;Holt & Miller, 2011).This includes diversity and abundance metrics observed from plants, algae, worms and other invertebrates (Dean, 2008;Fott et al., 1999;Jones et al., 1989;Marbà et al., 2012;Sommer et al., 1993;Wessels et al., 1999).Particularly, the richness, diversity and abundance of planktonic organisms are well-established bioindicators because they quickly respond to changes in nutrient load due to their fast generation times, high growth rates and high dispersal potential (Smith et al., 2016;Sommer, 1986).
They are also commonly used for long-term palaeoecological studies, where organismal remains retrieved from dated sediment cores are used to reconstruct the historical anthropogenic impacts on ecosystems (Battarbee et al., 2001;Catalan et al., 2013;Flower & Battarbee, 1983;Willis et al., 2010;Zillén et al., 2008).Here, diatoms are commonly used as palaeo-bioindicators of the trophic state of lakes, as their siliceous valves are deposited in sediments and can often be well identified (Bennion, 1994;Juggins et al., 2013).Such an analysis on a sediment core from Lake Constance demonstrated that diatom remains could track the historical changes in trophic state, add information about community changes to premonitoring times and suggested a recovery of the diatom community to its pre-eutrophic state (Milan et al., unpublished;Wessels et al., 1999).
To date, DNA metabarcoding (Taberlet et al., 2012) has not been employed to compare the temporal patterns of resilience, response and reversibility of aquatic biota across different taxonomic groups using a multibiomarker approach targeting sediment DNA.
In the present study, we selectively amplified and sequenced the

| Study site
Lake Constance (geographic midpoint: 47°35′N 9°28′E) is a large peri-Alpine lake of glacial origin (surface area, 535 km 2 ).It is a drinking water source to millions in Switzerland and South Germany and a significant economic resource for Austria, Switzerland and Germany (Petri, 2006).
Lake Constance is classified into two major basins, western Lower Lake (maximum depth, 40 m) and Upper Lake Constance (maximum depth, 251 m).The total inflow from the Alps into this region amounts to 12 km 3 /year (Gilfedder et al., 2010).Historically, like many other Alpine lakes (Jenny et al., 2016;Vollenweider, 1970), Lake Constance underwent major eutrophication during the 20th century.Its ecological state was monitored since the 1950s chemically and biologically using longterm (seasonal and annual) plankton measurements (Güde et al., 1999;Sommer, 1986;Sommer et al., 1993).From an assumed steady-state condition in the early 1900s (Grim, 1939;Güde et al., 1999), at which time phosphorus levels were estimated to range < 6 µg/L, the combined effect of population growth, untreated sewage channelling and intensified agriculture during the impending anthropogenic heavy era, caused a shift in its trophic state.Most notably and documented, the 1950s saw an overall increase in the anthropogenic phosphorus loads by a factor of 15-20 fold (Güde et al., 1999).To counteract this, the International Commission for the Protection of Lake Constance (IGKB, 1959) was established in 1959, sewage treatment programmes were implemented around the entire catchment area, and a ban on phosphate-containing detergents was implemented in Austria, Switzerland and Germany.This allowed for the system's successive recovery in terms of phosphorus load, from a peak of 86 µg/L in late 1970s to about 20 µg/L in the 1990s and recent measurements between 6-8 µg/L (Güde et al., 1999;IGKB, 1959;Stich & Brinker, 2010).

| Sediment coring and subsampling
The sediment core studied (BO17/149) was obtained in November 2017 from the Friedrichshafen bay (UTM coordinates, 32T 5,334,876 5,274,382) at a water depth of 187.6 m using a gravity multicorer (IGKB, 2009).Based on previous studies (Kober et al., 1999;Wessels, 1995;Wessels et al., 1995), the site was chosen at a wellcharacterized location in terms of sedimentation rates, visible lamination caused by flooding events (used for dating purposes) and plankton characterization during the eutrophication and the subsequent reoligotrophication (Figure 1a) (Straile, 2015;Wessels et al., 1999).The sediment subsampling process was performed on the same day of field coring during which strict precautions were maintained to avoid contaminations with modern DNA and cross-contaminations from other samples (see Appendix S1 for a detailed description of the procedures).

| DNA extraction and amplification
DNA was extracted from a total of 105 samples using the DNeasy Powersoil kit (Qiagen) following the manufacturer's instructions and using approximately 0.25-0.50g of sediment.DNA concentration and integrity were measured using Qubit (ThermoFisher Scientific).Three sets of primers were used: for microbial eukaryotes (EUK), we used the V7 region of the general eukaryotes 18S rRNA gene, amplifying a 260bp fragment (Gast et al., 2004;Van de Peer et al., 2000); for diatoms (DIA), a variable region of the rbcL gene, amplifying a 310-bp fragment (Vasselon et al., 2017); for cyanobacteria (CYA), a cyanobacteria-specific fragment of the 16S rRNA gene of 400 bp (Monchamp et al., 2016;Nubel et al., 1997) (see Table S1).To build the NGS libraries, a two-step protocol was adapted (see Appendix S1 for a detailed description).The final amplification products were pooled at equimolar concentrations and sent to GATC/Eurofins (Konstanz, Germany; for EUK and CYA) and Fasteris SA (Geneva, Switzerland; for DIA) for sequencing on Illumina MiSeq platforms (2 x 300 bp).
For taxonomic assignment, the filtered-merged reads were aligned using different methods and against multiple databases; (a) the curated 18S /16S rRNA genes database Silva V.132 (Quast et al., 2013) and custom rbcL DIA specific NCBI database using the plug-in feature-classifier classify-consensus-vsearch in QIIME-II (Rognes et al., 2016), (b) NCBI (BLAST, 2008) (GenBank database release 230) using command-line blastn followed by MEGAN (Huson et al., 2016) processing to find the LCA (last common ancestor) and (c) the EMBL-EBI ENA (European Nucleotide Archive) database release 140 using the ecotag package in Obitools (Boyer et al., 2014).
Contrary to findings by Capo et al. (2016), a 97% identity cut-off (with high observed confidence levels in taxonomic classification) was found to be more fitting for these particular data sets, giving accurate and consistent taxonomic classification against the three databases.A total of EUK = 7,721, DIA = 1,272 and CYA = 16,863 molecular operational taxonomic unit (MOTUs) were successfully classified (Table S3).
Further filtering of data included the removal of MOTUs (a) that represented more than 20% of the reads assigned to the negative controls (taken during subsampling, DNA extraction and amplification) (as recommended in Wangensteen et al., 2018) and (b) that are unassigned across the three databases.To minimize the effect of PCR bias, the data were further reduced by removing any MOTUs that were present in only one of the triplicates in each sample (Lange et al., 2015).The triplicates of each sample were then combined to one representative sample by adding the remaining read counts of the nonsporadic MOTUs to increase the reliability of the resulting diversity metrics (Haegeman et al., 2013).The taxonomic assignments from the three databases were harmonized to one representative MOTU in the abundance tables displaying the associated number of times a MOTU was detected.
Unless otherwise mentioned, all successive filtering and statistical analysis steps were performed in R studio version 3.5.2(RStudio, 2015).The nonrarefied abundance tables of the MOTUs were used to test (a) the intrasample homogeneity in terms of consistency in the number of MOTUs and their respective abundances and (b) the standard deviation of the total abundance and MOTU means as calculated by the R package "plyr" (Wickham, 2011).
Similarity comparisons were done using a hierarchical clustering tree calculated from Bray-Curtis dissimilarity matrices using the hclust function from the R package "vegan" (Oksanen et al., 2018).From the normalized final abundance tables, MOTUs that were not targeted using our specific biomarker were removed.Within the EUK data set, we removed sequences classified as Embryophytes, Metazoa and as Bacillariophyta for the following reasons: Embryophytes may not be well identified with this marker and terrestrial vegetation is not the target of our analyses.Aquatic plants can theoretically be captured by these primers, but they were not retrieved at our deep water site.Metazoa were found to be sporadically spread and did not provide any discernible patterns.Bacillariophyta (i.e.diatoms) were analysed at a substantially higher taxonomic resolution using the DIA specific primers and were thus not considered further in this data set.Within the DIA data set, we removed the MOTUs classified as non-diatoms (Chlorophyceae and Eustigmatophyceae).
From the CYA data set, we removed bacteria (with the inclusion of Melainabacteria for illustration purposes only in Figure 2c) and any non-cyanobacteria classified to the Oxyphotobacteria.With regard to the Oxyphotobacteria group, we retained these if the MOTUs were classified as "unclassified" or "uncultured cyanobacteria" in comparison with the reference databases used and through manual checking (BLASTN) (Table S4).We based all further inferences on historical biodiversity analyses on these groups.
A rarefaction of all remaining MOTUs in the three data sets was performed, resulting in normalized final abundance tables of the genetic markers.This was done using the rarefy function from the R package "vegan" (Oksanen et al., 2018).The resulting rarefied tables were used to display rarefaction curves.Following rarefaction, normalization using the mean proportion method was implemented, in which all samples were normalized to total read numbers of EUK = 41,964, DIA = 50,914 and CYA = 49,982 (Table S4).At this cut-off, full saturation of taxonomic richness was achieved.Due to the relatively low sequencing depths of samples corresponding to years 1897, 1905 and 2010/ 1932 and 1897/ 1905 in EUK, DIA and CYA, respectively, they were identified as outliers and hence removed from further analysis.

F I G U R E 2
Stratigraphy plots illustrating the per cent contribution of the individual taxonomic groups detected in the three data sets, of (a) microbial eukaryotes, (b) diatoms and (c) cyanobacteria.Diatoms were grouped a taxonomically described previously (Round et al., 1990).
The solid grey lines illustrate each (a-c) two major breakpoints across the past century, based on the stratigraphy-constrained tree constructed from the differences in the Chord distances between the samples (i.e.years), as shown on the left of each graph.The dotted grey lines illustrate the two equivalent statistical breakpoints in time, based on the change in the community composition (=PC1 values).
To illustrate the annual trend of the change in community composition, the raw PC1 values were added to each group on the right of each graph.Across all three data sets, the trend illustrates a community that has responded to the historic abiotic disturbances, culminating with a shift to a new community structure with the exception of the microbial eukaryote group, for which the community is suggested to have recovered to its 1950s pattern.).This was achieved using the derivative function with 10,000 simulations of finite distance of 1e−07 in "vegan" R package (Simpson, 2019).The samples taken throughout the core were classified according to their trophic phase, that is as oligotrophic (54 cm-39.5 cm), mesotrophic (OE) (38 cm-34.5 cm), eutrophic (32.5 cm-19 cm), mesotrophic (ER) (16 cm-14 cm) and re-oligotrophic (12.5 cm-Surface) (as illustrated in Figure 1a).This was done based on the trophic state index for lakes (TSI) inferred according to the phosphorus loads (Carlson, 1977) and the past characterization of the trophic state of Upper Lake Constance based on deposited sediment diatom frustules (Wessels et al., 1999).

| Alpha diversity
To test the effect of both phosphorus concentrations and air temperature on the general alpha diversity trends, we generated GAMMS from the "mgcv" package in R (Wood et al., 2016).
Subsequently, the MOTUs' Hill numbers, the effective species number Shannon (ENS-Shannon) and the Berger-Parker dominancy indices were calculated using the diversity function in the R package "fossil" (Vavrek, 2011) and the diversityresult function in the R package "biodiversityr" (Kindt & Coe, 2005) to detect the period with the less diverse community.The significance of variation amongst and between the different genomic markers' representative ENS-Shannon diversities, Bartlett and one-way ANOVA (analysis of variance) tests were carried out using the bartlett.testand aov functions.All regression curves were generated with the R package "ggplot2" (Wickham, 2016) to display the temporal contribution to the fitted response variables using generalized additive models (GAMs) (Hastie & Tibshirani, 1986) as the regression method of choice.Henceforth, for the choice of best fit models, the functions gamm and predict.gamfrom the "mgcv" package in R (Wood et al., 2016) were used and GAM parameterization was performed following technical recommendations as in Simpson & Anderson (2009).For all GAMs analyses, only the correlations with p-values < .05were considered significant.Diagnostic plots based on GAMs' diagnostics and ANOVA (Chambers & Hastie, 1992) results, that is; residuals, Akaike information criterion (AIC) values and both complete autocorrelation functions (ACF) and partial autocorrelation functions (pACF), were compared to check the robustness of the models generated.The most abundant groups (≥10% contribution) in the three data sets were plotted stratigraphically along the core, using the percentage of their reads.
This was done using the function stratiplot in the R package "analogue" (Simpson, 2007).To look at the community changes in the overall data set, histograms were generated in "ggplot" R package (Wickham, 2016) using the relative abundances of the taxonomic groups.

| Beta diversity
In order to track changes in the community structure of the three molecular markers studied, ordinations were carried out applying rank orders on the three data sets with the R function metaMDS using Bray-Curtis pairwise distance measures to fit nonmetric MDS (multidimensional Scaling) in the R package "vegan" (Oksanen et al., 2018).The nMDS plots were performed using the functions ordisurf and dendapply.To compare the effect of different data manipulation methods when analysing the changes in the community structure, the data were then first square-root-transformed according to Legendre & Gallagher, 2001 and a principal component analysis (PCA) was performed using the R function prcomp.
Because both nMDS and PCA showed a relatively similar overall change in the community structure, the principal component one (PC1) (showing the highest variation) and two (PC2) were taken as the main representative variation indices of the community structure change.
To investigate whether changes in phosphorus concentrations and air temperature values can be ascertained to drive community changes, PC1 was chosen as a response variable in the GAMMS models generated.For technical purposes, PC2 was also considered to highlight the differences in using different variation percentages.As above, the best fit models were selected based on the AIC, ACF and pACF.In all cases, the values of the phosphorus concentrations were logarithmically transformed for a more linear relationship and the smoothening of the correlation was done using penalized cubic regression splines ("cr") following (Wood, 2017).
Because the time points studied are of variable intervals, the cor-CAR1 correlation method (Singer & Willett, 2003) was used to describe the correlation, with year used as the correction factor.
Additionally, a sliding window analysis of a 13 data points break was constructed to follow the onset of temperature effect on the community change.

| Temporal breakpoints
We analysed the changes in the community composition captured by mainly the PC1 by segmented regression analysis to determine major breakpoints following (Toms & Villard, 2015).This was done using three years as the breakpoint estimators (1920, 1950 and 2000) in the R package "segmented" (Muggeo & Also, 2008) S2).The molecular operational taxonomic units (MOTUs) abundance table obtained from the three data sets were composed of a total of 7,721 (EUK), 1,272 (DIA) and 16,863 (CYA) MOTUs (Table S3).

| RE SULTS
We first investigated the consistency of the DNA signals obtained from the DNA extracted in triplicates from each sediment sample and confirmed that replicates clustered using the intrasample comparison for each sediment layer.This was verified by community structure similarity of intra-versus intersamples (Table S2, Figure S1).Additionally, there was a relatively consistent variation within the number of MOTUs and their relative abundances across intrasample replicates as represented by their standard deviations from intrasample replicates mean (Table S2, Figure S2).Normalized molecular inventories obtained for each sample were found to achieve full saturation of the taxonomic richness as shown by rarefactions curves (Figure S3).

| Temporal changes in the community composition and structure
In as well as the establishment of Aulacoseira, was observed (Figure 2b).
Almost all cyanobacterial groups were not detectable at a high read number in our data set prior to 1930s, at which time some exhibited an abundance peak such as members of Cyanobium and Synechococcus.The temporal trends of the members of those two genera were also marked by an increase in early 1970s and a decrease in the late 1990s (Figure 2c).Within the Oscillatoriophycideae and Nostocales groups, Dolichospermum and Aphanizomenon numbers surged with the onset of the 1970s to decrease again to undetectable levels in the mid-1990s and in the 1980s, respectively.The Aphanizomenon levels rose again only in the early 2000s.In the case of the nonphotosynthetic cyanobacteria groups of an ecological interest, there was a notable increase in numbers of Melainabacteria and Sericytochromatia reads in the 1930-1940s (Figure 2c).
Across the three data sets, a shift was observed in the community using both univariate and multivariate analyses (Figures 1 and   3).The stratigraphically constrained cluster analysis of both chord Figure 1a-c).All segmented regression models were highly significant (p < .001);however, in some cases, regression models with only two breakpoints received similar support as those with three breakpoints.

| Response of the community to environmental changes
Diversity measures calculated across all three data sets showed that there is a community response in terms of the alpha diversity to the change in the trophic status of Lake Constance.(Figure 1).Regarding changes in the abiotic factors "phosphorus loads" and "temperature", the DIA diversity values showed a significant correlation with both, while the EUK diversity values correlated significantly only with phosphorus level changes, and no significant correlation with either could be established with the CYA diversity values (Figure S9a).
Rank-ordered ordination plots illustrate unique communities in each of the five main trophic states of Lake Constance.Twodimensional nMDS and PCA were constructed using Bray-Curtis and Euclidean distances, respectively.Five main clusters were observed with the community structures belonging to each trophic state, while synonymously clustering together with no overlap between the historic oligotrophic and modern re-oligotrophic states (Figure 3, Figure S10a).Across the three markers, different taxonomic groups were shown to occupy the axis, respectively (Figure S10b).Meanwhile, using factor averages of the environmental variables within the nMDS ordinations, the eutrophic samples showed a higher affinity to high phosphorus values (p < .001,r 2 = 0.64), while the re-oligotrophic samples trended towards high temperatures (p < .001,r 2 = 0.65) (Figure 3).
Varying trends were observed when PC1 (37%-48% variation explained) and PC2 (18%-32% variation) were individually tested for strength of correlation with decreasing age (Figure 4).Using PC1, the change in community composition across the DIA and CYA groups is suggested to have altered from its starting point in the early 1900s, while the EUK group's modern-day community change resembles that of the 1950s (Figure 4).Whereas with PC2, values resembling that observed during the historic predisturbance state were reached since the 2000s, (see Figure 4a).
Henceforth, in relation to the abiotic factors studied, both PC1 and PC2 axes were used as representatives of the change in the Best fit GAMMS models representing each of the three data sets illustrating the regression plots of the first axis principal component (PC1) and the second (PC2) across the samples (i.e.years).Thicker lines represent PC1 and thinner lines the PC2.Across all three data sets, the trend illustrates a community that has responded to the historic abiotic disturbances that culminates with a shift to a new community structure with the exception of the microbial eukaryote group, in which the community is suggested to have recovered to its 1950s state.While a unanimous trend can be observed using PC1, with PC2, in the case of CYA and DIA, a somewhat reversible trend is depicted.(b, c) Models representing the rate of change in community compositions using the derivatives of PC1 and PC2 against time.The bold yellow represents the only significant instances of increase and decrease in the rate of change in community composition.In contrast to the temporal trends observed within the GAMMS regression models, the rate of change within the communities as described by both PC1 and PC2 is consistent community composition (Figure S9b).Even though GAMM analysis showed significant correlations between (a) PC1 and PC2 against phosphorus within EUK (p = .016and 0.077, respectively) and (b) PC2 against phosphorus within DIA (p = .032)and CYA (p < .001),no correlations were found against the general increase in temperatures except in the case of PC2 within CYA (p =.039) (Figure S9b).Due to the low R squared values across the models, patterns observed will not be discussed.According to the sliding window GAMM analysis, the contributions of air temperature changes to modifications in the scores of PC1 and PC2 axes appeared to be higher during the time period 1940 to 2002.Meanwhile, phosphorus concentrations remained significantly correlated to the change in community structure (as measured by PC1 and 2) of the three communities throughout the core depths (Figure S1a).

| Lake community changes driven by water pollution
The composition of the species assemblages as recovered with all three markers from the sedimentary archive in Upper Lake Constance showed an overall similar temporal progression throughout the past century.Both our uni-and multivariate analyses imply that the most significant shift in the community occurred already in the 1930s to 1940s, thus before the major eutrophication had been recognized by the classical limnological methods, and that this shift was followed by a second shift, apparent in the community structure (Figure 3), that coincided with the conspicuous mid-century peak of phosphorus.Continuous long-term monitoring data for Upper Lake Constance has been gathered since the 1950s, after signs of eutrophication became apparent (Straile & Geller, 1998), and studies investigating these dynamics thus report changes only from then on (Jochimsen et al., 2013;Sommer et al., 1993).A previous palaeolimnological investigation based on diatom frustules (Wessels et al., 1999) emphasized the importance of eutrophic conditions, but also noted mesotrophic conditions at the end of the 1930s, based mainly on the appearance of Tabellaria fenestrata, that was also regarded as dominating the living community in 1939 (Grim, 1939).A similar early response is known to have occurred also in other European lakes with similar histories (Jenny et al., 2016;Klee & Schmidt, 1987;Lotter, 1998;Smith, 2003), and our multitaxa data suggest that the initial shift in the community of Upper Lake Constance was more severe than the shift to full eutrophic conditions.
Analogous to other peri-Alpine lakes (Jenny et al., 2016), the population history of the region surrounding Lake Constance depicts a spike in the local human population from a historical 0%-15% of the total population change in 1900 to a staggering increase of 30%-45% from 192530%-45% from -195030%-45% from and of 45% from 195030%-45% from -197530%-45% from (Güde et al., 1999)).The early spike in population during 1925-1950 was accompanied by the construction of sewage systems that carried untreated household waste into the lake (Güde et al., 1999) and later the introduction of large industrial facilities in the Friedrichshafen region (Regionalverband Bodensee-Oberschwaben, 2001).Such an increase in urbanization accompanied by the release of untreated sewage is known to be "the" urban point source of nutrients that lead to early hypolimnetic hypoxia in other lakes (Jenny et al., 2016).We therefore hypothesize that the historical urbanization between 1925 and 1950 caused the initial 20th-century disturbance of the community.
The second breakpoint occurred approximately 15-20 years ago, that is between 1996 and 2002, coinciding with the eventual decrease in the measurable phosphorus concentrations in Lake Constance by means of the enforced regulations by the states to improve water quality.Delayed responses to the decline of phosphorus concentrations have been also observed for the phytoplankton biomass in Lake Constance (Jochimsen et al., 2013) suggesting that such a delayed response can also be observed with genomic tools in sediment cores.

| Comparison of early 20th-century communities with their modern counterparts: resilience and reversibility?
Ecological resilience and reversibility denote the ability of an ecosystem to resist change through self-organization to retain its essential structure and to recover from a disturbed state back to its original condition, respectively (Holling, 1973;Quinlan et al., 2016).According to our alpha diversity analysis, the onset of eutrophication (phosphorus > 8.0 µg/L) brought on a collective response of the total community in terms of loss of diversity coupled with an increase of dominancy of certain MOTUs between the 1930-1940s (Figure 1d).An equivalent response to eutrophication, along with an eventual recovery in the early 2000s, has been observed in other lakes of similar history, that is Lakes Baldeggersee (Lotter, 1998), Bourget (Capo, Debroas et al., 2017) and both Greifensee and Gorgova (Monchamp et al., 2017;Monchamp, Spaak, Domaizon et al., 2018) for certain organismal groups.Overall, this points to a lack of resilience, but to a certain degree of reversibility of alpha diversity in Lake Constance.The affinity of some taxonomic groups to low phosphorus inputs is commonly observed in many lakes (Berthon et al., 2014;Cohen et al., 1993;Lefranc et al., 2005;Wessels et al., 1999), where these groups occupy novel trophic niches that became available by the decrease and the disappearance of those taxa that were unable to adapt.
Meanwhile, on the level of community composition, cumulative multivariate analysis depicted a modern community dissimilar to that observed in the early 1900s.After each of the determined breakpoints, 1925-1940and 1996-2002 (Figure 2) (Figure 2), a new distinctive community structure was observed, tracing the major trophic states; oligotrophic (1890-1940), meso-to eutrophic (1940-1990s) and meso-to re-oligotrophic (1990s onwards).In fact, nonmetric multidimensional scaling analysis showed no overlap between the dispersion of the historical oligotrophic and the modern re-oligotrophic communities (Figure 3).Simultaneously, there was a lower and much shorter significant increase in the rate of change in the community composition during the second breakpoint compared to the first as was suggested by Jenny et al., 2016 (Figure 4b).This unique structuring of the communities in the different periods could be due to a number of factors.First, the observed recent state could be transient and a manifestation of the differential ability and speed of each of the taxonomic groups (Figure 2, Figure S4-S6) to react to changing conditions (Donangelo et al., 2012).Second, although phosphorus values drastically decreased again in Lake Constance, the overall chemical composition of the water remains altered due to additional anthropogenic forcing, with, for example, high nitrogen and chloride concentrations (Bericht über den limnologischen Zustand des Bodensees im Jahr 2018, 65.Kommissionstagung).In addition, the current biotic communities may be altered by novel interactions with recent neobiota, some of which are invasive (Hesselschwerdt & Wantzen, 2018).In summary, in the course of the last century, Lake Constance was altered decidedly, while the current communities are putatively shaped both by the history of eutrophication and re-oligotrophication and by additional factors.
In the principal component analyses, we see different temporal patterns, depending on which axes we focussed on.Values of the principal component axis 1 (PC1, explaining 37%-48% variation) are in most cases highly divergent before and after the eutrophication, while PC2 (explaining 18%-32% variation) values do show reversibility (Figure 4a).This discrepancy can be explained through biases introduced by the variables that maximally drive the variation (Chung & Storey, 2015).In our case, changes in the relative abundance of certain MOTUs may be responsible for the disparity in the temporal patterns we observe using the two different PCs (Figure S10b).Within the EUK community, while Fungi and Rhizaria dominate the first axes, Ciliates and the diverse group of Alveolata group along the PC2.Similarly, within the diatoms, while PC1 is driven by MOTUs with affinity to oligtrophic and mesotrophic conditions (i.e. Aulacoseira, Lindavia, Tabellaria and Nitschia), PC2 is dominated by Staurosira.In the CYA data set, the PC1 was dominated specifically by the unidentified Oxyphotobacteria and Synechoccales, which are found to increase in abundance with the onset of eutrophication, followed by a gradual decrease.

| Specific sensitivity of the analysed taxa and their utility and robustness as eDNA bioindicators
Each molecular marker used targeted a different taxonomic group with different ecological characteristics affecting its response rate.Comparing the timing of response to the anthropogenic effects, we detected an earlier response of both the complete microbial eukaryote community and the cyanobacteria than of the traditionally used diatoms.The comparison of DNA data across different taxonomic groups thus demonstrates that a particular strength of eDNA analyses lies in its independence of visual identification-making it possible to analyse groups that have previously not been targeted.
The apparent later response of the diatoms could in part be due to (a) the trophic interactions in play, that is mesozooplankton predation (Chen et al., 2017), or (b) their relatively large biomass that might affect their response rate (Fiore-Donno et al., 2019;Lara & Acosta-Mercado, 2012;Pawlowski et al., 2016;Wickham et al., 2015).However, we cannot speculate further on the exact reason behind this slight shift in diatom response due to a lack of statistical power (Figure 2).The question of whether there was indeed a later response merits further investigation.Ciliates in particular are small, widespread, sensitive to eutrophication (Madoni & Bassanini, 1999) and environmental pollutants (Madoni, 2005) and have a short response time (Pawlowski et al., 2016).In our data set, they showed a strong response to the rise in nutrient loads (Figure 2a).However, it has been suggested that the use of ciliates as bioindicators is hampered by their high rates of evolution, possibly causing an overestimation of unique ecological genotypes (Weber & Pawlowski, 2014).Another group that showed a particularly strong response were the fungi, possibly due to the abundance of zoosporic parasites that are sensitive to pollutants in freshwater systems and are highly diverse (Cudowski et al., 2015;Peay et al., 2016).Both Ciliates and Fungi are thus valuable as sensitive bioindicators that can be targeted in sediment-core DNA (Tikhonov et al., 2015), but their future use would benefit from a better taxonomic and biochemical characterization and understanding of their ecology (Krauss et al., 2011).
Overall, our data also exhibit a substantial robustness, highlighted by the degree of similarity across replicates, such as previous data sets (Lanze et al., 2017).However, community reconstruction using sedimentary eDNA metabarcoding shows some distinct limitations and biases (Domaizon et al., 2017),  et al., 2013;Sjögren et al., 2017).In our case, we could identify many sequences only to the level of genus or family, and this can lead to discrepancies with classical phenotype-based inves- tigations.An example of such a case is the very rare occurrence of sequences identified as the diatom Tabellaria in our data set, compared to its dominance in classical palaeolimnological investigation during the 1950s (Wessels et al., 1999).We furthermore cannot rule out that the overall increase in diversity observed for diatom eDNA, with a particular increase since the beginning of the 1990s, is not caused by DNA degradation.However, regarding our inferences, the breakpoints we identified were not related to a particularly strong change in the number of MOTUs in the core.All samples in the sediment core are subject to similar sources of bias, and the overall temporal trends seem to be reliable and robust.
Essentially, environmental monitoring of ecosystems using either morphological or genetic approaches each have their own limitations and strengths and provide complimentary data.

| Dependence of biotic changes on phosphorus and temperature
In our analyses, only phosphorus concentration in Lake Constance water was found to show an overall significant "relation" to community changes across the past century, while this did not apply to temperature.Multivariate analysis of the data set showed that there is an "affinity" of eutrophic samples towards high phosphorus levels, whereas the re-oligotrophic samples occupied high temperatures.
Previous studies have suggested that the interaction between rising temperatures and re-oligotrophic conditions may be the confounding factors shaping the modern peri-Alpine lake communities (Capo et al., 2016;Jeppesen et al., 2010;Monchamp, Spaak, Domaizon et al., 2018;Monchamp Spaak et al., 2019).The concomitance of both environmental conditions could explain the overall nonresilience of lake community to posteutrophication conditions.
In our investigations, we implemented an autocorrelation with a time covariate (i.e.year) in the multivariate models, and this showed nonsignificant correlations between the "overall" change in the community composition and the rising temperatures except for cyanobacteria.However, when such correction for the time factor-years-was not implemented, a significant correlation could be observed with temperature across all groups.Cyanobacteria's consistent affinity towards climatic warming has also been recently suggested in terms of the last century changes in the phylogenetic structure of cyanobacteria in numerous peri-alpine lakes (Monchamp et al., 2019).The observed lack of significance towards rising air temperatures is supported by measurements of chlorophyll a, in which the nutrient availability carried more significant effect over seasonality in terms of temperatures (Jochimsen et al., 2013;Stich & Brinker, 2010).

| CON CLUS IONS
In conclusion, molecular inventories obtained from a Lake Constance sediment core spanning the last century showed a high degree of protists and fungi are particularly interesting in this regard.This highlights the need for populating the protist and fungal phylogenetic sequence databases (Ishida et al., 2015;Pawlowski et al., 2012Pawlowski et al., , 2016;;Tedersoo et al., 2018;Weisse, 2017) for a better description and investigation of their diversity.For the time period studied here for Lake Constance, our data showed a significant correlation between the overall biotic fauna with changing phosphorus loads.
We cannot, however, substantiate an effect of temperature yet as driver for the general community assembly.The congruence of the three molecular markers in identifying the response, resilience and reversibility of the different lake communities compared to (palaeo-) limnological (counting) data confirm the importance and potential of sediment DNA in environmental monitoring.

ACK N OWLED G EM ENTS
This study was funded by the Deutsche Forschungsgemeinschaft DNA from microbial eukaryotes (phytoplankton and heterotrophic protists) (EUK), diatoms (DIA) and cyanobacteria (CYA), as archived in a finely laminated, well-dated sediment core from Upper Lake Constance spanning the last 110 years.At a high temporal resolution and using available phosphorus concentration and air temperature data, we investigated (a) the timing of group-specific community responses, with special attention to the early 20th century without monitoring data, (b) how resilient, responsive and reversible the freshwater communities have been to 20th century's environmental perturbations, (c) the group-specific sensitivity of the microbial lineages to the anthropogenic modifications and their utility as eDNA bioindicators and (d) the role of changing phosphorus loads and rising air temperatures in lake biodiversity change.
were removed at a 25 base cut-off and a minimum sequencing length of EUK = F:223 and R:179, DIA = F:273 and R:268 and CYA = F:267 and R:255 maintained, (b) discarding of sequencing errors based on a 2M read simulated error model, (c) merging of forward and reverse reads with a maximum accepted error of three bases and (d) putative chimeras removal.The combination of these cleaning procedures allowed for an approximate merging overlap of ± 100 bp.Hereby, any suspected sequencing errors were minimized.A total of 4,090,919 reads for the EUK data set with a mean length of 228.53 bp (std 26.37), 5,085,700 reads for the DIA data set with a mean length of 266.54 bp (std 25.79) and 5,076,107 reads for the CYA data set with a mean length of 393.273 bp (std 13.69) were successfully retained (Table b a e n a O th e r c y a n o b a c te ri u m (O x y p h o to b a c te ri a ) S y n e c h o c o c c u s D o li c h o s p e rm u m C y a n o b iu m P C C -6 3 0 7 G a s tr a n a e r o p h il a le s A p h a n iz o m e n o n S y n e c h o c y s ti s M ic ro c y s ti s c y to c h ro m a ti a O b s c u ri b a c te ra le s O th e r M e la in a b a c te ri a e c h o c c a le s O s c il la t o r io p h y c id e a e N o s t o c a le s M e la in a b a c t e r ia 6000 Changes in phosphorus levels across decades were reconstructed using calculated values during turnover (measurements from biweekly and monthly measured phosphorus loadings in the water column) and the diatom-inferred phosphorus in sediments (Milan, personal communication).Annual monthly mean air temperature anomalies were obtained from the Bregenz climate monitoring station using the Historical Instrumental Climatological Surface Time Series Of The Greater Alpine Region (HISTALP) online database of the mean homogenized temperatures (Auer et al., 2007).The rate of change in both measured phosphorus and air temperatures across the core depths were calculated based on the simultaneous confidence intervals generated from the posterior simulations of the fitted GAMMS (generalized additive models) using the first derivative of a spline (Figure 1a(I) and b(II) . The fitted piecewise linear regression models were further scrutinized based on the p-values, standard deviations and the R 2 values for the integrity of the statistically estimated change points.For an alternative study of the breakpoints in community change, stratigraphically constrained dendrograms were plotted based on Bray-Curtis distance measures (as calculated by metaMDS in the R package "vegan"(Oksanen et al., 2018) and squared chord distances using the mat function in fossil (R package;Vavrek, 2011).
Molecular inventories were obtained from 35 sediment samples spanning the time period from 1897 to 2017.These inventories covered the complex diversity of microbial lineages including microbial eukaryotes (EUK data set, excluding diatoms), diatoms (DIA data set) and cyanobacteria (CYA data set), when using PCR primers targeting the (EUK) V7 region of the general eukaryotes 18S rRNA gene, (DIA) variant region of the rbcL gene and (CYA) cyanobacteria-specific fragment of the 16S rRNA gene, respectively.The DNA concentrations as extracted from the sediment samples (0.25-0.5 g wet weight) ranged from 1.6 ng/µl−6.2ng/ µl (Table Figures S4 and S6).Major shifts in composition in the diatom assemblage included (a) a disappearance of Nitzschia from the 1920s, (b) a subsequent disappearance of Lindavia from 1940 and (c) a decline in a number of other taxa (Diploneis, Amphora, Asterionella and Fragilaria) in the 1930s and 1940s.From the 1930s, we observed a simultaneous rise in Stephanodiscus, Cyclotella and Staurosira, as well as a first appearance of Tabellaria in 1934.The molecular signal of members from the genus Discostella, Cyclostephanos and Sellaphora was not detected prior to 1950s.The latter increased decidedly from the 1990s, during which time a renewed increase of Cyclotella, Amphora and Fragilaria,

(
Figure S1) and Bray-Curtis distances (Figure S7) across the sediment record revealed two major breakpoints in each of the data sets (Figure 2a-c, horizontal solid grey lines).The earlier breakpoint for the EUK and CYA data sets was inferred-surprisingly early-between the years 1932 and 1934 (where a new branch was formed), while for the DIA data set, it was inferred between the years 1940 and 1949.The later breakpoint for the EUK data set was inferred between the years 1996 and 1997, while for both CYA and DIA data sets, it was inferred between 1990 and 1995.Segmented regression analyses of PCA scores with 3 breakpoints suggests for the 1st PCA axis a first breakpoint in the 1920s/1930s for all three groups (CYA: 1925 ± 3 (SE) years, EUK: 1928 ± 3 years and DIA: 1932 ± 2 years).The second breakpoint was inferred in the 1940s for DIA (1943 ± 3 years) and CYA (1949 ± 9 years) and in 1970 ± 5 years for EUK.The final, most recent breakpoints, were inferred earliest for the EUK data set at 1994 ± 4 years, followed by CYA at 2001 ± 5 years and DIA at 2002 ± 6 years (dashed grey lines in The total taxonomic richness and evenness in abundance within both EUK and DIA communities diminished during the years corresponding to the eutrophication of the lake (±1948-1990) and reached the F I G U R E 3 Ordination plots (stress <= 0.06) constructed using Bray-Curtis distances to illustrate the distinction between samples obtained from the different trophic periods.A clustering of samples within their respective trophic state can be observed, while the samples representing the re-oligotrophic state continue to change further away from that of the original oligotrophic composition.Significant affinity was observed with both annually measured phosphorus (p < .05r 2 = 0.639) and the diatom-inferred phosphorus values (not shown) (p < .05r 2 = 0.702) towards eutrophic samples, while temperature (p < .05r 2 = 0.652) towards the re-oligotrophic state c i h p o r t o g i l oe during the peak of eutrophication(1960-  1980s) (Figure1d).Since the 1990s, corresponding to the period of re-oligotrophication, the eDNA suggests a decided increase in diversity, which for the diatoms reached values four to eightfold higher compared to the 1980s.Contrary to this, the CYA community expressed a slight increase in diversity, which dropped with the start of re-oligotrophication (late 1990s).This drop was even more pronounced when no smoothing of the ENS-Shannon number was performed (FigureS8).Conjointly, the dominancy index prescribed by the Berger-Parker values depicted a less homogenous community composition across the three groups during the eutrophic years with a recovery observed with receding phosphorus levels (Figure1e).A nonequal and nonhomogeneous variation was detected between the ENS-Shannon diversity values represented by the individual genomic markers used, as estimated by Bartlett's test (p < .001)and one-way ANOVA (p < .001).A general trend of decreasing diversity was detected and coupled with an increase of phosphorus levels especially > 20 µg/L, which also coincides with the significant rise in temperatures since the 1950s caused in particular by (a) the taxonomic resolution of the chosen molecular marker coupled with the completeness of reference databases, (b) the integration of DNA in the sediment record and its preservation and degradation causing changes in the community of reconstructed sedimentary eDNA compared to living communities and (c) biases caused by experimental manipulation (DNA extraction, PCR, sequencing) causing further differences in the sequence counts compared to the original populations (Kermarrec congruency across the three groups microbial eukaryotes, diatoms and cyanobacteria, but with temporal offsets.Multivariate analyses depict distinct communities in each of the trophic phases, with a first major change in the 1930s-1940s, and an additional, second change until to full eutrophic conditions in the 1970s.After re-oligotrophication, communities shifted again, but these modern communities are different to those of the early 1900s.All analyses revealed two major breakpoints, at which the community composition shifted.The first breakpoint was determined for the microbial eukaryote and cyanobacteria data sets in the 1930s, while the diatom data set showed this shift in the early 1940s.The second breakpoint was determined for all taxa at about 1994-1997.We hypothesize that the first disturbance may reflect the historical increase in human population and urbanization during the 1920s, to which the diatoms, traditionally used as palaeo-bioindicators, showed the slowest response.The use of novel molecular bioindicators, which are not limited to organisms with visible remains such as diatom frustules, is thus a promising tool to extend the time frames of lacustrine biomonitoring, and certain

(
DFG, German Research Foundation) as part of the Research Training Group (GRK2272; project A1 to AM and DS) R3-Responses to biotic and abiotic changes, Resilience and Reversibility of lake ecosystems.We would like to thank all the different laboratory teams involved, of the AGs Meyer, Schleheck, Schink and Epp, for their support and to our graduate school coordinator, Dr. Tina Romer, for her excellent organizational assistance.We also would like to thank both Alexander Fiedler and Ralf Schneider for their technical advice and to Alexander especially for his help in overnight subsampling.We thank Manuela Milan for providing the time series of diatom-inferred total phosphorus concentrations and two anonymous reviewers for their helpful advice.Open access funding enabled and organized by Projekt DEAL.