Amplicon-based and metagenomic approaches provide insights into toxigenic potential in understudied Atlantic Canadian lakes

Samples were collected bi-weekly from three sites around each of Harvey Lake and Washademoak Lake from 30 July to 10 November 2018 (Fig. 1). Harvey Lake is located in York County, New Brunswick, and is approximately 6 km long and 1.25 km wide, with a maximum depth of 13.6 m. It is classified as an oligotrophic lake, based on a median total phosphorus concentration of 5 μg·L⁻¹ from 2005 to 2016 (New Brunswick Department of Environment and Local Government 2019). It is included in the Government of New Brunswick fish-stocking program, receiving landlocked Atlantic salmon. Washademoak Lake is a widening of the Canaan River upstream of the Saint John River in Queen’s County, New Brunswick. It has a maximal depth of 28.7 m, and is approximately 22 km long and 2.5 km at the widest point. It is mesotrophic, based on a median total phosphorus concentration of 17 μg·L⁻¹ from 2005 to 2016 (New Brunswick Department of Environment and Local Government 2019).

Surface water temperature was recorded, and surface samples were collected in 4-L jugs by wading, then transported back to the lab at in situ temperatures for immediate processing. Maximal volumes of each sample were filtered through 3-μm nominal pore-size polycarbonate filters (MilliporeSigma, Burlington, MA, USA) (Table S1). Total DNA was extracted from filters using MoBio’s PowerWater DNA Isolation Kit (now QIAGEN, Germantown, MD, USA) following the manufacturer’s instructions. Extracts were eluted in 100 μL Tris and stored at −20 °C until further analysis.

Three strains of cyanobacteria were used as positive controls for the PCR-based assays. Microcystis aeruginosa CPCC 300 (Canadian Phycological Culture Centre, Waterloo, Canada) was grown on a 14:10 light–dark cycle in BG11 medium at 20 °C. Dolichospermum lemmermanii NIVA-CYA83/1 (syn. Anabaena lemmermanii NIVA-CYA83/1) and Planktothrix agardhii NIVA-CYA126 (Norwegian Institute for Water Research, Oslo, Norway) were grown on a 14:10 light–dark cycle in Z8 medium at 18 °C. Cultures were either extracted using MoBio’s Microbial DNA Isolation Kit and eluted in 100 μL of TrisEDTA (Microcystis), or MilliporeSigma’s Gene Elute Bacterial Extraction Kit and eluted in 200 μL of TrisEDTA (Dolichospermum and Planktothrix). DNA extracts were stored at −20 °C until further analysis.

The presence/absence of the general microcystin synthetase E gene (mcyE) was determined using PCR-based detection with mcyE-F2 (Vaitomaa et al. 2003) and mcyE-R4 (Rantala et al. 2004) primers (Table S2). All reactions were performed in 20-μL volumes with 1 μL of DNA template per reaction, and positive (appropriate culture extracts) and negative (ddH₂O) controls were included in each PCR run. For the initial screening with mcyE-F2 and mcyE-R4, the PCR mixture contained 1 × Mg (+) PCR buffer, 250 μM dNTPs, 0.5 μM of forward and reverse primers, 0.5 U of Ex Taq DNA (Takara Bio, San Jose, CA) polymerase and sterile water. An MJ Mini Personal Thermal Cycler (BioRad, Hercules, CA, USA) was used for the PCRs cycles with initial denaturation step at 95 °C for 3 min, followed by 24 cycles of 94 °C for 30 s, 60 °C for 30 s, and 72 °C for 1 min, and a final extension at 72 °C for 10 min. Samples that tested positive for mcyE were then assessed with genus-specific assays to identify the source of the toxin genes. For Microcystis-specific mcyE, 20-μL PCRs were set-up with 1X Ex Taq Buffer Mg2+ plus (Takara Bio), 250 μM dNTPs, 0.5 μM of mcyE-F2 and MicmcyE-R8 primers (Table S2), 1.25 U of Ex Taq polymerase (Takara Bio), and sterile water (Vaitomaa et al. 2003). The PCR amplification included an initial denaturation step at 95 °C for 3 min, followed by 38 cycles of 94 °C for 30 s, 60 °C for 30 s, 72 °C for 60 s, and a final extension at 72 °C for 10 min (Vaitomaa et al. 2003). For Planktothrix-specific mcyE, PCRs were set-up with 1X Standard Taq Reaction Buffer (New England Biolabs, Ipswich, MA, USA), 62.5 μM dNTPs, 0.5 μM of mcyE-F2 and PlamcyE-R3 primers (Table S2), 1.25 U of Taq DNA polymerase (New England Biolabs), and sterile water to a final volume of 20 μL (Rantala et al. 2006). The PCR amplification was as follows: initial denaturation step at 95 °C for 3 min, followed by 35 cycles of 94 °C for 30 s, 55 °C for 45 s, 68 °C for 90 s, and a final extension at 68 °C for 5 min (Rantala et al. 2006). The Dolichospermum-specific mcyE reaction was set-up in 20 μL 1X Ex Taq Buffer Mg2+ plus (Takara Bio Clontech), 250 μM dNTPs, 5.0 μM of mcyE-F2 and 0.5 μM AnamcyE-R12 primers (Table S2), 1.25 U of Ex Taq polymerase (Takara Bio), and sterile water (Doblin et al. 2007). A touchdown PCR program was used, consisting of an initial denaturation step at 98 °C for 3 min, followed by 25 cycles of 94 °C for 30 s, 65 °C for 30 s decreasing by 0.5 °C per cycle, 72 °C for 60 s, and a final extension at 72 °C for 5 min (Doblin et al. 2007). PCR products were run at 70 V for 40 min on 1–1.5% TrisBorateEDTA agarose gels, stained with 1X SYBRSafe Gel stain, (Invitrogen, Waltham, MA, USA), and visualized using a BioRad gel documentation system.

The presence/absence of the anatoxin-a biosynthesis C gene (anaC) was assessed by PCR amplification with the anaC-genF and anaC-genR primers (Table S2; Rantala-Ylinen et al. 2011) as follows: a 20 μL PCR was set-up with 1 μL of DNA, 1X Ex Taq Buffer Mg2+ plus (Takara Bio), 200 μM dNTPs, 0.5 μM of anaC-genF and anaC-genR primers (Table S2), 1.25 U of Ex Taq polymerase (Takara Bio), and sterile water (Rantala-Ylinen et al. 2011). PCR conditions for this reaction started with an initial denaturation step at 94 °C for 2 min, followed by 35 cycles of 94 °C for 30 s, 55.5 °C for 30 s, 72 °C for 30 s, and a final extension at 72 °C for 5 min (Rantala-Ylinen et al. 2011).

All samples negative for mcyE and anaC were assessed for PCR inhibition using 16 S rRNA gene amplification with the universal bacterial primers 27 F and 1492 R (Lane 1991). The PCR reaction mix (25 μL) consisted of: 1 μL of extracted genomic DNA, 1x reaction buffer, 2 mM MgCl₂, 0.2 mM dNTP, 1 μM of each primer, Taq polymerase (1 U, New England Biolabs) and sterile water. The cycling conditions were as follows: initial denaturation at 94 °C for 2 min followed by 30 cycles of 94 °C for 1 min, 55 °C for 2 min, and 72 °C for 2 min. The final cycle was followed by extension at 72 °C for 5 min.

To investigate seasonal changes in the microbial community composition and its potential to produce cyanotoxins and other bioactive metabolites, we selected 8 DNA samples from each lake (Harvey (site A) and Washademoak (site B)) for shotgun Illumina sequencing (see Table S1 for sampling dates details). The 16 libraries (∼450 base pairs insert length) were prepared with the Illumina Nextera Flex kit (Illumina, San Diego, CA, USA) for MiSeq+NextSeq platforms and sequenced (150 base pairs paired-end reads) at the Integrated Microbiology Resource of Dalhousie University (imr.bio/) using the Illumina NextSeq 550 system.

Raw Illumina reads were uploaded to MG-RAST 4.0.3 (Meyer et al. 2008) (mg-rast.org) for quality control. Then, reads were assigned to taxonomic units using BLAT (Kent 2002) and the M5NR database (Wilke et al. 2012) as reference, considering 60% sequence similarity in a protein sequence (conceptual translations) of at least 15 amino acids and considering an E-value cut-off ≤ 1e–5.

We processed the raw Illumina reads for quality control with bbduk (jgi.doe.gov/data-and-tools/bb-tools/) considering the parameters ktrim = r k = 23 mink = 11 hdist = 1. Sanitized reads were then initially assembled with the metaSPAdes pipeline (SPAdes v3.12.0) (Bankevich et al. 2012) considering the k-mer options -k21,33,55,77,99. Then we used SqueezeMeta v1.1.0 (Tamames and Puente-Sánchez 2019) to identify and retrieve individual metagenome-assembled genomes (MAGs) from each DNA sample. The SqueezeMeta pipeline first bins the assembled sequences with both MaxBin (Wu et al. 2015) and MetaBat (Kang et al. 2015), then merges equivalent pairs of binned sequences with DAS Tool (Sieber et al. 2018), and evaluates for quality assembly with CheckM (Parks et al. 2015). Each merged bin represents an individual MAG.

To investigate the phylogenetic affiliation of the recovered good quality MAGs (i.e., completeness ≥ 90% and contamination < 5%), we used the Insert Set of Genomes Into Species Tree 2.1.10 application of KBase (Arkin et al. 2018), which estimates an approximately maximum-likelihood phylogenetic tree with FastTree-2 (Price et al. 2010) considering a reference of 49 Clusters of Orthologs Groups from closely related genomes selected based on sequence identity. The phylogenetic analysis allowed the preliminary identification of good quality cyanobacterial MAGs unambiguously affiliated to the genus Microcystis (details provided in the Results section).

To improve the assembly (scaffolding) of the Microcystis MAG, likely linked to microcystin production (e.g., 27 August 2018 sample) in Washademoak Lake, we first performed independent assemblies of the 27 August read collection with MEGAHIT v1.2.9 (de novo assembly) (Li et al. 2015), MetaCompass v2.0.0 (reference-guided assembly with the complete genome of Microcystis aeruginosa NIES-843 (accession number NC_010296) as reference) (Cepeda et al. 2017), and SPAdes v3.12.0 (with the “trusted contigs” option, using a reference-guided assembly) (Bankevich et al. 2012). Then, we used Metassembler v1.5 (Wences and Schatz 2015) to merge the Microcystis MAGs recovered from our four alternative independent assembling tries, which included the initial metaSPAdes try plus the MEGAHIT, MetaCompass, and SPAdes trusted-contigs assemblies. After merging the four Microcystis MAGs, we used GenomeFinisher v1.4 (Guizelini et al. 2016) and the complete genome of M. aeruginosa NIES-843 as reference for a final scaffolding step.

We used the stand-alone software package of the NCBI Prokaryotic Genome Annotation Pipeline (Tatusova et al. 2016) to annotate the resulting Microcystis scaffolds. Then, to identify homologous genes shared between our Microcystis scaffolds and the complete genomes of the closely related M. aeruginosa strains NIES-843 and NIES-2549 (Kaneko et al. 2007; Yamaguchi et al. 2015), we used the software GET_HOMOLOGUES v3.3.2 (Contreras-Moreira and Vinuesa 2013) with the orthoMCL algorithm (75% minimum coverage in BLAST pairwise alignments and max E-value 1e–05) (Li et al. 2003). The sequence Average Nucleotide Identity (ANI) between our scaffolded Microcystis genome and data from different Microcystis species, including isolates from freshwater bodies in Canada (Pérez-Carrascal et al. 2019), was calculated with OrthoANI v1.2 with default parameters (Lee et al. 2016). Contig mapping to reference sequences was carried out with BWA-MEM (Li 2013) with default parameters.

To investigate the genetic potential of the Microcystis MAG from Washademoak Lake to produce secondary metabolites, we used the DeepBGC v0.1.18 tool (Hannigan et al. 2019) with a 0.9 threshold score to predict the presence of biosynthetic gene clusters (BGCs) and antiSMASH online bacterial version (Blin et al. 2019) to get the most similar known BGC according to MIBiG 2.0 database (Kautsar et al. 2019). We also used DeepBGC to predict the “activity” and “product class” associated to each detected BGC.

Chalky blue-green clumps of cyanobacteria (∼1 mm in diameter) were visible in Washademoak Lake from 30 July (the first day of sampling) through to 23 October 2018, with an intense bloom visible at site B on 30 July. Dolichospermum sp. was identified via microscopy as the dominant species in the 30 July sample. Similar clumps were also observed at Harvey Lake from 13 August to 25 September, though their abundance was very low relative to Washademoak Lake, and this water body was not considered to be blooming during our study period.

Surface water temperature was highest in Washademoak Lake on 30 July (31 °C) and on 13 August in Harvey Lake (28 °C) (Fig. S1). By 10 November, the temperature had decreased to 8 °C in Washademoak, and 6 °C in Harvey, with a brief warm period elevating the temperature/decreasing cooling around 8 October.

More than 70% (34/48) of the surface water samples from Harvey and Washademoak lakes were PCR positive for the presence of the mcyE gene, with positive samples obtained from the beginning of the sampling period in late July through late October. No mcyE-positive samples were detected in November. Interestingly, the mcyE assay was negative for Washademoak site B on 30 July, when an intense cyanobacterial bloom was observed. All mcyE-positive samples were also positive for the Microcystis-specific mcyE primers. Only one sample (Humphrey’s Wharf (site B), Washademoak Lake) was positive for mcyE amplification using the Planktothrix-specific primers, and all evaluated samples were negative for the Dolichospermum-specific mcyE priming sequences (summary of PCR results in Table 1).

PCR screening for the presence of the anaC gene resulted in positives for 18 of the 48 surface water samples from Harvey and Washademoak Lakes (Table 1) collected throughout the summer and fall. No samples tested positive in November.

Samples that tested negative for the mcyE and anaC gene were assayed with a universal primer pair (27F1 and 1492Rl; Table S2) that amplifies the bacterial 16S rRNA gene to determine if the negative results were due to PCR inhibition. All tested samples were positive for 16S rRNA amplicons.

The Illumina reads produced from the last sampling date (10 November 2018) from both lakes were discarded because of low quality. Taxonomic binning of the raw reads with MG-RAST indicates that the bacterial communities in samples collected in Washademoak and Harvey lakes were mainly composed of Proteobacteria, Cyanobacteria, and Actinobacteria, which together account for at least 71% of all recovered major taxonomic groups in any DNA sample (Fig. 2a).

Cyanobacteria dominated in Washademoak Lake on 30 July and 27 August, where their relative abundance was more than 88%. This coincides with peak surface water temperatures (29–31 °C), and an intense cyanobacterial bloom that was visible on 30 July at this location. Bacteroidetes, Firmicutes, Verrucomicrobia, Planctomycetes, Chloroflexi, and Acidobacteria were other bacterial lineages present, but their accumulated relative abundance was not higher than 19% and 26% in Washademoak and Harvey lakes, respectively (Fig. 2a).

The MG-RAST analysis also revealed that members of the Synechococcus, Cyanothece, Microcystis, Nostoc, and Dolichospermum (planktonic Anabaena species have been transferred to the genus Dolichospermum (Wacklin et al. 2009)) genera were the more abundant cyanobacteria in both lakes (65%–79% in Washademoak and 65%–68% in Harvey) (Fig. 2b). Well-known bloom-forming genera (Microcystis, Nostoc, and Dolichospermum) were more prevalent in Washademoak with higher relative abundances (>45% in 4 samples) than in Harvey Lake (not higher than 29% in any sample) (Fig. 2b). In particular, the Washademoak Lake sample collected on 30 July was dominated by Dolichospermum and Nostoc, accounting for circa 70% of the total cyanobacterial abundance; yet, the PCR assays indicate there was no microcystin potential. Microcystis was relatively abundant in the same lake from August through October, peaking on 27 August (>35%) and 8 October (>25%), which coincides with the general and Microcystis-specific PCR assays indicating microcystin production potential from August through October.

Our assembly and binning strategy recovered three good quality MAGs from Washademoak Lake, but it failed to identify assemblies of comparable quality in Harvey Lake samples. The three MAGs assembled from Washademoak Lake DNA samples (from 27 August and 8 October 2018) revealed a good degree of completeness (from 90% to 95%) and low levels of sequence contamination (from 0.5% to 4.5%). The phylogenetic analysis of a concatenated alignment of 49 core genes (Clusters of Orthologous Groups families preselected by KBase) shows that two of three Washademoak Lake MAGs belong to cyanobacteria of the genus Microcystis, whereas the third recovered genome is a β-proteobacteria of the Order Burkholderiales (Fig. 3).

The additional scaffolding steps allowed us to assemble a nearly complete cyanobacterial genome from the 27 August sample (hereafter referred as Microcystis aeruginosa WS75) (Table S3) highly similar (ANI > 97.4%) to Microcystis aeruginosa strains isolated from lakes Champlain and Caron and the reservoir Choinière in Quebec, Canada (Pérez-Carrascal et al. 2019) (Fig. S2), suggesting a close phylogenetic relationship between these lineages from New Brunswick and Quebec (Lee et al. 2016; Jain et al. 2018). The sequence identity of M. aeruginosa WS75 to genomes of other Microcystis species (e.g., M. flos-aquae TF09 and M. wesenbergii TW10) is circa 95% ANI (Fig. S2).

Regardless of differences in size and total number of genes identified, the G+C content and number of CRISPR arrays detected in the M. aeruginosa WS75 assembly (circa 4.9 Mbp) are comparable to those from M. aeruginosa NIES-2549 and M. aeruginosa NIES-843 (Kaneko et al. 2007; Yamaguchi et al. 2015) (Table 2). The repertoire of RNA-coding genes (rRNAs and tRNAs) of M. aeruginosa WS75 is identical to the NIES-2549 strain. In contrast, only 65%–71% of the 4914 predicted genes in M. aeruginosa WS75 are shared with both NIES-843 (3504 shared) and NIES-2549 (3209 shared). The reciprocal ortholog detection also revealed that M. aeruginosa WS75 contains a set of 765 exclusive genes not shared with the other two investigated M. aeruginosa genomes.

The genome-mining survey detected 29 BGCs in M. aeruginosa WS75 accounting for more than 5% of the total assembled genomic sequences of this cyanobacterium (Table S4). According to DeepBGC, a “product activity” label was assigned to 20 of the clusters, whereas only 13 were assigned to a “product class” (Table S4). Although the antiSMASH survey indicates that most of these BGCs have no evident homologous counterparts in the MIBiG repository, 11 of them are similar to known clusters available in this database. These include complete, or near complete, BGCs involved in the biosynthesis of piricyclamides (5.7 Kbp, 91% of the BGC in M. aeruginosa PCC 7005; Fig. S3), the trypsin inhibitor aeruginosin 98- A,B,C (21.5 Kbp, 78% of the BGC in M. aeruginosa NIES-98; Fig. S4), micropeptin K139 (100% of the BGC in M. aeruginosa K-139)/cyanopeptolin (100% of the BGC in Microcystis sp. NIVA-CYA 172/5) (38 Kbp; Fig. S5) and, important to this work, the production of microcystin (60 Kpb; Fig. 4a). The BGC predicted with cytotoxic activity (hereafter WS75 MCY BGC) was very similar in order and orientation to the microcystin BGC of M. aeruginosa PCC 7806 (82% of the BGC) (Tillett et al. 2000) (Fig. 4a). We did not identify evidence of an anatoxin-a BGC in M. aeruginosa WS75.

Our PCR analysis detected the mcyE gene in most DNA samples from both lakes, but besides the WS75 MCY BGC, no other contigs containing the complete MCY cluster were recovered from other samples. However, when contigs assembled with metaSPAdes were mapped against WS75 MCY cluster, fragments (contigs) of the MCY cluster were found in samples from 8 October (WS78), 13 August (WS74), 11 September (WS76), and 25 September (WS77), which matches with the PCR results (Fig. 4b). The contig mapping coverage against the WS75 MCY cluster varies from the almost complete reconstruction of the MCY BGC in samples WS78 and WS74 to the highly fragmented assemblies from WS76 and WS77 (Fig. 4b). Importantly, there is no evidence of the MCY BGC in the Washademoak Lake samples collected on 30 July and 23 October, which is consistent with the negative mcyE PCR results for those same dates (Table 1). No contigs from Harvey Lake mapped to WS75 MCY BGC. We did not find evidence of additional gene clusters associated with production of other cyanotoxins in the metagenomic contigs from both lakes.

The PCR screening indicates that during the 2018 sampling period, the potential for microcystin production, as inferred from the presence of the mcyE gene, was present throughout the summer and early fall of the year in both Washademoak and Harvey lakes. The use of different pairs of genus-specific mcyE primers allowed us to preliminarily identify Microcystis as the main contributor to the latent production of microcystin in these lakes. Despite the high relative abundance of Dolichospermum sequences in both lakes, which has been associated with microcystin production elsewhere (Huisman et al. 2018), no Dolichospermum-specific mcyE sequences were amplified by PCR. This result suggests that toxic genotypes of this taxa were either not present or below the limit of our PCR protocol during the period of study. The PCR survey also indicates that Planktothrix was a rare contributor to toxin production in certain locations. The primers used in our anaC PCR survey were designed to amplify sequences from diverse cyanobacterial genera (e.g., Anabaena, Aphanizomenon, Oscillatoria, Dolichospermum, and Microcoleus); therefore, we were unable to distinguish the particular cyanobacterial taxa associated with anatoxin-a potential. Recent reports have identified anatoxin-a, likely produced by Microcoleus, as a threat for human and animal health in Wolastoq/Saint John River in New Brunswick (McCarron et al. 2019), and our results indicate that this potent neurotoxin is also a concern for lakes in the same region.

The PCR approach permitted the detection of microcystin toxin potential in Harvey Lake, which has had blooms in the past but did not during this study period. This highlights the sensitivity of this approach and its ability to identify lakes at risk for toxic blooms, even under nonblooming conditions. As the frequency and severity of cyanobacterial blooms is anticipated to increase with climate change and cultural eutrophication (Huisman et al. 2018), the sensitivity of PCR provides an early-warning indicator for studying how lakes such as Harvey Lake, which have no or sporadic toxic blooms now, may transition to having more severe or frequent toxic blooms in the future.

The microbial community in both lakes was dominated by Cyanobacteria, Proteobacteria, Actinobacteria, Bacteroidetes, and Verrucomicrobia, which is consistent with previous metagenomic surveys in other lacustrine environments (Eiler and Bertilsson 2004; Pope and Patel 2008; Steffen et al. 2012). Unlike Harvey Lake, typical bloom-forming cyanobacteria (Microcystis, Dolichospermum, and Nostoc) dominated in Washademoak Lake during 2018, and indeed, a surface cyanobacterial bloom was present during the July sampling period. Dominance by bloom-forming cyanobacteria is typically associated with specific environmental conditions (Zohary and Breen 1989; de J. Magalhães et al. 2019). For example, Microcystis presence is commonly correlated with high retention times in the water column, low water transparency, and more stable and prolonged periods of thermal stratification (Soares et al. 2009; Rangel et al. 2016; Guedes et al. 2018), while being more sensitive to mixing and turbulence (Reynolds 2006). Total cyanobacterial biomass is also positively correlated with water temperature and nutrient concentration (i.e., phosphorus) (Beaulieu et al. 2013; Smith 2019). Certainly, the surface water temperatures in Washademoak Lake reached a higher maximum than in Harvey Lake in 2018 (31 °C vs 28 °C) (Fig. S1), and Washademoak Lake is mesotrophic whereas Harvey Lake is oligotrophic (New Brunswick Department of Environment and Local Government 2019).

We assembled a near-complete genome of a toxigenic Microcystis aeruginosa strain inhabiting Washademoak Lake. The difference in genome size and presence of unique genes in M. aeruginosa WS75 relative to other strains of the same species (i.e., NIES-843 and NIES-2549) can be explained, in part, not only by disparities in the number of genes involved in DNA replication, recombination, and repair (Yamaguchi et al. 2018), but also by the high number of taxon-specific genes likely involved in the biosynthesis of toxins and other secondary metabolites, which is an important characteristic of this cyanobacterial lineage (Humbert et al. 2013; Huang and Zimba 2019). The high sequence identity at the genomic level between M. aeruginosa WS75 and M. aeruginosa strains characterized in fresh water bodies in Quebec strongly suggests these isolates are part of the same locally adapted lineage of the globally distributed M. aeruginosa (Pérez-Carrascal et al. 2019).

The recovery of MAGs from environmental samples has opened the possibility of screening for new genes and gene clusters, unlocking the previously underestimated metabolic potential of bacterial communities from lakes and other aquatic systems (McMahon 2015; Cuadrat et al. 2018). Metagenomics has been particularly useful for studying potentially toxic bloom-forming cyanobacteria, avoiding the challenges and limitations of culture-based studies (Alvarenga et al. 2017; Semedo-Aguiar et al. 2018; Pérez-Carrascal et al. 2019). Additionally, genome-scale investigations of bloom-forming and toxin-producing cyanobacteria have been key to appreciating their genetic diversity and capabilities (Pérez-Carrascal et al. 2019), genome plasticity (Meyer et al. 2017), and to improve our capacity to manage freshwater ecosystems during blooming events (Yamaguchi et al. 2020).

Our genome mining analysis revealed a high potential for M. aeruginosa WS75 to produce diverse secondary metabolites. That 5% of the M. aeruginosa WS75 genome involved in the production of secondary metabolites is consistent with previous observations that suggest cyanobacteria dedicate between 5% and 6% of their coding capacity to produce secondary metabolites, with some genomes encoding up to 23 BGCs (Shih et al. 2013; Calteau et al. 2014). Identifying the MYC BGC in the M. aeruginosa WS75 genome is not surprising given that members of this genus are recognized as producers of microcystins (Pineda-Mendoza et al. 2016), but contributes to understanding the genomic footprint of this toxic genus in the Atlantic region of Canada. The MYC BGC of M. aeruginosa WS75 lacks the uma5 and uma6 genes, but the uma coding regions (uma1 to 6) have been reported in genomes of both toxic and nontoxic strains and are likely not involved in the biosynthesis of microcystins (Tillett et al. 2000). Consistent with the PCR survey, we found evidence of fragmented MYC BGCs in all PCR-positive samples from Washademoak Lake.

In addition to the MYC BGC, we also identified diverse gene sets associated with the production of other bioactive secondary metabolites that have not received the same attention as microcystins (Pearson et al. 2019). Cyanobacteria, and species of Microcystis in particular, are known to produce potentially bioactive peptides that are underexplored with respect to their relevance in human health and ecological effects (Janssen 2019). The M. aeruginosa WS75 biosynthetic repertoire includes the ability to produce aeruginosin, which is a potent inhibitor of serine proteases such as thrombin (involved in the human blood coagulation enzymatic cascade) and trypsin (part of the digestive enzymatic battery) (Ersmark et al. 2008; Pearson et al. 2019), with some variants being highly toxic to invertebrates (Scherer et al. 2016). We also identified BCGs involved in the production of cyanopeptolins and micropeptins (CTPs), which can act as inhibitors of the coagulation Factor Xla, human kallikrein, plasmin, thrombin, trypsin, and chymotrypsin. Due to their protease inhibition activity, CTPs cause high mortality rates in freshwater water flea (Daphnia sp.) and fairy shrimp (Thamnocephalus platyurus) (Huang and Zimba 2019).

Secondary metabolites also play important roles in microbial community dynamics and are thought to be produced in response to abiotic and biotic stress, thus providing advantages to the holder over other organisms coexisting in the same environment (Kultschar and Llewellyn 2018). For example, microcystin produced during a Microcystis bloom has an allelopathic effect, inhibiting the growth of Synechococcus elongatus (Hu et al. 2004). In general, bioactive metabolites can modify the abiotic and biotic conditions of aquatic environments to favor dominance of bloom-forming cyanobacteria (Kurmayer et al. 2016).

Our metagenomic analysis failed to recover a MAG from Harvey Lake or genome fragments coding for the anatoxin-a BGC in either lake, despite the anaC gene being detected by PCR analysis. In contrast to PCR detection, shotgun sequencing of environmental DNA is typically biased towards the dominant organisms in each sample, with low probabilities of detecting rare members of the community (Bharagava et al. 2019). Accordingly, abundance of anatoxin-a-producing genotypes in both lakes is expected to be low, and microcystin producing cyanobacteria would not be dominant in Harvey Lake.

PCR-based surveys targeting specific genes involved in toxin biosynthesis are widely used, rapid, and a cost-effective method for qualitative (conventional PCR) and (or) quantitative (qPCR) identification of toxic cyanobacterial genotypes (Dittmann et al. 2017; Rantala-Ylinen et al. 2017). The sensitivity of PCR enables the detection of individual cells even before a cyanobacterial bloom occurs, providing time to initiate appropriate managements actions such as public warnings and to employ more costly chemical analyses for identification and quantification of cyanotoxins (Kurmayer et al. 2017).

The sensitivity of PCR allowed us to identify, but with limited taxonomic precision, the genetic potential for microcystin and anatoxin-a production in two New Brunswick lakes that have experienced recurrent cyanobacterial blooms. This rapid and effective approach requires limited expertise to undertake, making it an accessible tool for other understudied regions where toxin potential is not well-characterized. Findings can then be used to advocate for more comprehensive monitoring programs and employing quantitative approaches, such as qPCR for gene abundance and chemical analyses for toxin concentration. In addition, the PCR results allowed us to direct our metagenomic survey to recover the almost complete genome of the microcystin-producing cyanobacteria inhabiting Washademoak Lake. Our analyses revealed that M. aeruginosa WS75 from Washademoak Lake belongs to a lineage with toxigenic potential widely distributed in other lakes across eastern Canada.

Detailed characterization of cyanobacteria inhabiting Atlantic Canada freshwater systems is important for a region that is relatively understudied and becoming increasingly impacted by cyanobacterial blooms. Although our PCR-based survey was targeted to genes involved in microcystin and anatoxin biosynthesis and the metagenomic approach did not recover BGCs for the production of other cyanotoxins, investigations into the cyanobacterial genetic potential to produce other, but less common, toxic secondary metabolites (e.g., cylindrospermopsins and saxitoxins) should be explored in freshwater systems of Atlantic Canada.

We hope that surveys combining PCR and metagenomic approaches, which also are becoming a relatively low-cost survey option, are regularly employed to explore the toxin and other bioactive secondary metabolite potential of cyanobacteria communities in other understudied regions.