Lake Pounui is a small (46 ha), shallow (max. depth 9.6 m), lowland coastal lake situated 14 m above sea level, 30 km northeast of Wellington, New Zealand (41°20'34"S, 175°6'48"E). Lake Pounui’s catchment (627 ha) reaches an elevation of 470 m in the foot hills of the Rimutaka Ranges and is covered by unmodified indigenous beech-podocarp forest (96%), with the remainder in pastoral land cover (Perrie and Milne 2012). Although Lake Pounui has a forested catchment, it now regularly experiences cyanobacterial blooms dominated by Microcystis spp. and Dolichospermum spp.

Triplicate ponar grab samples were taken at a central point in the Lake. Using a sterile spatula, the top 1 cm of the undisturbed surface sediment layer of each grab was placed in separate sterile 400 ml containers and transported chilled to the laboratory. Samples were stored chilled for two days and then processed. The samples were homogenised by gentle hand shaking and 53 aliquots (ca. 0.22 g) from each container were transferred into pre-weighed microtubes (1.7 ml; Qiagen) with a 1000 µl mechanical pipette set to 200 µl and using a sterile wide bore pipette tip (step 1 in Fig. 1). The microtubes were then weighed and sediment weights calculated by subtracting the weight of empty microtubes from the filled tube weight (target weight 0.22 g +/- 0.03 g).

Figure 1.

Schematic diagram of the experimental set up and analytical process used in this study, broken down into six steps. 1) Aliquoting of three independent surface sediment grabs, 2) Addition of DNA treatments, 3) Chilling samples in sets of 3 reps × 13 time points (or 14 time points for control samples, which include a time 0 control), 4) The storage phase where triplicates of each treatment are removed from fridge at specific intervals and frozen for later analysis, 5) DNA extraction of all treatments including controls, using a QiaCube (QIAGEN, Germany), 6) DNA amplification of N. spumigina-spiked treatments and a specific investigation of Microcystis abundance via droplet digital PCR (ddPCR) and metabarcoding of the non-spiked control samples for assessing community structure.

DNA and cell lysate from the brackish-water cyanobacterium Nodularia spumigina were used for the experiment as this species is not present in Lake Pounui. Three spiking treatments were used: 100 ng and 200 ng of DNA extract and 10 µl cell lysate (step 2 in Fig. 1).

To obtain extracellular DNA (exDNA) for spiking, DNA was extracted from a non-axenic culture of N. spumigina (CAWBG21; http://cultures.cawthron.org.nz/). Cultured material (1 ml) was siphoned off into three microtubes and centrifuged (10 mins, 3,000 × g). The remaining liquid media was decanted and the cells transferred into six bead tubes of a DNEasy PowerSoil DNA Isolation Kit (QIAGEN, USA). DNA extraction was automated using a QIAcube (QIAGEN, USA), with DNA eluted in 100 µl volumes. DNA extracts were pooled and quantified with a NanoPhotometer NP80 (Implen GmbH, Munich, Germany) to calculate the required volume for spiking to achieve concentrations of 100 ng/µl (low concentration spike) and 200 ng/µl (high concentration spike). To obtain cell lysate for spiking, 10 ml of N. spumigina (CAWBG21) culture with liquid media was sonicated (30 secs) and filtered (3 µm; filter type/brand/brand) to remove excess non-DNA cellular contents.

For each treatment, 13 sediment aliquots were spiked in triplicate (step 3 in Fig. 1). The DNA or cell lysate was mixed with the sediment by pipetting the mixtures up and down three times. The samples were stored in the dark in a polystyrene container at 4 °C (n = 39 samples per treatment). Triplicate non-spiked control samples for each time point were chilled alongside the treatments, including a set of triplicate samples for time zero that were immediately frozen at -20 °C (n = 49 control samples). Triplicate samples were removed from chilled storage and frozen (-20 °C) at 1, 2, 3, 4, 5, 6, 8, 10, 12, 14, 16, 20 and 24 weeks (step 4 in Fig. 1).

Each step of the molecular analysis, including DNA extraction, PCR set-up, template addition, PCR amplification, amplicon purification and quantification, was conducted at the Cawthron Institute (Nelson, New Zealand) in separate sterile rooms dedicated to that step, with sequential workflow to prevent cross-contamination. Rooms dedicated to DNA extraction, PCR set-up and template addition were sterilised with ultra-violet light for a minimum of 15 mins before and after each use. The PCR set-up and template addition were undertaken in laminar flow cabinets with HEPA filtration. Aerosol barrier tips (Eppendorf, Germany) were used throughout.

The sediment from each time point (n = 53) was transferred from microtubes into bead tubes of a DNEasy PowerSoil DNA Isolation Kit (QIAGEN, USA) using a pipette and the bead tube liquid to assist collection of sediment from the microtubes. DNA was then extracted using a QIAcube (QIAGEN, USA), with DNA eluted in 100 µl volumes (step 5 in Fig. 1). Extractions were performed in batches of ten to twelve samples, including a negative control (i.e. extraction blank containing all the reagents, but no sediment) every two batches. DNA concentration and quality were measured using a NanoPhotometer NP80 (Implen GmbH, Munich, Germany).

The genus Nodularia is the only known producer of the hepatotoxin nodularin, which is produced via the nodularin synthetase­­­ enzyme complex (Koskenniemi et al. 2007). To quantify spiked extracellular DNA and cell lysate, absolute concentrations of the nodularin synthetase gene subunit F (ndaF) were measured in all samples using a BioRad QX200 Droplet Digital PCR (ddPCR) system (step 6 in Fig. 1). Primers used were the ndaF8452 forward primer 5′-GTG ATT GAA TTT CTT GGT CG-3′ and the ndaF8640 reverse primer 5′-GGA AAT TTC TAT GTC TGA CTC AG-3′ (Koskenniemi et al. 2007). Each ddPCR reaction included 450 nM of forward and reverse primer, 10 µl of 2 × BioRad ddPCR EvaGreen Supermix, 3 µl DNA and sterile water for a total reaction volume of 22 µl. The BioRad QX200 droplet generator partitioned each reaction mixture into nanodroplets by combining 20 µl of the reaction mixture with 70 µl of BioRad droplet oil. After processing, this resulted in a total nanodroplet volume of 40 µl, which was transferred to a PCR plate for amplification using the following cycling protocol; hold at 95 °C for 5 min, 40 cycles of 95 °C for 30 s and 60 °C 1 min, a cooling step of 4 °C for 5 min and a final enzyme deactivation step at 98 °C for 5 min. The plate was then analysed on the QX200 Droplet Reader. For each ddPCR plate run, one negative control (RNA/DNA-free water; Life Technologies) per column of samples and one positive control (genomic DNA extracted from a sample known to contain N. spumigina) per plate were included. When samples were too concentrated (lacking clear negative droplets), these were diluted (1:10) with RNA/DNA-free water (UltraPure^TM, Life Technologies, CA, USA) and re-analysed. The results where then converted to copies per gram using the following formula: number of copies per µl × 22 µl (the initial volume of the PCR reaction) × 33.333 (proportion of total eluted DNA from the extraction step; 100 µl / 3 µl) × 10 (if 1:10 dilution was used in PCR reaction) × 0.22 (sediment wet weight in grams).

To characterise the bacterial prokaryote and eukaryote composition over the 24 weeks of storage, the V3-V4 regions of the prokaryote 16S rRNA gene and the V4 region of the eukaryote 18S rRNA gene were amplified by PCR, using the prokaryote-specific primers 341F: 5’-CCT ACG GGN GGC WGC AG-3’ and 805R: 5’-GAC TAC HVG GGT ATC TAA TCC-3 (Herlemann et al. 2011; Klindworth et al. 2013) and eukaryote primers Uni18SF: 5’-AGG GCA AKY CTG GTG CCA GC-3’ and Uni18SR: 5’-GRC GGT ATC TRA TCG YCT T-3’ (Zhan et al. 2013), each including adapter tails for dual-index sequencing (Kozich, Westcott, Baxter, Highlander and Schloss 2013). PCR reactions were performed in 25 µl volumes for each sample with the reaction mixture containing 12 µl of MyFi Mix, 5 µl of RNA/DNA free water (UltraPure, ThermoFisher Scientific), 1 µl each of forward and reverse primer (10 µM) and 6 µl of template DNA. The reaction mixture was held at 95 °C for 1 min followed by 35 cycles of 95 °C for 15 s, 50 °C for 15 s, 72 °C for 15 s, with a final extension step at 72 °C for 7 min. Negative PCR controls and DNA extraction controls were run alongside the samples for each PCR analysis. PCR products from the triplicates were pooled and 25 µl of the pooled product was cleaned and normalised using SequalPrep Normalization plates (ThermoFisher Scientific, USA). This resulted in a cleaned DNA concentration of ~ 1 ng µl^-1. Samples were sent to Auckland Genomics Facility for sequencing on an Illumina Miseq^TM platform (step 6 in Fig. 1). The samples were prepared for sequencing following the Illumina 16S metagenomics library prep manual with the exception that, after the indexing PCR, 5 ml of each sample (including water samples acting as sequencing blank) was pooled and a single clean-up was undertaken on the pool instead of samples being individually cleaned. The library pool quality control was undertaken on a bioanalyser before the library was diluted to 4 nM, denatured and diluted to a final loading concentration of 7 pM with a 15% PhiX spike. Raw sequence reads are deposited in the NCBI short read archive under the accession number: PRJNA780583.

Bioinformatic pipelines for both the 16S and 18S rRNA genes were identical unless otherwise stated. Primers were removed from the raw reads using cutadapt with one mismatch allowed (Martin 2011) and subsequently processed using the DADA2 package (Callahan et al. 2016) within R (R Core Team 2018). Reads were truncated (230 and 228 base pairs (bp) for forward and reverse reads, respectively) and filtered with a maximum number of “expected errors” (maxEE) threshold of two forward reads and four reverse reads; reads not matching this criterion were discarded from further analysis. A parametric error matrix was constructed, based on the first 10⁸ bp of the sequences. Sequence variants for the forward and reverse reads were inferred, based on this error matrix following sequence dereplication. Singletons were discarded and remaining paired-end reads were merged with a maximum mismatch of 1 bp and a required minimum overlap of 10 bp. The chimeric sequences were removed using the removeBimeraDenovo command within the DADA2 package.

The resulting Amplicon Sequence Variants (ASVs) were used for taxonomic classification against the SILVA 138 (Pruesse et al. 2007) and PR2 (Guillou et al. 2013) databases for the 16S and 18S rRNA gene datasets, respectively. We chose to use the PR2 database for the classification of eukaryotes because it provides better classification of protists, which are likely to be the most important and abundant organisms in these samples, compared to other available databases, for example, SILVA. The sequences were classified, based on the rdp classifier (Wang et al. 2007) with a bootstrap of 50 to allow classifications at higher taxonomic levels. The results were combined into a phyloseq object (McMurdie and Holmes 2013) and sequences assigned as eukaryotes, chloroplasts and mitochondria were removed from the 16S rRNA dataset, while bacterial sequences were removed from the 18S rRNA dataset. Read numbers for ASVs found in negative control samples were removed from the dataset via subtraction. For comparisons between samples, subsampling to an even depth was undertaken: 13,000 and 8,000 reads for bacterial and eukaryotic datasets, respectively.

To investigate abundances of the cyanobacterial genus Microcystis in the non-spiked samples, a ddPCR assay, targeting the 16S rRNA gene, was used. Absolute concentrations of the 16S rRNA gene were measured in all samples using a BioRad QX200 ddPCR system, Cyanobacteria-specific primers (MICR-184 forward primer, 5′- GCC GCR AGG TGA AAM CTA A-3′; MICR-431 reverse primer, 5′- AAT CCA AAR ACC TTC CTC CC-3′; Neilan et al. 1997) and a Microcystis-specific probe (MICR-228 probe FAM-AAG AGC TTG CGT CTG ATT AGC TAG T-BHQ-1; Rinta-Kanto et al. 2005). Each ddPCR reaction included 450 nM (1 µl) of each primer and probe, 10 µl of 2 × BioRad ddPCR Probe Supermix, 2 µl DNA and sterile water for a total reaction volume of 22 µl. PCR set-up and droplet reading settings and equations used to determine concentrations were as described above for the nodularin assay. The PCR conditions were as follows: hold at 95 °C for 10 min, 40 cycles of 94 °C for 30 s and 60 °C 1 min and a final enzyme deactivation step at 98 °C for 10 min. All samples were diluted (1:100) with RNA/DNA-free water (UltraPureTM, Life Technologies, CA, USA). The plate was then analysed on the QX200 Droplet Reader. For each ddPCR plate run, one negative control (RNA/DNA-free water; Life Technologies) per column of samples and one positive control (genomic DNA extracted from a sample known to contain Microcystis) per plate were included. The results where then converted to copies per gram using the following formula: number of copies per µl × 22 µl (the initial volume of the PCR reaction) × 50 (proportion of total eluted DNA from the extraction step; 100 µl / 2 µl) × 100 (dilution used in PCR reaction) × 0.22 (sediment wet weight in grams).

Statistical analysis of the data was undertaken within R (R Core Team 2018). Linear regression was used to assess the relationship between concentrations of the ndaF gene (copies per gram as determined with ddPCR) and storage time (weeks), for each treatment (spiking with 100 ng, 200 ng and cell lysate). A separate linear regression was used to assess the relationship between concentrations of Microcystis-specific 16S rRNA gene (copies per gram as determined with ddPCR) and storage time (weeks).

Community diversity and composition analyses were conducted for both bacteria (16S rRNA) and eukaryote (18S rRNA) metabarcoding datasets from the non-spiked samples. Observed ASV richness, Shannon and Simpson Indices were considered and linear regressions were used to assess the relationship between each diversity measure and storage time (weeks). Bray-Curtis dissimilarities of individual samples were compared using Principal Components Analysis (PCoA), plotted using the plot_ordination function within the phyloseq package in R (McMurdie and Holmes 2013). Compositions of bacterial phyla, eukaryotic phyla, cercozoan classes and cyanobacterial genera for each week of storage were visualised using the plot_bar function within phyloseq. Relative abundances were calculated for each week of storage by summing the three replicates. Statistical differences in taxonomic composition grouped by storage time (weeks) were tested via permutational multivariate analysis of the Bray-Curtis dissimilarities, using the adonis function within the vegan package (Oksanen et al. 2019).

Copy numbers of the ndaF gene did not increase or decline over the 24 weeks of storage for the 100 ng and 200 ng spiked treatments (R² = 0.067, p = 0.11 and R² = 0.058, p = 0.13, Fig. 2a, b). However, average copy numbers declined in the cell lysate spiked treatment from 315,777 (SD = 90,737, SE = 52,387) copies per gram to 151,360 (SD = 42,334, SE = 24,442) at week 12 and 40,325 (SD = 37,542, SE = 21,675) by week 24, equating to 52% and 87% reductions from initial yields, respectively. (R² = 0.42, p < 0.001, Fig. 2c).

Figure 2.

Linear regression between concentration (copies per g of wet weight) and storage time for the ndaF gene. Surface sediment samples were spiked with 100 ng (a) or 200 ng (b) of DNA or cell lysate (c). Samples were stored chilled (4 °C) in the dark and sampled at weekly to bi-weekly intervals and assessed using droplet digital PCR over a 24-week period.

A total of 17,548 ASVs were reduced to 13,096 ASVs (ranging from 1,987–2,569 ASVs across 37 samples) after quality control, taxonomic filtering and rarefaction. Unassigned sequences at phylum level were retained and comprised 4.1% of remaining ASVs. Over the 24-week period, mean observed bacterial richness decreased from 2,364 (SD = 73.1, SE = 42.2, n = 3) to 2,081 (SD = 62, SE = 35.8, n = 3). There was a weak negative relationship between observed richness and sampling time (R² = 0.23, p = 0.0025, Suppl. material 1: Fig. S1). In contrast, Shannon and Simpson Indices did not change significantly over the study period (R² = 0.041, p = 0.23, R² = 0.007, p = 0.61, respectively, Suppl. material 1: Fig. S1).

Bacterial phyla comprising > 10% of retained sequences (summed across replicates) were Proteobacteria, Verrucomicrobiota, Bacteroidota (Fig. 3). The PCoA analysis showed a clear shift in the bacterial community during the experiment and replicates clustered closely together (Fig. 4a). Multivariate analysis identified significant differences in the bacterial composition across 24 weeks (F = 1.366, DF = 13, 23, R² = 0.435, p = 0.001; Suppl. material 2: Table S1). Changes in phyla relative abundances were visible with a noticeable decrease in cyanobacteria, which dropped from 9.6% at the beginning of the experiment to 0.4% at week 8 and ≤ 0.3% from week 10 to 24 (Fig. 3). Cyanobacteria taxa predominantly comprised the genus Microcystis (57–99% of cyanobacterial sequence reads between weeks 0 and 14 and 24–57% between weeks 16 and 24) and read numbers decreased from 3,738 to 115 between weeks 0 and 8 (Suppl. material 1: Fig. S2). Proteobacteria noticeably increased from 15% at week 0 to 33% at week 24 (Fig. 3). Within this phylum, Gammaproteobacteria were the dominant class (≥ 84%) throughout the 24 weeks and, together with Alphaproteobacteria, these comprised ≥ 99% of sequence reads. Of the Gammaproteobacteria, Burkholderiales (47–73%) was consistently dominant throughout the experiment and Xanthomonadales increased from week 12 to dominate the Proteobacteria sequence counts in week 24 (48%; Suppl. material 1: Fig. S3).

Figure 3.

Bacterial phyla (relative abundance of 16S rRNA gene sequences, summed across replicates, n = 3) of non-spiked samples over 24 weeks of chilled (4 °C, dark) storage determined using metabarcoding. Sequences were rarefied to 13,000 reads. The 15 most abundant phyla are shown with the remainder grouped as ‘Other’.

Figure 4.

Principal Coordinates Analysis (PCoA) of bacterial 16S rRNA gene (a) and eukaryotic 18S rRNA gene (b) metabarcoding data using Bray-Curtis dissimilarity coloured by storage time (weeks).

Microcystis 16S rRNA gene copies declined in a non-linear fashion over time (R² = 0.78, p < 0.001; Fig. 5). From week 0, average gene copies (2.2 million, SD = 725,736, SE = 419,004, n = 3) declined by 60% to 911,150 after 3 weeks and 99% to 12,736 after 24 weeks.

Figure 5.

Concentrations (copies per g of wet weight) of the Microcystis 16S rRNA gene in non-spiked samples over 24 weeks of chilled (4 °C, dark) storage (Loess smoothed trend-line with 95% confidence intervals).

Linear regression of eukaryote taxa in non-spiked samples did not show any relationship between observed richness and time, but there was a very weak positive relationship when using Shannon and Simpson Indices (R² = 0.2, p = 0.0053 and R² = 0.22, p = 0.0033, respectively, Suppl. material 1: Fig. S1). Mean Shannon Index increased by 20% from 2.75 at week 0 to 3.29 at week 24 and mean Simpson Index increased 13% from 0.77 at week 0 to 0.87 at week 24.

A total of 3,462 ASVs decreased to 2,736 ASVs (ranging from 108–369 ASVs across 37 samples) following quality control, taxonomic filtering and rarefaction. Unassigned sequences at phylum level were retained and comprised 13.4% of remaining ASVs. Eukaryote phyla comprising > 20% of retained sequences (summed across replicates) were Annelida, Arthropoda, Cercozoa and Dinoflagellata (Fig. 6). The PCoA analysis showed convergence of the eukaryotic community in samples stored for longer periods (Fig. 4b). Permutational multivariate analysis showed significant differences (F = 1.314, DF = 13, 23, R² = 0.426, p = 0.012; Suppl. material 2: Table S1) in eukaryotic composition across 24 weeks of storage. Changes in the relative abundances of eukaryotic phyla were visible with a noticeable increase in cercozoan sequences, which increased from 2.95% at the onset of the experiment to 53.4% at week 10 and 30.5% by week 24 (Fig. 6). Relative abundance of phylum Annelida and Arthropoda fluctuated across the 24 weeks of storage, ranging from 0.03–59% for the first 5 weeks and reducing to 0–43% for the remainder of the experiment. Within Cercozoa, Sarcomonadea and Thecofilosea were the dominant classes throughout the experiment (Fig. 7). Sarcomonadea comprised > 50% of reads for all weeks, except weeks 5, 6, 8 and 14. Thecofilosea increased in relative abundance from 9% at week 0 to 54% at week 6 and remained above 34% for the rest of the experiment (excluding week 20). Haplotaxida were the dominant order of Annelida (> 99% across all time points) and Calanoida was the dominant order of Arthropoda (> 99%) for nine of the thirteen time points (Suppl. material 1: Fig. S4).

Figure 6.

Eukaryote phyla (relative abundance of 18S rRNA gene sequences, summed across replicates, n = 3) of non-spiked samples over 24 weeks of chilled (4 °C, dark) storage determined using metabarcoding. Sequences were rarefied to 8,000 reads. The 15 most abundant phyla are shown with the remainder grouped as ‘Other’.

Figure 7.

Cercozoa composition and number of reads as determined by metabarcoding of the eukaryotic 18S rRNA of non-spiked samples over 24 weeks of chilled (4 °C, dark) storage.

We found differential effects of chilled storage on DNA yields from sediments artificially spiked with purified DNA and DNA from freshly lysed cells. Consistent with our first hypothesis, we observed no change to DNA yields for purified DNA when spiked at either 100 or 200 ng per sample during the 24 weeks of storage. Similar results have been found in previous studies with little change observed in recovered environmental DNA (eDNA) from soil samples stored for 14 days (Lauber et al. 2010). Our results suggest the purified extracellular DNA quickly became sediment-bound, likely reducing DNA degradation. This is in agreement with experimental studies that found: i) when DNA was added to non-sterile soils, 80% of it bound to soil within 1 hr (Blum et al. 1997) and ii) maximum DNA absorption on to quartz sand was reached after 2 h (Lorenz and Wackernagel 1987). Once DNA molecules are bound to sediment, they can be difficult to cleave off (Romanowski et al. 1991).

Consistent with our second hypothesis, yields of DNA from freshly lysed cells declined steadily across the 24 weeks, with concentrations reducing 87% over the experiment. Compared to the purified DNA-spiked samples, higher starting volumes of DNA in the lysate-spiked samples may have led to an overabundance of DNA versus sediment binding sites. Additionally, cell lysate also contains cellular nucleases that may contribute to DNA degradation before some molecules become sediment bound. Bacteria surviving in chilled sediments can also produce or take advantage of nucleases in the sediment matrix, allowing them to feed on depurinated DNA and subsequent degradation products (Wasmund et al. 2021). Proteobacteria, which increased in abundance during the experiment (see discussion below) are known to have genes that enable consumption of DNA (Finkle and Kolter 2001) and this could be another possible mechanism for reduction in DNA of lysate-spiked samples. These cell lysate samples represented a realistic scenario of a senescing algal bloom, where cells gradually die and settle on the sediment surface. Our results highlight the need for rapid sampling, particularly for eutrophic lakes where blooms are a regular occurrence.

Biological communities, detected in sediment samples, can be a combination of living cells, dormant cells and extracellular DNA (Taberlet et al. 2012). Some bacterial and eukaryotic microorganisms can survive in low energy, anoxic conditions and proliferate by feeding on each other or nutrient compounds in the sediment matrix (Koizumi et al. 2003; Compte-Port et al. 2020; Pearman et al. 2021a). Based on short-term storage studies, we hypothesised that maintaining samples in chilled and dark storage conditions would result in limited shifts in biological communities because the environment is similar to that from which they were collected. The number of bacterial taxa (observed ASV richness) showed a weakly associated decrease with time, but there was no significant change in either diversity index (Shannon and Simpson). This aligns with Lauber et al. (2010), who found similar phylogenetic structure and diversity of soil bacterial communities across storage temperatures.

However, there were several exceptions to this pattern with the most notable being the cyanobacterial genus Microcystis; a bloom-forming microorganism able to produce potent hepatotoxins (Harke et al. 2016). The relative abundance of this genus halved by week 3. To confirm that the reduction of 16S rRNA sequence abundances in the metabarcoding data was not due to sequencing bias, which can occur when community composition changes over time (Lee et al. 2012; Pochon et al. 2013), a Microcystis-specific ddPCR assay was applied and an equivalent decrease in copy numbers was detected. Our results also showed a shift in eukaryotic community composition around 5–6 weeks of storage. This was mainly driven by an increase in Cercozoa; an ubiquitous phylum of single-celled eukaryotes that are common in sediment and functionally and morphologically diverse, including both autotrophs and free-living predators (Bass and Cavalier-Smith 2004; Dumack et al. 2020). The observed change in composition indicates the persistence of live organisms in sediment samples at 4 °C and suggests that they may feed on other taxa and proliferate over time. Although, to the best of our knowledge, Cercozoa have not specifically been shown to feed on Microcystis, some constituents can be bacterivorous or omnivorous (Flues et al. 2017; Dumack et al. 2020). Other unicellular eukaryotes have been shown to feed on Microcystis blooms in both field and laboratory experiments (Zhang et al. 2008; Van Wichelen et al. 2010). The detected increase in Cercozoa suggests a biotic mechanism for the decrease in Microcystis sequences by preferential grazing. Unlike some Cyanobacteria, Microcystis does not produce akinetes, which can be larger than vegetative cells and have thickened cell walls and this may make Microcystis cells more vulnerable to grazing. For studies investigating the history of Cyanobacteria proliferations (for example, Pal et al. 2015 and Monchamp et al. 2016), their abundance could be under-represented if samples are chill-stored for more than two weeks.

Gammaproteobacteria were the dominant class of proteobacterial phyla throughout the experiment, with the order Xanthomonadales notably increasing in abundance after week 12. Similar patterns of increased Gammaproteobacteria abundance were found in experimental liquid microcosms (Flues et al. 2017), where Gammaproteobacteria were positively affected by Cercomonads grazing on the bacterial community and the authors attributed this to the distinct changes in bacterial community composition. We suggest a similar positive effect on Xanthomonadales occurred in our samples, as members of this order are reported to possess anti-predation cell structures that confer a competitive advantage over other bacteria (Bayer-Santos et al. 2019).

Amongst the eukaryote data, there were several other notable patterns including variable numbers of Annelida sequences. As larger organisms, they are likely to have a patchier distribution and it is probable they were living in the sediment samples at the time of collection. Given that there were no discernible trends, we speculate that uneven numbers of them were contained in the sediment during the subsampling process. Calanoida, commonly found as zooplankton, were the predominant order of arthropod detected in most samples and their abundance was also uneven. It is unlikely these arthropods were alive in the sediment samples; therefore, most of the DNA we detected was likely extracellular. The patchy nature of these organisms highlights the need for larger sediment samples when working with larger organisms for precise biodiversity characterisation.

Considering the clear changes in bacterial and eukaryotic community compositions observed in this study, we recommend subsampling and freezing sediment core samples, particularly those near the surface where microbes may be more active, within 2–4 weeks of collection and chilled storage. Immediate freezing of subsamples within 1–2 weeks is recommended for studies where low abundance or rare taxa are important and for assessments of microorganisms (bacteria, cyanobacteria and phytoplankton). Time-frames for broad-scale community structure assessments may still be possible up to 5 weeks after sampling with chilled storage.

Our analysis of surface sediment is relevant for deeper core deposits as these surface processes (such as DNA binding, degradation or consumption) are part of the wider burial process. However, the small (0.22 g) samples were stored in isolation rather than as whole core samples and may not be affected by processes occurring in a larger sediment matrix, such as microbial recruitment from the adjacent sediment. Our interpretation is limited to lakes of similar trophic level and sediment type. The interplay of sediment binding sites and lower algal productivity in oligotrophic lakes may produce differing scenarios and further studies are required.