Sedimentary DNA Reveals Cyanobacterial Community Diversity over 200 Years in Two Perialpine Lakes

Greifensee (47°20′N, 8°40′E) and Lake Zurich (47°13′N, 8°45′E) are two natural perialpine lakes located near the city of Zurich in northeastern Switzerland (Fig. 1). These lakes were chosen on the basis of the availability of long-term data on phytoplankton communities and because recent sediments of both lakes are characterized by the formation of annual varves (Fig. 2). These consist of a pale summer layer and a dark winter layer (27) that allow high temporal resolution dating. Greifensee (surface area, 8.45 km²; maximum depth, 32 m) is monomictic, with one complete mixing event occurring every winter. It is currently classified as eutrophic (average phosphorus of 52 μg/liter in 2015) and has anoxic deep water layers over summer when the lake is strongly thermally stratified (from June to December). Lake Zurich is mesotrophic (average phosphorus of 14.8 μg/liter in 2010) and has a surface area of 65 km² and a maximum depth of 136 m; it is one of the largest perialpine lakes. The lake is divided into two basins by a natural dam, and the focus of this study is on the lower lake, located near the city of Zurich. The lake is considered to be monomictic or dimictic, but the increase of strength of thermal stratification in the last few decades as a consequence of climate warming impedes complete mixing of the water column (28, 29).

The phytoplankton community in Greifensee has been monitored over the last five decades by the Swiss Federal Institute of Aquatic Science and Technology (Eawag). In the present study, we used a data set consisting of species compositions and counts as measured by microscopy from integrated water samples collected over the upper 20 m with a Schröder sampler (30) every month from 1974 to 2010. For Lake Zurich, we used a data set consisting of phytoplankton samples collected by the Zurich drinking water company (Wasserversorgung Zürich [WVZ]) at 14 discrete depths of the water column at monthly sampling intervals from 1976 to 2010. In the two lakes, phytoplankton identification and cell counts were performed using the Utermöhl method (31). The taxonomic affiliations of species were harmonized over the entire data set according to reference 32, and species that were not found at five time points were excluded. This was done to remove taxa that are potentially the result of misidentification or inconsistency in taxonomic identification through time.

Three sediment cores that were 63 mm in diameter and approximately 1 m in length were collected in 2013 using a gravity corer in the deepest part of Greifensee (32 m) and at a 98-m depth in the center of lower Lake Zurich. The cores were sealed and stored in a vertical position in a dark room at 4°C until processing within the following months. One core from Greifensee and two cores from Lake Zurich were opened longitudinally and photographed in a room where no DNA work or PCR amplification had been performed before. The opened cores were visually inspected to identify disturbances that may affect the temporal reconstruction. One half of each core was used as a reference for counting the annual laminations in the sediments. For Greifensee, the reference half of the core was also used for radiometric measurements (see Fig. S1 in the supplemental material) to build an age model. The other half of the Greifensee core and one half of each Lake Zurich core were used for DNA isolation. In the case of Lake Zurich, we used the second core later on to collect a few additional samples between the years 1800 and 1935 in order to complete the time series. All cores were processed at different times to avoid cross contamination. In addition, older sediments were processed independently from recent sediments to minimize contamination.

Varves were formed from the 1930s onward in Greifensee (33) and from the early 1900s onward in Lake Zurich (34) (Fig. 2). The undisturbed varves indicate the absence of benthic bioturbation or sediment mixing (33). In order to select layers of sediments corresponding to specific years over the recent history of the lakes, we built a high-resolution dating profile by varve counting. The Lake Zurich sediment core profiles were aligned with a recent high-resolution age model of a core taken at the same location (A. Gilli, ETH Zürich, personal communication). In Greifensee, we performed radiometric analyses (²¹⁰Pb and ¹³⁷Cs; see Fig. S1 in the supplemental material) (35) to corroborate visual varve counting. Briefly, 20 samples were collected from the half core at 1-cm intervals and were freeze-dried. The sediments were then homogenized, and 3 to 6 g of sediments was used for radiometric measurements on a high-purity germanium (HPGe) well detector (gamma spectrometer) at Eawag facilities. The upper 1.5 cm of each core (corresponding approximately to the years 2009 to 2013) was discarded from the richness comparison because of possible mixing in the surface sediments likely caused during core collection and handling.

On the basis of the age model, we collected sediment strata that corresponded to 1 or 2 years spanning the period for which annual varves were visible (Fig. 2). We collected additional samples from older sediments between the years ∼1840 and 1930 in Greifensee and between ∼1800 and 1900 in Lake Zurich for a total of 20 time points between ∼1840 and 2012 in Greifensee and 23 time points between ∼1800 and 2009 in Lake Zurich. Sediment samples for DNA analysis were taken from the center of the cores while carefully avoiding the sediments in contact with the plastic liner to prevent contamination. Samples were immediately transferred to the DNA clean facility dedicated to eDNA work at Eawag in Dübendorf. Samples were either immediately processed or kept at −20°C until DNA isolation. Strict ancient DNA work protocols were followed in order to prevent contamination with modern DNA (following protocols from reference 36). All instruments used for sedDNA isolation were placed in a laminar flow hood under a UV lamp for a minimum of 20 min prior to use. From each stratum, sedDNA was extracted from approximately 1 g of wet sediments using the PowerSoil DNA isolation kit (Mo Bio Laboratories, Inc., CA, USA) according to the manufacturer's recommendations. The extractions were performed in batches of seven samples, and one negative extraction control (which was treated in the same way as the samples except it did not contain sediments) was added for every seven samples. Each sedDNA extract was quantified using a Qubit (1.0) fluorometer (Thermo Fisher Scientific) following the manufacturer protocol for the double-stranded DNA high-sensitivity assay (dsDNA HS). The genomic sedDNA was inspected for degradation on agarose gels that were dyed with SYBR Safe DNA gel stain. A PCR amplification test was also performed on all sedDNA extracts using published cyanobacterium-specific primers, amplifying an ∼800-nt-long fragment of the 16S rRNA gene (37, 38) to assess the level of DNA preservation in the sediments.

Library preparation for Illumina high-throughput sequencing (HTS) was performed following the optimized protocol of the Genetic Diversity Centre (GDC; ETH Zürich, Switzerland). The protocol uses a two-step PCR approach. A first PCR was performed using previously published cyanobacterium-specific primers CYA359-F, 5′-GGGGAATYTTCCGCAATGGG-3′, and CYA784-R, 5′-ACTACWGGGGTATCTAATCCC-3′ (39). These amplify an approximately 400-nt-long fragment of the V3 and V4 regions of the 16S rRNA gene. The primers were modified for Illumina sequencing by adding overhanging adaptors (see Table S1 in the supplemental material) and inserting up to three random nucleotides between the adaptor and the primer sequence to increase cluster complexity on the flow cell. All modified primer pairs were pooled for the PCR to avoid amplification bias. The specificity of the primers was verified in silico against the Greengenes database (40) to ensure the primer target a wide range of taxa all over the Cyanobacteria phylum, which was the most important characteristic considered.

Sedimentary DNA samples (including negative extraction controls) were amplified in two separate PCRs of 20 μl each, containing 10× Roche FastStart PCR buffer (Roche, Inc., Basel, Switzerland), 25 mM MgCl, 0.2 mM deoxynucleoside triphosphate (dNTP) mix, 0.2 mM each primer (Microsynth, Balgach, Switzerland), and 0.05 U of FastStart Taq polymerase (Roche). The volume of the template DNA used in each reaction varied between 1 and 4 μl depending on the sedDNA concentration. The thermal cycler PCR program included a first denaturation step at 95°C for 4 min, followed by 30 cycles at 95°C for 20 s, annealing at 62°C for 30 s, extension at 72°C for 60 s and a final extension step at 72°C for 5 min. PCR products were purified with 0.8× Agencourt AMPure XP beads (Beckman Coulter, Nyon, Switzerland) and resuspended in 20 μl of 10 mM Tris. The purified product was quantified with a Qubit (1.0) fluorometer.

Depending on the yield, which varied between <0.2 and 7 ng/μl, 2 to 15 μl of purified PCR product was used as the template for a second round of low-cycle PCR to add the Illumina Nextera XT index kit with 1× Kapa HiFi HotStart ReadyMix (Kapa Biosystems, MA, USA). This second PCR was performed with an initial denaturation step at 95°C for 3 min, followed by 8 cycles at 95°C for 30 s, 55°C for 30 s, and 72°C for 30 s and a final extension step at 72°C for 5 min. Samples were purified once more with 0.8× AMPure beads, and the final product was quantified on a 7500 Fast real-time PCR system (ABI) using the KAPA library quantification kit (Illumina).

The indexed samples were pooled at equimolar concentrations into one library. One negative PCR control was also added to the pool. Paired-end 2 × 300 bp sequencing was performed on an Illumina MiSeq (software v.2.5.0.5; Illumina Inc.) at the Genetic Diversity Centre (GDC), ETH Zürich. A total of 6,445,171 raw reads were obtained from 39 samples from Greifensee, and 6,192,583 raw reads were obtained from 43 samples from Lake Zurich.

The sequencing data were quality controlled and processed using an in-house workflow developed at the GDC (see Fig. S2 in the supplemental material). A data quality check of the raw reads was done using FastQC v0.11.2 (41). First, ambiguous and low-quality nucleotides (false-discovery rate [q value], <10) at the end of the reads were removed to improve overlap recognition and error correction. Next, forward and reverse reads were merged into amplicons using USEARCH (v8.0.1623_i86linux64) (42), allowing a minimum overlap of 15 and a minimum merging length of 300 nucleotides. Cutadapt v1. 5 (43) was used to trim the full-length forward and reverse primer sites of the merged reads. Mismatches were not allowed, except for wildcards to compensate for wobble bases. Last, quality filtering and amplicon size selection were performed using Prinseq-lite v0.20.4 (41). The quality filtering, merging, primer trimming, and size selection steps removed about 30% of the data. The negative control had less than 0.14% of reads compared to the average of the other samples and was therefore removed.

The primer-trimmed, quality-filtered, and size-selected amplicons were clustered into operational taxonomic units (OTUs) using the UPARSE workflow (42). The clustering was based on a minimum identity of 97% sequence similarity and an abundance size threshold of 5. De novo and reference based chimera filtering were also applied. All of the steps and parameters applied are detailed in Fig. S2 in the supplemental material and reference 41.

The taxonomic assignment of OTUs was done following UTAX with our own database constructed from the Greengenes database (40) with the addition of a few decoy sequences (see Fig. S2 in the supplemental material) (http://www.drive5.com/usearch/manual/utax_user_train.html). The assignment was done with a confidence threshold of 0.85. OTUs that were assigned to chloroplasts or heterotrophic bacteria were discarded. All except for two photosynthetic cyanobacteria OTUs from Greifensee were assigned to the order (and class) level with high confidence. At the family level, 51 OTUs were assigned, and 27 OTUs were identified to the genus level. Similarly, in Lake Zurich, only one photosynthetic cyanobacteria OTU could not be assigned to the order and class levels. Sixty-two OTUs were assigned a family name, and 31 were assigned to the genus level. None of the OTUs could be assigned a species name with high confidence.

All subsequent data analysis steps were carried out with R software v3.0.1 (44) and the support of various libraries and packages for R and from Bioconductor (http://www.bioconductor.org). First, the OTU count table, including the taxonomic assignment alongside the OTU FASTA sequences and a project-specific metafile, were imported in R using the Bioconductor package phyloseq (45) (see workflow in Fig. S2 in the supplemental material). For comparison with species richness estimated from microscopy data, only OTUs assigned to photosynthetic cyanobacteria were used.

To account for differences in sequencing depth among samples (see Fig. S3 in the supplemental material), the abundance of sequencing reads (cyanobacteria only) for each sample was rarefied to an even sampling depth (10,789 in Greifensee and 10,857 in Lake Zurich) using the rarefy_even_depth function in phyloseq.

The OTU composition of duplicated samples was compared by permutational multivariate analysis of variance using distance matrices (PERMANOVA) using the vegan R package (46). The adonis function was run using the OTU table transformed into a matrix of Jaccard distances based on the incidence of OTUs. The results of the PERMANOVA confirmed our assumption that the sample replicates were not significantly different from one another. For a comparison of richness estimated from water and sediments, linear models were run in R on the rarefied sedDNA samples and on the annual estimates of pelagic species richness in water.

Samples for microscopy were collected monthly between 1975 and 2010 in both lakes. However, in Greifensee, some samples were not collected; therefore, the number of months sampled varied between 8 and 12 per year. To correct for the bias introduced by uneven sampling effort, we estimated the annual species and genus richness by randomly selecting eight time points per year. This random selection was made for 1,001 permutations for each year, and the median was used as our richness estimate. Even though the sampling effort was constant (n = 12) for Lake Zurich over the entire sampling period, the same randomization procedure with permutations to calculate cyanobacterial richness over eight sampling dates per year was performed in order to have comparable estimates of annual richness in the two lakes. To compare annual richness in pelagic and sediment samples, the mean annual richness in water was calculated by averaging the values of three pelagic samples (i.e., year before, year, and year after) around the corresponding year in sedDNA samples. This last step was done to correct for potential inaccuracies in the chronology of the sediment samples.

We used linear models to explain the OTU richness using species richness in water and lake identity as predictors at 12 time points over the years for which both sequencing and microscopy data were available (between 1975 and 2010). We tested for the significance of these explanatory variables using Akaike information criterion corrected for small sample size (AICc). A delta AICc value of ≥2 was considered to be a significant difference.

To verify the accuracy of the taxonomic assignment based on UTAX, we used the OTU FASTA sequences from the two lakes to construct a phylogeny. Sequences were aligned in the software Geneious (47) using the Geneious multiple alignment tool with default settings. Chloroflexus aurantiacus was added as an outgroup, and 13 cyanobacteria reference sequences from GenBank or CyanoBase (http://genome.kazusa.or.jp/cyanobase) were also added. The alignment was used to build a tree based on Bayesian inference (Fig. 3) under 10,000,000 generations using MrBAYES v3.2 (48). The first 25,000,000 generations were excluded at the burn-in step, and the tree's standard deviation of split frequencies was below 0.05.

To confirm the presence of cyanobacteria that can potentially produce microcystin, we used the published primers mcyA-Cd 1R and mcyA-Cd 1F (49). These amplify an ∼300-nt-long amplicon in the mcyA condensation domain of microcystin-producing strains of the genera Planktothrix, Microcystis, and Anabaena. PCR was carried out in a final volume of 50 μl containing 1× PCR Buffer I (Applied Biosystems), 2.4 U AmpliTaq Gold DNA polymerase (Applied Biosystems), 0.5 mM MgCl₂ (in addition to the 1.5 mM contained in buffer), 0.2 mM each dNTP, 0.2 μM of both primers, and 4.8 μg of bovine serum albumin (BSA) (GeneOn, Germany). The PCR program consisted of a first step of polymerase activation at 95°C for 10 min, followed by 35 cycles at 95°C for 15 s, annealing at 56°C for 30 s, and extension at 72°C for 45 s and a final extension step at 72°C for 5 min. The PCR products were visualized on a 2% agarose gel stained with SYBR Safe DNA gel stain. PCR products were purified using the Illustra GFX PCR DNA and gel bands purification kit (GE Healthcare, Little Chalfont, United Kingdom) and directly sequenced (Microsynth, Balgach, Switzerland) to verify the specificity of the PCR product.

The sequence reads obtained in this study were deposited in the European Nucleotide Archive (ENA) under project number PRJEB13044 (http://www.ebi.ac.uk/ena). Unique mcyA sequences from this study have been deposited in GenBank under accession numbers KX437768 and KX437769.

A total of 163 OTUs spanning the phylum Cyanobacteria were recovered from the sediments of Greifensee and Lake Zurich (Fig. 3). In Greifensee, 78 OTUs were obtained from sequencing the 39 sedDNA samples (20 time points from ∼1840 to 2012), while 85 OTUs were obtained in the 41 Lake Zurich sedDNA samples (23 time points from ∼1800 to 2009). Most OTUs were assigned to photosynthetic cyanobacteria, 74% (58 OTUs) in Greifensee and 81% (69 OTUs) in Lake Zurich (Fig. 3). About 25% (20 OTUs) in Greifensee and 19% (16 OTUs) in Lake Zurich were classified as Melainabacteria or ML635J-21, deep-branching groups of nonphotosynthetic cyanobacteria. Interestingly, the two lakes displayed contrasting communities of nonphotosynthetic cyanobacteria, with only 7 (24%) shared OTUs (Fig. 3). Among photosynthetic cyanobacteria, 25 OTUs were unique to Greifensee and 36 OTUs were unique to Lake Zurich (Fig. 3). Thirty-five percent (33 OTUs) of the 95 photosynthetic cyanobacteria OTUs detected in this study were shared (i.e., identified in at least one sample from both lakes).

In order to have a comparable estimate of OTU and species richness, we calculated the annual richness of cyanobacterial orders. The OTU richness within each cyanobacterial order estimated from the sedDNA samples was similar within and between the lakes over the six time points spanning the last 40 years (Fig. 4A and B). The order Chroococcales accounted for most of the annual richness in all samples (40% to 56%). Synechococcales accounted for 13% to 26% of the annual OTU richness, whereas Nostocales, Oscillatoriales, and Pseudanabaenales each accounted for less than 20% of the annual richness estimates in all of the samples.

A total of 42 cyanobacterial species in 26 genera were identified in the pelagic samples from Greifensee between 1974 and 2010 (n = 481), and 37 species belonging to 21 genera were identified in Lake Zurich between 1976 and 2010 (n = 420). Nineteen genera (68%) were found in both lakes. Only 2 genera (7%) were unique to Lake Zurich, whereas 7 genera (25%) were only present in Greifensee. About half of the species (53%) found in Greifensee were unique to that lake, whereas 15 out of 37 species (40.5%) composing the Lake Zurich data set were solely found in that lake. Of all the species listed in the two data sets, 35% were found in both lakes, which is a similar proportion as that for the shared OTUs in the sediments (Fig. 3).

The annual species richness within each cyanobacterial order varied over the years (Fig. 4C and D). In Greifensee, the orders Chroococcales and Synechococcales generally comprised the highest proportion of species annual richness (mean = 32% in both groups), whereas in Lake Zurich, Chroococcales and Pseudanabaenales constituted, on average, 50% of all species each year. Synechococcales, Nostocales, and Oscillatoriales never individually comprised more than 25% of the annual species richness.

Five OTUs predominated (in terms of abundance of sequencing reads) in Greifensee over the six time points between 1974 and 2010 (Fig. 5A), accounting for approximately 96% of the total number of reads. All other OTUs each contributed less than 1% of the reads in that lake. In Lake Zurich, three OTUs were most abundant over the six time points between 1975 and 2010. They accounted for >80% of the reads, and all of the other OTUs each contributed less than 4% of the reads. The first two (Synechococcus spp.) and third (Anabaena sp.) most abundant OTUs in Lake Zurich were also abundant in Greifensee.

The annual species abundances (cells/liter) were calculated from the time series of pelagic samples for the 6 years investigated between 1975 and 2010. In Greifensee, Anathece bachmannii and Aphanocapsa delicatissima (Synechococcales) were the two most abundant species, followed by an Aphanothece sp. (Chroococcales) (Fig. 5C). Planktothrix rubescens (Oscillatoriales) accounted for more than 80% of cell counts at all time points in Lake Zurich, and an Aphanothece sp. (Synechococcales) and Aphanizomenon flos-aquae (Nostocales) were the second and third most abundant species, respectively (Fig. 5D).

The annual cyanobacterial richness varied in each lake, both in the sediments and the pelagic samples. Between 12 and 37 OTUs were found in the sedDNA samples from Greifensee, whereas in the microscopy data sets, the estimated annual richness varied between 6 and 20 species.

In Lake Zurich, the annual richness recovered from the sediments ranged from 12 to 52 OTUs, whereas the estimated annual richness in water varied between 7 and 25 species. There was a strong and significant correlation between the annual species richness in water and the OTU richness in the sediment from the corresponding years (r = 0.81; n = 12). The linear model (Fig. 6) shows a strong and significant positive relationship between pelagic species richness and OTU richness in sediments in both lakes, with the intercepts differing between lakes (adjusted R² = 0.89; n = 12; P < 0.001).

In another linear model, we used genus richness to explain OTU richness in sediments. This was to verify that the relationship between species and OTUs was not caused by biases in fine taxonomic classification. The best linear model explaining OTU richness included both genus richness and lake identity as explanatory variables (adjusted R² = 0.81; n = 12; P < 0.001) (see Fig. S4 in the supplemental matieral).

To better elucidate the differences between the richness of cyanobacteria estimated from sediment and water samples, we used linear models on the annual richness of species grouped by order (Fig. 7). The water-sediment richness relationship was significant in Chroococcales, Oscillatoriales, Pseudanabaenales, and Synechococcales (P < 0.02). Lake identity was a significant factor in explaining OTU richness in all orders except for Pseudanabaenales. The linear model was not significant for the order Nostocales (P = 0.09).

Because potentially microcystin-producing cyanobacteria, such as Anabaena sp., Microcystis sp., and Planktothrix rubescens, have been present in Greifensee and Lake Zurich for several decades, we screened the sequencing data for OTUs assigned to potentially toxic cyanobacterial taxa and verified the presence of mcyA genes in the sedDNA samples by PCR. The conservative taxonomic assignment, based on a high confidence threshold (85%), did not allow us to determine OTUs to a species level, but we identified OTUs corresponding to the aforementioned genera. In Greifensee, two OTUs that are assigned to the genus Microcystis were found sporadically between the 1930s and the 1970s, and either one or both OTUs were detected in all eight samples between 1984 and 2012. In Lake Zurich, a single OTU reference sequence that is assigned to Planktothrix sp. was found in a majority of samples between ∼1800 and 2009.

PCR amplification of mcyA genes confirmed the presence of potentially microcystin-producing cyanobacteria in 4 of 15 sediment layers from Greifensee and in 13 of 19 layers in Lake Zurich (Table 1). Interestingly, in Greifensee, only the most recent sedDNA samples (years 2006, 2009, 2011, and 2012) tested positive for the presence of the mcyA gene. Direct sequencing of the amplicon revealed that the genes found in Greifensee were all related to the same single species of Microcystis. In Lake Zurich, mcyA genes were detected in sediments between the years 1912 and 1962 and in samples between the years 1993 and 2010 (Table 1). A single mcyA sequence related to Planktothrix rubescens was present in 12 of 13 samples from Lake Zurich (the sequencing of the sample dated to 1922 failed).

Our work shows that it is possible to study cyanobacterial communities by sequencing DNA from lake sediment cores. We successfully sequenced amplicons recovered from DNA archived in the sediments of two lakes over the last two centuries, and we were able to validate the data with two independent time series, consisting of 40 years of phytoplankton microscopic identification from the same lakes. We also reconstructed the history of potentially microcystin-producing cyanobacteria over the last century. Our results are consistent with the historical information describing the cyanobacterial community composition in the two lakes.

One of the main limitations in sedimentary DNA studies is the degradation of DNA over time (50). However, the cold and anoxic/hypoxic conditions at the bottom of the two deep and stratified lakes studied here are ideal for DNA preservation (51). We first verified the quality of the DNA preserved over the past 200 years in the sediments of Greifensee and Lake Zurich by amplifying an 800-nt-long fragment of the 16S rRNA gene. This test confirmed the possibility of sequencing a shorter DNA fragment of ∼400 bp. Another important limitation of cyanobacterial investigations using sequencing technologies is the lack of exhaustive and well-curated reference databases, which limits the taxonomic assignment of OTUs. While the existing reference databases (40, 52, 53) are well developed for microbial 16S rRNA analysis, the coverage of cyanobacteria, especially the freshwater taxa, needs to be improved. In this study, however, sequencing a relatively large DNA fragment (400 nt) allowed us to use the OTU reference sequences in a more detailed phylogeny analysis based on Bayesian inferences (Fig. 3).This opens up the possibility of investigating cyanobacterial phylogenetic diversity and community structure.

Because the diversity of whole natural cyanobacterial communities had never been assessed using HTS technologies on sedimentary records, we had no clear expectations regarding the recovery efficiency of our approach. However, based on the time series of pelagic observations, we had some prior knowledge of the cyanobacterial community composition over the past 5 decades in the two lakes. The sequencing of circa 40 sedDNA samples per lake spanning 200 years yielded a similarly high diversity, covering all major clades of cyanobacteria in the two lakes (Fig. 3). Even though the PCR primers were thought to be cyanobacterium-specific, they coamplified chloroplasts and heterotrophic bacteria. This coamplification did not have an impact on community reconstruction because the Illumina sequencing run produced millions of amplicons, which were sufficient for an optimal coverage of the cyanobacterial diversity in most samples (see Fig. S3 in the supplemental material). In general, older samples contained less DNA, which led to a lower number of amplicons sequenced. This problem can be solved by pooling DNA extracts and increasing the template DNA for PCR amplifications or by pooling PCRs. In this study, we were mainly interested in the diversity of recent samples (i.e., between 1975 and 2010) for comparison with pelagic samples; therefore, we did not attempt to optimize the coverage in the older sample.

Interestingly, sequencing revealed the presence of unexpected deep-branching groups of cyanobacteria termed Melainabacteria and ML635J-21 in the sediments of the two lakes (Fig. 3). Recent evidence from whole-genome sequencing confirmed that Melainabacteria constitute a class within the Cyanobacteria phylum because the two groups share common ancestral traits, such as the cell envelope structure and the presence of putative circadian rhythms (54). Taxa within this group have been detected in various environments, including groundwater, drinking water and wastewater treatment plants (54, 55), terrestrial plants, and animal guts (56). To our knowledge, their presence has not been previously reported in lakes. The sequencing of these clades from the sediments of our two lakes shows that the cyanobacterium-specific primers used can target taxa over the whole phylum. From our results, we cannot conclude whether these nonphotosynthetic cyanobacteria live in the water column or colonize the sediments, but the contrasting diversity in the two lakes (Fig. 3) may reflect the adaptation of these taxa to specific local conditions. More eDNA studies using HTS technologies may help to elucidate the ecological roles and the sensitivity of these clades to environmental changes.

The strong and significant relationship observed between the annual cyanobacterial richness estimated at 12 time points from sedDNA and from pelagic samples between the mid-1970s and 2010 in the two lakes (Fig. 6) reinforces the validity of our reconstruction approach. The greater cyanobacterial richness observed in the sediments of the two lakes compared to the microscopic estimates in the pelagic samples (Fig. 6) can be partially explained by differences in the detection limit of the two methods. Other studies have shown that diversity estimation based on morphology generally underestimates the true diversity of cyanobacteria, which emerges from genetic methods (57, 58). Our results suggest that the richness of Chroococcales and Synechococcales was widely underestimated in the microscopy data compared to the genetic data from sedDNA (Fig. 7). This is probably because several taxa within the two groups are unicellular picocyanobacteria (<2 μm in diameter), which are difficult to classify based on morphology (59, 60).

With the sequencing approach, we were able to verify the presence of potentially toxic cyanobacterial taxa that have been observed in the lakes, like Microcystis aeruginosa and Planktothrix rubescens. Regular blooms of M. aeruginosa have been reported in Greifensee over the past 15 years (Eawag, unpublished), and the phytoplankton community of Lake Zurich has been largely dominated by P. rubescens over the last 3 decades (28). Our data confirmed the presence of two OTUs that are assigned to Microcystis species in Greifensee and a single OTU sequence that is related to P. rubescens in Lake Zurich. The detection of the mcyA genes started at the same time that there was an increase in the abundance of one of the OTUs assigned to Microcystis sp. in Greifensee. Although our results show that the number of sequencing reads does not reflect the number of cells counted in pelagic samples (Fig. 5), it is likely that an increase in the relative OTU abundance reflects a change in the pelagic community, in this particular case, the dominance of a microcystin-producing M. aeruginosa genotype.

In Lake Zurich, the detection of an OTU assimilated to Planktothrix rubescens was supported by the amplification of the mcyA gene related to the same taxa. The latter finding is consistent with early reports of the presence of P. rubescens (formerly Oscillatoria rubescens) forming large reddish blooms known as the Burgundy blood phenomenon at the surface of Lake Zurich in 1897 (61, 62) and with a recent study showing that a single genotype of P. rubescens constituted almost 100% of the lake's population over almost 30 years (1980 to 2008) (63).

The relative abundances of sequencing reads did not match the relative annual species densities estimated in the water by microscopy (Fig. 5). Several methodological and biological explanations for this result can be hypothesized. First, PCR and high-throughput sequencing of bacterial 16S rRNA genes introduce biases that can lead to inaccurate population data (64). Also, the number of 16S rRNA gene copies per cell is known to vary among cyanobacterial taxa (65). Traditional microscopy methods are also not free of biases, as plankton identification and cell estimations can vary greatly from one person to another. Rare phytoplankton taxa have been shown to be severely underestimated by traditional sampling methods (66), and it is known that many cyanobacterial taxa are impossible to distinguish on the basis of their morphology only (60). Furthermore, local phenomena in the lake, such as the presence of grazers, buoyant cells, and water currents, affect the sedimentation rate of plankton and may influence the proportion of phytoplankton cells that can be found in the sediments. For the above-mentioned reasons and others, abundance data should be interpreted with extreme caution in sedDNA studies. Nonetheless, the high richness recovered in this study and the strong relationship observed between independent data sets of cyanobacterial community composition (sediments versus water) confirm that the approach for reconstructing past diversity was successful in both Greifensee and Lake Zurich.

This study presents a validated approach to characterize the past composition of cyanobacterial communities archived in lake sediments. Our results demonstrate that amplicon sequencing of a relatively large DNA fragment is useful for investigating the richness and phylogenetic relatedness of cyanobacteria in lakes where the sediments are undisturbed. Our approach, applied to varved sediments, will allow us to explore phylogenetic diversity and community assembly of cyanobacteria over centuries with high temporal resolution. In a forthcoming paper, we investigate the impact of human-induced environmental changes on cyanobacterial phylogenetic diversity and community structure. The ability to recover and sequence important functional genes, like those underpinning the production of secondary metabolites, will assist us in studying the factors that favored toxic cyanobacterial taxa. This approach can, in principle, be extended to other planktonic organisms to help address ecological questions, such as those related to eutrophication, climate change, colonization processes, and invasive species, which are all relevant to the assessment and management of ecosystem processes and services.