Metabarcoding Reveals Lacustrine Picocyanobacteria Respond to Environmental Change Through Adaptive Community Structuring

Picocyanobacteria (Pcy) are important yet understudied components of lake foodwebs. While phylogenetic studies of isolated strains reveal a high diversity of freshwater genotypes, little is known about abiotic drivers associated with Pcy in different lakes. Due to methodological limitations, most previous studies assess potential drivers using total cell abundances as a response, with often conflicting and inconsistent results. In the present study, we explored how picocyanobacterial communities respond to environmental change using a combination of epifluorescence microscopy and community data determined using 16S rRNA gene metabarcoding. Temporal shifts in picocyanobacterial abundance, diversity and community dynamics were assessed in relation to potential environmental drivers in five contrasting lakes over 1year. Cell abundances alone were not consistently related to environmental variables across lakes. However, the addition of metabarcoding data revealed diverse picocyanobacterial communities that differed significantly between lakes, driven by environmental variables related to trophic state. Within each lake, communities were temporally dynamic and certain amplicon sequence variants (ASVs) were strongly associated with specific environmental drivers. Rapid shifts in community structure and composition were often related to environmental changes, indicating that lacustrine Pcy can persist at high abundances through collective community adaptation. These results demonstrate that a combination of microscopy and metabarcoding enables an in-depth characterisation of picocyanobacterial communities and reveals strain-specific drivers. We recommend that future studies cease referring to picocyanobacterial as one functional group and take strain specific variability into consideration.


INTRODUCTION
Picocyanobacteria (Pcy) are increasingly being recognised as important components of freshwater ecosystems (Stockner, 1991;Caroppo, 2015). Forming the base of lacustrine microbial foodwebs, Pcy contribute significantly to primary productivity and their relative contribution is thought to increase with decreasing trophic state (Petersen, 1991;Stockner and Shortreed, 1991;Callieri and Stockner, 2002). Most work to date has focused on populations in oligo-and mesotrophic lakes, however recent studies suggest that the importance of Pcy in eutrophic and hypertrophic lakes may have been overlooked (Carrick and Schelske, 1997;Jasser et al., 2010;Wood et al., 2017;Schallenberg et al., 2021). Pcy are now known to inhabit a wide range of freshwater and brackish environments, likely due to their high adaptability and phenotypic plasticity (Callieri, 2017;Huber et al., 2017).
While the diversity and dynamics of marine Pcy (largely marine Synechococcus and Prochlorococcus) have been wellstudied, freshwater Pcy have received less attention (Callieri and Stockner, 2002). New technologies are revealing the presence and abundance of diverse Pcy strains in a range of freshwater environments (Ruber et al., 2016;Wood et al., 2017). However, little is known about the temporal diversity of Pcy across lakes of varying type and trophic state, and whether the environmental drivers that control Pcy dynamics differ among lakes.
Unravelling the complex abiotic drivers of species composition and density is difficult. To date, most studies exploring environmental drivers have used data on cell abundances and pigment type and not considered changes in community composition. These studies often produce inconsistent and conflicting results, adding to the uncertainty around which factors are important in controlling freshwater Pcy abundances and shaping communities. A review of studies which have assessed environmental drivers of Pcy community dynamics illustrates the difficulty in determining not only broad scale but also local population drivers ( Table 1). For example, within oligotrophic lakes, total phosphorus (TP) is both positively (Stockner and Shortreed, 1991) and negatively Zhong et al., 2013) correlated with Pcy abundances (Table 1). Similarly, within eutrophic and hypertrophic lakes, temperature has been shown to be both positively (Szelag-Wasielewska, 2004;Pulina et al., 2017) and negatively (Huber et al., 2017) correlated with Pcy abundance (Table 1). Therefore, it is not surprising that large-scale spatial studies, including those that assess lakes across a trophic gradient, often find insignificant correlations (e.g., Lavallée and Pick, 2002; Table 1).
This variation in identified drivers could be due to differences in the types and locations of lakes, although research on larger planktonic cyanobacteria, especially bloom forming cyanobacteria, suggests drivers are similar at a global scale (Harke et al., 2016;Huisman et al., 2018). A more likely scenario is that due to their small size and associated issues with identification beyond family level, Pcy have historically been studied and reported as one group ("Picocyanobacteria" or "autotrophic picoplankton"), failing to acknowledge the diversity and varied functional abilities of Pcy (Cabello-Yeves et al., 2018;Sánchez-Baracaldo et al., 2019;Di Cesare et al., 2020). While total cell abundances are important for understanding the contribution of Pcy to primary productivity, shifts or differences in community composition might better reflect changes in biotic and abiotic drivers.
Diverse communities of Pcy have now been found in freshwater and brackish lakes using genetic fingerprinting and barcoding techniques (Sanchez-Baracaldo et al., 2008;Callieri et al., 2012a;Ruber et al., 2016;Pulina et al., 2017;Wood et al., 2017;Schallenberg et al., 2021), suggesting that freshwater Pcy diversity is as high or higher than that of their marine counterparts (Sanchez-Baracaldo et al., 2008). Studies employing these techniques have found that different Pcy operational taxonomic units (OTU's) can respond differently to the same driver, indicating the potential for community level responses to drivers. For example, Felföldi et al. (2011) found different Pcy OTU's responded differently to drivers such as temperature (Table 1). Further, Bertos-Fortis et al. (2016) found that the most abundant Synechococcus OTU's identified in Baltic Sea samples correlated in often-opposite ways with the same environmental variable. Studies from the marine environment also highlight distinct niche specificity of Pcy strains across horizontal and vertical transects (Stockner and Shortreed, 1988;West and Scanlan, 1999;Johnson et al., 2006;Six et al., 2007). Callieri et al. (2012a) showed that this was also the case in Lake Maggiore, where distinct vertical partitioning of Pcy OTU's was found during stratification. With an often-greater diversity of strains inhabiting lacustrine compared to marine environments, it is possible that spatial and temporal community shifts may exhibit a stronger response to environmental change than overall cell abundance.
In the present study we used epifluorescence microscopy and community composition data from DNA metabarcoding (targeting the 16S ribosomal RNA gene; 16S rDNA) to assess temporal Pcy abundances, diversity and dynamics in five contrasting lakes over a 12-month period. This dataset allowed us to explore drivers of spatial and temporal changes in Pcy across a diversity of lake types while investigating the complementarity and effectiveness of these two methodologies. We hypothesised that: (i) Cell abundances within lakes will be temporally dynamic but not consistently related to environmental drivers across the lakes or trophic states, suggesting that strain-level responses may be occurring; (ii) Lacustrine Pcy will be diverse, with different communities inhabiting different lakes due to drivers related to trophic state; (iii) Pcy communities will differ temporally within lakes driven by shifts in physico-chemical variables; and (iv) Pcy communities in deep, stratified lakes will show vertical niche partitioning due to differing epi-and hypolimnetic conditions, suggesting vertical niche partitioning. Overall, we theorised that the combined use of microscopy and metabarcoding techniques will provide an enhanced picture of Pcy dynamics and responses to physicochemical change.

Study Lakes
Five lakes in southern New Zealand known to contain Pcy (Wood et al., 2017) were sampled monthly for 1 year (February 2018 to February 2019; Figure 1). These lakes represent different lake types, trophic states, morphologies, and salinities ( Table 2). Lakes Wanaka and Wakatipu are large, deep, oligotrophic, Frontiers in Microbiology | www.frontiersin.org typically monomictic peri-alpine lakes ( Table 2). Lake Hayes is a small, typically monomictic, eutrophic lake that experiences both clear and turbid water phases. Tomahawk Lagoon is a small, eutrophic, urban coastal lagoon without direct connection to the sea. Lake Ellesmere/Te Waihora is a large, supertrophic, Intermittently Closed and Open Lagoon Lake (ICOLL) which is opened to the sea multiple times per year ( Table 2).

Sample Collection
All lakes were sampled monthly in the main basin. In the monomictic lakes (Wanaka, Wakatipu, and Hayes), both epilimnion and hypolimnion samples were taken. Composite samples were taken from the epilimnion in Lakes Wanaka and Wakatipu that were comprised of equal volumes of water from 0.5, 10, 15, 30, and 45 m depths, and hypolimnion samples were taken from 150 m deep. In Lake Hayes, epilimnion samples were combined samples from 0.5, 2.5, 5, and 7.5 m depths while the hypolimnion samples were from 22 m. Samples from Lakes Wanaka, Wakatipu and Hayes were collected using a clean Van Dorn sampler before being transferred into sterile 1 L Nalgene™ bottles. In shallow Lake Ellesmere/Te Waihora and Tomahawk Lagoon, samples were taken from 0.5 m depth using a sterile, submerged 1 L Nalgene™ bottle. All samples were placed in the dark immediately and transported on ice to the laboratory within 24 h and subsampled for nutrient analysis, epifluorescence microscopy and DNA extraction.

Water Quality
Temperature (°C), dissolved oxygen (DO %) and salinity (ppt) depth profiles were recorded from Lakes Wanaka, Wakatipu, and Hayes using a SeaCAT CTD (SBE 19plus V2, Sea-Bird Electronics, WA, United States). The average value taken across the combined epilimnion depth range was used, to correspond with the depth of the epilimnion water sample taken. In Lake Ellesmere/Te Waihora and Tomahawk Lagoon, temperature, DO and salinity were measured using a water quality sonde (YSI, OH, United States), at 0.5 m depth. Secchi depth was measured at all lakes except for Tomahawk Lagoon. Nutrients measured were total nitrogen (TN), TP, dissolved reactive phosphorus (DRP), nitrate-nitrite nitrogen (NNN) and ammoniacal nitrogen (NH 4 -N). Samples (50 ml) for measurement of dissolved nutrients (NNN, NH 4 -N, and DRP) were filtered (GF/F, Whatman, United Kingdom) using an acid washed syringe. Nutrient samples from Lakes Ellesmere/Te Waihora, Wanaka, Wakatipu, and Hayes were analysed using a Continuous Flow Analyser at Hill Laboratories following APHA 4500 methods (Hamilton, New Zealand) with the following limits of detection; 1 μg L −1 (DRP and NNN), 4 μg L −1 (TP), 5 μg L −1 (NH 4 -N), and 10 μg L −1 (TN). Chlorophyll-a (Chl-a) analysis for these lakes was undertaken using subsamples of up to 500 ml of water following APHA 10200 H methods at the same laboratory with a detection limit of 0.00002 g/m 3 . Nutrients in Tomahawk lagoon samples were measured on a Continuous Flow Analyser (San++, Skalar Analytical, Breda, The Netherlands) using standard colorimetric methods with the following limits of detection; 0.1 μg L −1 (TN) and 0.5 μg L −1 (TP, DRP, NH 4 -N, and NNN). Samples for Chl-a analysis in Tomahawk Lagoon were filtered (<200 ml; GF/F, Whatman, United Kingdom) and the filters immediately frozen (−20°C) until later extraction using 90% alkaline aqueous acetone. Chl-a was measured spectrophotometrically following Wetzel and Likens (1991) using a UV mini 1,240 Spectrophotometer (Shimadzu, Japan).

Epifluorescence Microscopy
Cell abundances were enumerated using epifluorescence microscopy following the protocol described in Schallenberg et al. (2021). Briefly, 15 ml subsamples of water from Lakes Wanaka, Wakatipu, Hayes and Tomahawk and 0.1 ml from Lake Ellesmere/Te Waihora (diluted with 9.9 ml Milli-Q water) were each filtered onto black polycarbonate filters (0.2 μm Cyclopore, Whatman, United States) using a 1.2 μm backing filter to ensure even distribution. Filters were viewed using a Zeiss Axiophot microscope (Carl Zeiss Jena GmbH, Jena, Germany) equipped with a narrow-banded green filter set CY3 (Ex 520-570, DM 565, BA 535-675). Using this filter set, phycoerythrin-rich picocyanobacteria autofluoresce bright yellow while phycocyanin-rich picocyanobacteria autofluoresce bright red (Jasser and Callieri, 2017). Individual cells, colonies, and the number of cells per colony were counted from at least 20 fields of view or 400 individual cells at 1000× magnification.

DNA Extraction and 16S Ribosomal RNA Gene Sequencing
Water samples were filtered in triplicate (150 ml for Lakes Wanaka, Wakatipu, and Hayes; 30 ml for Lake Ellesmere/Te Waihora; 50 ml for Tomahawk Lagoon) using sterile filters (0.22 μm, S-Pak, Sigma, United States) and immediately frozen (−20°C). DNA extraction, PCR amplification and sequencing of Pcy sequences was carried out following the methodology described in Schallenberg et al. (2021). Briefly, DNA was extracted from filters using a PowerSoil ® DNA Isolation Kit (QIAGEN, United States) including one negative extraction control every 23rd sample. PCR was performed using cyanobacteria-specific primers CYB359-F 5'-GGGGAATYTTC CGCAATGGG-3' and CYB784-R 5'-ACTACWGGGGTATCT AATCCC-3' (Nübel et al., 1997) which target the V3-V4 region [~380 base pairs (bp)] of the cyanobacterial 16S rRNA gene. PCR conditions are outlined in Schallenberg et al. (2021). PCR products were purified using an Agencourt ® AMPure ® XP Kit (Beckman Coulter, CA, United States) before being diluted to 5 ng μl −1 and submitted to Auckland Genomics Limited (Auckland, New Zealand) for library preparation and sequencing. Library preparation followed the Illumina 16S metagenomics library prep manual. A second round of PCR added dual indexes to the samples and 5 μl of each subsequent PCR reaction was pooled before a single clean-up was undertaken. The pooled library was quality controlled using a bioanalyzer and the library was diluted to 4 nM, denatured, and diluted to a final concentration of 7 ρM with a 15% PhiX spike. An Illumina MiSeq was used to produce paired end sequences (2 × 250 bp) which have been deposited in the NCBI short read database under BioProject accession number: PRJNA749894.

Bioinformatics
Forward and reverse reads were assigned to samples after automatic demultiplexing using MiSeq Reporter (v2). Primers were eliminated from sequences allowing for one mismatch per sequence using Cutadapt (Martin, 2011). The remaining bioinformatic analysis was conducted in R version 4.0.4 (R Core Team, 2021). The DADA2 package version 1.16.0 (Callahan et al., 2016) was used to filter and trim forward and reverse reads to 230 and 228 bp with 2 and 4 maximum errors, respectively. Error profiles for the forward and reverse reads were assessed using a parametric model utilising 1 × 10 8 bases. Sequences were then dereplicated and a pseudo-pooling approach was used to infer amplicon sequence variants (ASVs). Sequences were merged, allowing one mismatch, with a minimum overlap of 10 bp followed by chimera removal using the default consensus approach. Taxonomy was assigned using the SILVA 138 database (Quast et al., 2013) by employing an RDP Naïve Bayesian Classifier with a bootstrap threshold of 50 (Wang et al., 2007). The phyloseq package version 1.32.0 (McMurdie and Holmes, 2013) was used for data pre-processing and statistical analysis. Chloroplasts and non-cyanobacterial sequences were removed, and the resulting cyanobacterial sequences filtered to include only ASVs with more than five reads that were present in more than one sample. This approach was chosen as low abundance ASVs were not the focus of this study and given the inclusion of triplicate samples, if an ASV was present in only one triplicate, it is likely to be erroneous. The distribution of sample sequencing depth at various processing steps is provided in Supplementary Figure S1 and a table of sequencing depths for each sample is provided in Supplementary Table S1. Triplicate samples were then merged after confirming they clustered together on a non-metric multidimensional scaling (NDMS) ordination.

Statistical Analysis
All statistical analysis was carried out using R version 4.0.4 (R Core Team, 2021) and Microsoft Excel (Excel version 2016). Bioinformatic and statistical analysis scripts can be found at https://github.com/lenaschall/Temporal_Picocyanobacteria_Scripts. Spearman rank correlations were calculated between total cell abundance and individual environmental variables for each sample. To test for significant differences in colonial cell contributions between the epi-and hypolimnion of the monomictic lakes, t-tests were conducted. Metabarcoding samples were normalised by randomly subsampling without replacement, to an even sequencing depth (4,000 reads) following assessment of a rarefaction curve plotted using the ggrare function from the ranacapa package (Kandlikar et al., 2018), confirming most samples had reached a plateau at this level (Supplementary Figure S2). Due to the low number of cyanobacterial reads in some samples, 24 samples and an associated 12 ASVs were discarded from the statistical analysis. For the remaining analysis of Pcy communities, rarefied data were used and only ASVs from the Order Synechococcales were retained.
Observed richness along with Shannon and Simpson alpha diversity metrics were calculated for each lake using the phyloseq package and differences in Shannon diversity between lakes was assessed using ANOVA. A permutational multivariate analysis of variance (PerMANOVA) was used to test for significant differences in picocyanobacterial communities between lakes Frontiers in Microbiology | www.frontiersin.org 6 November 2021 | Volume 12 | Article 757929 using the adonis function from the vegan package (Oksanen et al., 2020). The assumption of homogeneous dispersions between groups was tested between lakes using the betadisper function from the vegan package. To assess temporal community structure between months within each lake, reads were converted to relative abundances and rare ASVs (contributing on average <1% of sample reads) were combined into one category. Pcy community structure was analysed using stacked bar plots in the R package ggplot2 version 3.3.5 (Wickham, 2016). Canonical analyses of principal coordinates (CAP) were performed on the square-root transformed ASV data using both the presence/absence of ASVs (binary Jaccard) and total read abundance (weighted Jaccard) to assess differences in Pcy communities between the five lakes, and within each lake. For analysis between the five lakes, environmental variables related to all lakes were included together in the analysis, while the remaining analyses were conducted separately for each lake. Environmental variables were first analysed for collinearity and pairs returning a correlation coefficient >0.7 were scrutinised (Dormann et al., 2013), with the least biologically relevant variables removed from further analysis (Supplementary Table S2). The bioenv function from the package vegan was used to find the optimal subset of environmental variables that maximise the rank correlation with community dissimilarities. The Variance Inflation Factors (VIF) of the Bioenv selected variables were calculated to confirm variable independence before inclusion in the CAP model (Supplementary Table S2). The variables were then included as CAP constraints to determine the influence of the optimised set of environmental variables on Pcy community structure across months, using the weighted Jaccard dissimilarity matrices. Individual ASVs were plotted on the CAP, with abundant ASVs labelled, to assess relationships between ASVs and environmental constraints. To assess significance of the models, a PerMANOVA was calculated on the CAP ordinations using the Adonis function in vegan to determine if the constraints, environmental variables and CAP axes explained a significant amount of the variation in community dissimilarities.
To test for significant differences in the vertical structuring of Pcy communities in Lakes Wanaka and Wakatipu, a PerMANOVA was designed with depth as a factor. To determine if nutrient concentrations differed significantly with depth, nutrients measured from corresponding epi-and hypolimnetic samples from each monomictic lake were analysed using two-tailed t-tests.

Temporal Abundance and Environmental Correlations
Abundances of single celled and colonial Pcy were variable over the sampling period in all lakes (Figure 2;  Supplementary Table 3). In oligotrophic Lake Wanaka, peak abundances were recorded in spring (September, October, and November) in both the epi-and hypolimnion (Figure 2A). In oligotrophic Lake Wakatipu, bimodal autumn (March, April, and May) and winter (June, July, and August) peaks occurred in the epilimnion with peak hypolimnetic abundances occurring in Winter (August; Figure 2B). Winter destratification is reflected by a large increase in cell abundances in the hypolimnion (150 m depth) of Lakes Wanaka and Wakatipu in June and August, respectively (Supplementary Figure S3). Colonial cells in Lake Wanaka accounted for an average of 25% of total cells over the year, while colonial cells in Lake Wakatipu accounted for 7% on average. A t-test revealed no significant difference between the percent contribution of colonial cells in the epi-and hypolimnion of Lakes Wanaka (F = 0.54, p = 0.47) and Wakatipu (F = 2.15, p = 0.15).
In the two eutrophic lakes (Lake Hayes and Tomahawk Lagoon), both PE and PC-rich Pcy were found, with PE cells more common in Lake Hayes than Tomahawk Lagoon (Figures 2C,D). Pcy were not consistently present in these two lakes, with cell numbers in Lake Hayes peaking in Summer 2018 (January and February) but absent from July 2018 to January 2019 ( Figure 2C). No colonial cells were found in Lake Hayes. Pcy in Tomahawk lagoon reached a peak abundance in February 2018 (47,543 cells ml −1 ; Figure 2D). Tomahawk Lagoon had the highest variation in colonial cells, accounting for between 0 and 100% of total cells throughout the year. Supertrophic Lake Ellesmere/Te Waihora had a consistently high number of Pcy cells throughout the year, with a peak in February 2019 of 4.5 × 10 7 cells ml −1 (Figure 2E). The cells in Lake Ellesmere/Te Waihora consisted of single-cell and colonial PC-rich morphotypes, with colonial cells occurring in groups of four. Colonial cells contributed an average of 51% to the total cell count over the year, however a sharp decline in colonial cells occurred in February 2019 correlated with the largest peak in overall abundance ( Figure 2E).
Spearman's rank correlations between cell abundances (both total cells and the percentage of colonial cells) and corresponding water quality parameters revealed a number of significant correlations, both positive and negative (Table 3). No consistent correlations were found across all lakes or within oligotrophic or eutrophic lakes. Total hypolimnetic cell abundances in Lakes Wanaka and Wakatipu and colonial cells in Lake Wakatipu were significantly negatively correlated with NNN while in Lake Ellesmere/Te Waihora NNN was positively correlated with colonial cells (Table 3). In the Lake Wakatipu epilimnion, total cell abundances were negatively correlated with TP but positively correlated with DRP. The contribution of colonial cells was significantly positively correlated with TN and DO in Tomahawk Lagoon and Lake Ellesmere/Te Waihora, respectively.

Diversity and Dynamics of Picocyanobacterial Communities Between Lakes
A total of 13,996 ASVs were recovered across the five study lakes, 5,946 (42%) of which were cyanobacterial (class: Oxyphotobacteria). After removal of chloroplasts, filtering, and sample rarefaction, 206 cyanobacterial ASVs remained of which 68 ASVs (33%) were picocyanobacterial (order:  Figure 3). Lakes Wanaka, Wakatipu, and Hayes returned the lowest alpha diversity scores, regardless of the weighting imposed by the differing methods (Figure 3). Temporal shifts in Pcy observed richness occurred within each lake, with Lake Ellesmere/Te Waihora and Tomahawk lagoon exhibiting a decrease in richness over the sampling period (Supplementary Figure S4). Richness remained relatively stable in the other lakes except for the  Lake Hayes hypolimnion where there was a large increase in richness during January and February 2019 (Supplementary Figure S4). After filtering of potentially erroneous low-abundance ASVs and sample rarefaction, 68 Pcy ASVs were found across the five study lakes. Pcy communities were significantly different between the lakes using a weighted Jaccard dissimilarity matrix (PerMANOVA F = 31.51, R 2 = 0.73, p < 0.001; Figure 4). However, variances were not homogeneous between lakes (F = 2.68 p = 0.039), meaning within-lake community differences could in-part be driving this significant result. Nevertheless, CAP showed that samples from individual lakes clearly clustered together except for Lakes Wanaka and Wakatipu where communities showed some overlap (Figure 4). These patterns were consistent using dissimilarity matrices incorporating both the presence and absence of ASVs (binary Jaccard) and read abundance (weighted Jaccard; Supplementary Figure S5). Environmental conditions varied between lakes, as expected given the gradient of trophic states sampled (Supplementary Figure S6). Of note were the low DO concentrations recorded in the hypolimnion of Lake Hayes.
The constrained variables explained 59% of variation in community dissimilarity and a PerMANOVA of the model found the constraints to be highly significant (F = 14.70, p < 0.001) along with the individual terms DRP, NNN, and TP (p = 0.005 for each variable). CAP1 which explains 38.6% of the variation is related to trophic state, with TP and the removed collinear variables TN, Secchi depth and Chl-a associating strongly with that axis.

Oligotrophic Lakes
The Pcy community varied temporally in Lake Wanaka, with shifts in both the relative abundance of ASVs ( Figure 5A) and the presence and absence of abundant ASVs throughout the year. Five hypolimnion samples and a single epilimnion sample produced insufficient reads and were removed after rarefaction ( Figure 5A). Community composition and structure was highly similar in the epi-and hypolimnion from June 2018 to February 2019 when comparable samples were available, except for ASV CY.19 which remained abundant in the hypolimnion for 2 months (January and February 2019) after becoming rare (<1% of sample reads) in the epilimnion ( Figure 5A). A general decrease in the relative abundance of CY.1 and increase in CY.2 is seen in Winter and Spring (after lake turnover) in the epilimnion. In total, the environmental constraints explained 42% of the variation in the communities ( Figure 5B). A PerMANOVA conducted on the constrained model showed that the constraints were significant (F = 3.04, p = 0.004), the horizontal axis CAP1 was significant (F = 4.81, p = 0.010), and Secchi depth was a highly significant constraint (F = 3.23, p = 0.010), followed by TP (F = 2.43, p = 0.025).
In oligotrophic Lake Wakatipu, six samples (three from each the epilimnion and hypolimnion) were removed during rarefaction. Remaining samples showed variation in both community composition and structure across the study months ( Figure 5C). Relative abundances of ASVs CY.1 and CY.2 appear to shift in opposing ways, with lower relative abundances of CY.1 after

A B
FIGURE 3 | Shannon (A) and Inverse Simpson (B) alpha diversity estimates for the five study lakes, including samples from February 2018 to February 2019. Grey numbers above the plots indicate the total number of picocyanobacterial amplicon sequence variants (ASVs) in each lake over the 12-month period, prior to filtering and normalisation.
lake turnover, particularly in the hypolimnion. Community composition was otherwise similar in the epi-and hypolimnion over the months sampled, with temporal variations in structure seen across both layers, appearing to correspond with lake turnover in August ( Figure 5C; Supplementary Figure S3). A constrained ordination using weighted Jaccard dissimilarities showed that environmental constraints explained 31% of the variation in Pcy communities ( Figure 5D). However, PerMANOVA on the constrained model found the constraints to be insignificant in explaining community dissimilarities. In addition, most abundant ASVs cluster together near the origin of the CAP plot, further confirming that communities are similar across samples and the constraints are not sufficient at capturing the potential small amounts of variation observed.

Eutrophic Lakes
Collected samples from July 2018 to January 2019 in both the epi-and hypolimnion of Lake Hayes did not contain sufficient cyanobacterial reads and were removed after rarefaction along with the June hypolimnion sample. The resulting samples showed that the Pcy community in Lake Hayes was highly dynamic. The Pcy community in Autumn (May and June) was dominated by CY.1 with a shift to a CY.2 dominated community in the Summer (January and February; Figure 6A). Annual differences occurred in the epilimnion with CY.11 relatively abundant in January and February 2019 yet absent in February 2018. Epilimnetic communities also differ between May and June with CY.16, CY.31 and CY.48 abundant in May while in June CY.3 and CY.43 were abundant (Figure 6A). A PerMANOVA on the model found the constraints to be significant in explaining the community variation (F = 3.48 and p = 0.008), with Secchi depth highly significant (F = 4.57 and p = 0.005). Higher Secchi depths were related to the Autumn community, along with individual ASVs CY.2, CY.3, and CY.43, while CY.11, CY.18, and CY.4 were related to higher temperatures and CY.1 was related to lower temperatures ( Figure 6B). However, Secchi depth was positively significantly correlated with TP and NH 4 -N, and temperature was positively significantly correlated with DO and negatively with DRP, suggesting that these variables may also be involved in the resulting community structure.
Tomahawk Lagoon was not sampled in February and November 2018 or February 2019. Samples from March, April, and September to November 2018 had insufficient cyanobacterial reads and were removed from the analysis post-rarefaction. The Pcy community was relatively similar in Tomahawk Lagoon over the months when Pcy were present, with the exception of January 2019 when the ASV CY.11 became highly abundant ( Figure 6C). This sampling coincided with a large bloom of Aphanizomenon sp. and Dolichospermum sp. (Supplementary Figure S7), which together accounted for 96% of cyanobacterial reads. Lower-abundance ASVs (accounting for <10% of total reads) were variable throughout the year, with increases of CY.4 in July and CY.45 in August. The environmental constraints explained 67% of the variation (Figure 6D), however there were too few samples to confidently run a permutational analysis of the model. Nevertheless, both the bar plots and CAP suggest that the community in January 2019 diverges from other sampled months, correlated with high TN and associated with a large increase in the relative abundance of ASV CY.11. Interestingly, the increased abundance of CY.11 in Summer 2019 is seen in both eutrophic lakes.
The Pcy community in Lake Ellesmere/Te Waihora showed temporal succession between February 2018 and December 2019 with a gradual increase in CY.3 and CY.5 ( Figure 7A). Simultaneously, CY.0 and CY.12 decreased in relative abundance. An abrupt community shift occurred between December 2018 and January 2019 correlating with a lagoon opening event that took place from 12 December 2018 until 25 January 2019. DO and salinity explained 51% of the variation in community dissimilarities and were found to be significant (F = 4.67 and p = 0.001; Figure 7B). Despite three smaller lagoon opening events over the year, salinities remained between 5 and 8 ppt from February 2018 to December 2018 before rising to 13.5 ppt in February 2019. Subsequently, the community became heavily dominated by CY.0, with CY.12 and CY.9 also contributing high relative abundances. DO levels were highest from July to December 2018 (>100%), associated with ASVs CY.3, CY.5, CY.7, CY.10. CY.13, CY.20, and CY.21. ASVs CY.24 and CY.30 were abundant in February to April 2018, associated with lower DO levels.

Vertical Picocyanobacterial Dynamics in Monomictic Lakes
T-tests confirmed NNN (Wanaka: t = 2.08, p = <0.01, Wakatipu: t = 2.16, p = >0.01) and TN (Wanaka: t = 2.12, p = 0.01, Wakatipu: t = 2.08, p = >0.01) differed significantly between the epi-and hypolimnion of Lakes Wanaka and Wakatipu during the sampling period while other nutrients were similar with depth. Similarly, PAR differed between sampled depths in both lakes, with 1% surface PAR being reached between 22 and 39 m depths in both lakes over the course of the year. However, PerMANOVA analysis (factor = depth) confirmed communities were not significantly different between the two vertical layers in either lake (Wanaka: F = 1.07 and p = 0.40, Wakatipu: F = 1.44 and p = 0.11).
Nutrient concentrations between the epi-and hypolimnion in eutrophic Lake Hayes were significantly different over the sampling period (DRP: t = 2.08, p = <0.01, TP: t = 2.12, p = 0.04, and NH 4 -N: t = 2.09, p = 0.02). While a permutational test could not be performed due to the low number of samples, interesting differences in community structure and composition were observed between comparable epi-and hypolimnion samples ( Figure 6A). For example, in May, 1 month prior to lake turnover, the epilimnion contained ASVs (CY.16,CY.29,CY.31,and CY.48) not present in the hypolimnion sample. Conversely, during strong stratification in January and February 2019, CY.22 and CY.31 were abundant in the hypolimnion yet not present in the epilimnion (Figure 6A). During January and February 2019, DO at the hypolimnion sampling depth was <1.5%, indicating that these ASVs may be able to tolerate severe hypoxia.

Temporal Abundances and Environmental Correlations
Picocyanobacterial abundances differed between and within lakes, with supertrophic Lake Ellesmere/Te Waihora consistently returning the highest cell abundances. In general, abundances in the oligotrophic lakes were less variable compared to the eutrophic lakes. Trophic state appeared to be related to pigment type and overall temporal patterns of abundance. PE-rich cells were dominant in oligotrophic lakes while eutrophic lakes contained both PE and PC cells at differing ratios and supertrophic Lake Ellesmere/Te Waihora contained only PC cells. Previous studies support these findings (Vörös et al., 1998;Vila and Abella, 2001), however eutrophic lakes have typically been found to contain predominantly PC cells, whereas the present study found eutrophic Lake Hayes to be dominated by PE cells. Lake Hayes is a complex lake known to fluctuate between clear and turbid water states and is showing signs of slowly recovering from eutrophication (Schallenberg, 2020). This variable water quality with periods of high water clarity may explain the higher than predicted proportion of PE-rich cells found.

Single Cell Abundances and Environmental Correlations
Studies have found that single or bimodal abundance peaks occurring in Spring or early Summer and Summer-Autumn are common for temperate oligotrophic lakes while abundances in eutrophic lakes tend to peak during Summer (Callieri et al., 2012b). The present study supports this, with a single Spring peak and bimodal Autumn and Winter peaks recorded in the epilimnion of oligotrophic lakes Wanaka and Wakatipu, respectively. In these lakes, Pcy abundances in the hypolimnion were related to lake turnover with abundances increasing significantly once the lakes became isothermal. In Winter, cell abundances were higher in the hypolimnion than the epilimnion in Lake Wanaka. Interestingly, NNN was negatively correlated with total cell abundances in the hypolimnion of Lakes Wanaka and Wakatipu, while TP was negatively correlated with Wakatipu epilimnion abundances. Negative responses of Pcy to nutrient addition (N and P) have been previously reported in oligotrophic lakes and mesocosms (Rhew et al., 1999;Burns and Schallenberg, 2001;Schallenberg and Burns, 2001;Burns and Galbraith, 2007;Yang et al., 2019), however the reasons for this response remain unknown and require further study. Total cell abundances in the eutrophic and supertrophic lakes were not significantly correlated with any measured environmental variables. In the small, eutrophic lakes (Hayes and Tomahawk), strong Summer peaks were followed by low abundances and sometimes absences of Pcy in other seasons and months, a common pattern that has been reported from other eutrophic systems (Callieri et al., 2012b). Conversely, supertrophic Lake Ellesmere/Te Waihora had consistently high abundances throughout the year. While it is generally accepted that Pcy may reach high abundances in eutrophic lakes, their relative importance and contribution to total phytoplankton biomass is thought to decrease with increasing trophic state (Stockner, 1991;Bell and Kalff, 2001). Often, cell abundances decrease with increasing trophic state, accompanied by an increase in picoeukaryotes (Ruber et al., 2018). Pick (2000) showed that the relationship between cell abundances and TP or trophic state is not necessarily linear, instead following a polynomial function whereby Pcy abundances increase linearly up to TP concentrations of ~15 μg L −1 before declining. However, some supertrophic lakes and coastal lagoons such as Lake Ellesmere/ Te Waihora (average TP over the study period = 200 μg L −1 ) do not appear to adhere to this model, with abundances often high and particularly difficult to predict (Bec et al., 2011;Callieri et al., 2012b;Caroppo, 2015). Microbial studies in brackish, eutrophic lakes are increasingly revealing significant Pcy populations which may have been previously overlooked (Carrick and Schelske, 1997;Pulina et al., 2017;Schallenberg et al., 2021), demonstrating the need to improve our knowledge of microbial foodwebs in these systems. This is particularly pertinent with sea level rise and climate change predicted to significantly affect these habitats (Carrasco et al., 2016).

Colonies and Environmental Correlations
The relatively low number of colonial cells found in the oligotrophic lakes is not surprising given that eutrophic environments tend to favour these morphotypes (Callieri et al., 2012b). However, high numbers of colonies were recorded in certain months, e.g., Lake Wanaka in February 2019 >50% of cells formed colonies. Many hypotheses exist around the drivers of colony formation in temperate oligotrophic lakes including possible ease of nutrient transfer in times of nutrient depletion, UV protection and predator avoidance (Klut and Stockner, 1991;Crosbie et al., 2003;Callieri et al., 2016). Colonial peaks in the epilimnion occurred in Summer in both oligotrophic lakes when UV levels are typically highest. The February 2019 colonial peak in Lake Wanaka coincided with a 50% drop in NNN concentrations, while colonial cells in the hypolimnion of Lake Wakatipu were significantly negatively correlated with NNN. Interestingly, no colonies were found in eutrophic Lake Hayes during this study. In contrast, colonies in eutrophic Tomahawk Lagoon and supertrophic Lake Ellesmere/Te Waihora were common, accounting for up to 100% of Pcy cells suggesting that colonial morphotypes are commonly associated with brackish eutrophic coastal lakes and lagoons. Colonial cells in Tomahawk lagoon and Lake Ellesmere/Te Waihora were positively correlated with certain nutrients. This contrasts with the negative relationship found in Lake Wakatipu, indicating that Pcy in these lakes respond differently to drivers relating to trophic state. The abundance of colonial cells in Lake Ellesmere/Te Waihora dropped drastically in February 2019, 2 months after the lake was opened to the sea, implying that significant mixing with the marine environment may have reduced colonial morphotypes.
In line with our first hypothesis, cell abundances were shown to be temporally dynamic within the study lakes. This was also the case in both the epi-and hypolimnion of monomictic lakes. While cell abundances were correlated with environmental variables, responses were inconsistent between lakes and no environmental variables were significantly correlated with abundances across all lakes, or within trophic states.

Diversity and Dynamics of Picocyanobacterial Communities Between Lakes
Pcy communities were diverse in all study lakes and these communities differed significantly between most lakes, confirming our second hypothesis. Shannon diversity differed significantly between the lakes and overall patterns of diversity were consistent using both the Shannon and Inverse Simpson alpha diversity indices which are weighted towards rare and abundant species, respectively. Observed richness was highest in supertrophic lake Ellesmere/Te Waihora and varied temporally in the supertrophic and eutrophic lakes. The five lakes revealed diverse and distinct Pcy communities, except for oligotrophic Lakes Wanaka and Wakatipu which shared five abundant ASVs. These lakes are similar in trophic state, type, size and proximity ( Table 2). However, the results of this study suggest that proximity is less important in determining Pcy community structure than local conditions and trophic state. For example, eutrophic Lake Hayes and oligotrophic Lake Wakatipu are situated with 15 km of each other yet contain significantly different Pcy communities. The physico-chemical variables TP, NNN and DRP were identified as significant drivers differentiating Pcy communities among the five study lakes. These drivers along with their removed cross-correlated variables are related to lake trophic state, suggesting that trophic state may be a significant driver of Pcy community structure across lakes. The present findings contrast somewhat with those of Wood et al. (2017) who found several dominant Synechococcus OTU's persisted in lakes across a range of trophic states. Similarly, Ruber et al. (2016) found that Pcy communities were similar in trophically contrasting, connected lakes. One possible reason for this could be the use of ASVs in the present study while Wood et al. (2017) and Ruber et al. (2016) used OTU clustering, potentially allowing for numerous ASVs to group as one OTU. The connectedness of the lakes in Ruber et al. (2016) could also explain the similarities found, suggesting that dispersal is stronger than local selection in connected lakes. It is also likely other factors not included in this study are responsible for the remaining unexplained variation and a larger dataset including community analysis in lakes across the full trophic gradient might reveal more detailed community patterns.

Temporal Drivers of Picocyanobacterial Communities in Contrasting Lakes
Picocyanobacterial communities within lakes varied temporally and these shifts were often related to changing environmental conditions. Lakes Wanaka and Wakatipu communities revealed a seasonal succession of abundant ASVs which was related to Secchi depth and TP. Although the potential contribution of Chl-a and DRP in shaping these communities should not be ignored given these variables were removed due to significant negative and positive correlations with Secchi depth, respectively. Similarly in Lake Wakatipu, Secchi depth and TP were among the best predictors of community dissimilarities, however these were not statistically significant. ASVs shared between the two oligotrophic lakes were abundant at the same times of year except for CY.1 which peaked in relative abundance in Summer in Lake Wanaka while being a substantial component in Lake Wakatipu in Summer, Autumn, and Winter. Across both lakes, CY.1 and CY.2 alternate as being the most abundant ASV present. Phylogenetic analysis by Schallenberg et al. (2021) suggests that these ASVs are not closely related, belonging to Group A Cyanobium gracile and Clade III, respectively.
Temporal analysis of eutrophic Lake Hayes suggests a highly dynamic and responsive Pcy community, with large temporal shifts in the presence/absence and relative abundance of ASVs. These shifts occurred both between seasons and within seasons, related to Secchi depth and temperature. Lake Hayes is a eutrophic, monomictic lake that often experiences hypoxia in the epilimnion and anoxia in the hypolimnion over Summer (Bayer et al., 2008). While DO was not included in the CAP analysis due to multiple instances of collinearity with other variables, it was positively correlated with temperature suggesting that it is also related to community structuring. These results suggest niche specificity of certain Pcy ASVs in Lake Hayes as well as a rapid adaptation to changing environmental conditions. Pcy community structure in Tomahawk Lagoon was relatively stable between seasons until a major shift in Pcy dominance in January 2019. Within 1 month, two previously abundant ASVs disappeared and one ASV became highly abundant. This shift was related to a bloom of Dolichospermum sp. and Aphanizomenon sp., suggesting that co-occurring cyanobacteria may cause strain-specific Pcy responses. Certain Pcy have been found associated with blooms of larger filamentous cyanobacteria (Pick, 2016)  The Lake Ellesmere/Te Waihora Pcy community underwent gradual temporal succession for most of the year followed by a large shift in community structure after the ICOLL was opened to the ocean in December 2018. The corresponding increase in salinity was reflected by a major shift in Pcy dominance, suggesting halotolerance of certain ASVs. Analysis of corresponding environmental drivers confirmed that these shifts were related to salinity and DO levels. Given the highly dynamic nature of ICOLLS including temporal shifts in water quality, salinity and water level (Schallenberg et al., 2010), Pcy in these systems are subject to dramatic environmental change. The dominance of Pcy in some hypertrophic lakes and lagoons (Sorokin and Dallocchio, 2008;Schapira et al., 2010;Bec et al., 2011;Albrecht et al., 2017;Pulina et al., 2017) might be explained by the high adaptability of freshwater Pcy (Callieri, 2017;Huber et al., 2017). The present results suggest that temporal succession of different ASVs could explain this adaptability, allowing Pcy to persist at high abundances while community composition varies in response to environmental change.
Of the potential physico-chemical variables included in this study none were consistent in driving communities across all study lakes. This is not surprising given the different trophic conditions and individual ASVs dominating the Pcy communities. In support of our third hypothesis, we have shown that within-lake communities are temporally dynamic, and abundant ASVs are associated with specific environmental drivers. Different species were shown to fill specific niches as environmental conditions changed. In addition, certain ASVs were found to persist at rare levels (<1% of sample Pcy reads) before significantly increasing in relative abundance within months. Teasing apart the many potential abiotic drivers of Pcy is difficult due to the cross-correlated nature of variables and the effects of seasonality (Callieri et al., 2012b). However, we found patterns of gradual succession and rapid disturbance following disruptive events such as mixing in deep lakes and marine influence in a coastal lagoon. It is likely that drivers other than those measured in this analysis may explain some of the variation in community structure. For example, Callieri et al. (2016) and Huber et al. (2017) found that grazing pressure by zooplankton can exert preferential pressure on certain Pcy phenotypes. Future studies would benefit from incorporating potential biotic drivers which may help explain Pcy abundances and community structure in freshwater systems.

Picocyanobacteria Community Dynamics in Deep Stratified Lakes
Vertical partitioning of Pcy strains is common in the oligotrophic marine environment where significant temperature and light gradients create niches utilised by various clades of marine Pcy (Six et al., 2007;Scanlan et al., 2009;Bertos-Fortis et al., 2016;Farrant et al., 2016). Such depth partitioning was also found in Lake Maggiore where Pcy OTU's differed with depth, particularly during Summer stratification (Callieri et al., 2012a) and high altitude Andean lakes (Caravati et al., 2010). However, no significant vertical partitioning of Pcy communities was found in monomictic, oligotrophic Lakes Wanaka and Wakatipu over the course of the year. While community structure varied temporally, this structure was mirrored in both the epi-and hypolimnion of these lakes. This contrasts with findings from Lake Maggiore, a similarly large peri-alpine, oligotrophic lake, however these results are similar to those found by Ruber et al. (2018) in the Osterseen lakes of Bavaria. There, relatively small differences in Pcy OTU's were found with depth while temporal community shifts were greater (Ruber et al., 2018). While TN and NNN differed significantly with depth in these study lakes, along with PAR, it is possible that these differences are not strong enough to drive community shifts in these lakes. In the present study, a lack of comparable epi-and hypolimnetic samples during both isothermal and mixed periods, meant that that it was not possible to use a nested analysis of both depth and mixing phase to determine the influence of mixing phase on vertical communities.
In contrast to these large, monomictic, oligotrophic lakes, vertical partitioning of the Pcy community was strong in eutrophic, monomictic Lake Hayes although due to the limited sample size this could not be tested statistically. Physico-chemical conditions differed significantly between the epi-and hypolimnion over the study period. In all 3 months with comparable epiand hypolimnion samples, communities differed in both structure (presence/absence) and composition (relative abundance). Of particular interest was the presence of ASVs in the hypolimnion of Lake Hayes during Summer hypoxic conditions. This indicates that Pcy may not only survive under severe hypoxia, but certain ASVs that appear and become relatively abundant during this time may have a selective advantage under these conditions. While DNA metabarcoding alone cannot determine whether cells are alive and viable, survival under hypoxic and anoxic conditions may not be uncommon for Pcy, with the recent discovery of viable Pcy at 750 m depth in the anoxic Black Sea . The strains isolated from this depth had the ability to accumulate Chl-a under dark anoxic conditions and photosynthesise when exposed to oxic, light conditions . This finding provides further evidence of picocyanobacterial ubiquity and challenges our understanding of Pcy physiology by highlighting the extreme environments inhabited by certain strains.
In contrast to our fourth hypothesis, these results indicate that vertical partitioning of Pcy communities in deep lakes is not always present, can be variable and may be lake-specific. We recommend the use of larger datasets to statistically confirm differences between stratified and mixed periods, along with culturing of strains of interest from both the epi-and hypolimnion of stratified lakes to assess optimal growth conditions. This is of particular interest in low-oxygen environments such as Lake Hayes, where differing photosynthetic and metabolic pathways could be used. The use of metagenomics and metatranscriptomics would provide further information on the functions of these strains and how their functional profiles shift with changing environments.

CONCLUSION
The temporal and vertical dynamics of freshwater Pcy were explored in five lakes with contrasting trophic states using a combination of epifluorescence microscopy and environmental DNA metabarcoding. Abundances differed between and within lakes, and there were few relationships with environmental variables. However, these responses were not consistent across all lakes. The addition of community-scale analysis revealed Pcy communities were diverse and differed significantly between lakes. Pcy communities within the lakes were temporally dynamic and certain ASVs were strongly associated with lake-specific environmental drivers. Shifts in community composition were often related to environmental changes, indicating that lacustrine Pcy respond rapidly to change through adapting community structure. This study shows that community-scale responses allow Pcy to persist at high abundances through seasonal changes and disruptive events including isothermal mixing in deep lakes and a marine intrusion in a coastal lagoon. Finally, vertical structuring of Pcy communities was found only in one eutrophic lake, suggesting vertical niche partitioning may not be as common in lakes as in the marine environment. Together these results demonstrate that a combination of microscopy and DNA metabarcoding can lead to an improved understanding of freshwater Pcy dynamics and reveal strainspecific drivers. We recommend that future studies cease referring to Pcy as one functional group and utilise molecular techniques to allow for the identification of individual genotypes and their specific drivers.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: https://www.ncbi.nlm. nih.gov/bioproject/749894. Bioinformatic and statistical analysis scripts can be found at: https://github.com/lenaschall/ Temporal_Picocyanobacteria_Scripts.

AUTHOR CONTRIBUTIONS
LS, SW, and CB conceived and designed the study. LS analyzed the data, prepared figures and tables, and authored the manuscript draft. JP and SW helped with metabarcoding and statistical analysis. All authors contributed to the article and approved the submitted version.