A total of 285 water samples (93 in the Upper Thames River, 47 in the Lower Thames, 70 in Lake St. Clair, 14 in the Detroit River, and 61 in Lake Erie) were collected between January 2016 and October 2019 (sampling frequency detailed in Supplementary Table S1). Samples were taken repeatedly at the same sites over the years and seasons encompassing different systems defined as the Upper and Lower Thames River, Lake St. Clair, Detroit River, and the Western, Central and Eastern basins of Lake Erie (Figure 1). Surface samples were collected in triplicate using a Niskin bottle in a rosette configuration for Lake Erie, a Van Dorn bottle for Lake St. Clair and the Detroit River, and from the shore or from a bridge with a pole holding sterile bottles for the Thames River. All samples were kept in a cooler on ice and filtered within 24 h on board the Canadian Coast Guard Ship (CCGS) Limnos for Lake Erie and in the lab for the other systems. Temperature (Temp), specific conductivity (Cond), pH, and dissolved oxygen (DO) were measured at each site with a SBE 25plus Sealogger (Seabird) on board the CCGS Limnos, and with a Quatro Pro (YSI Inc.) in Lake St. Clair, Detroit and Thames Rivers. Water samples for total phosphorus (TP), total dissolved phosphorus (TDP), total nitrogen (TN), total dissolved nitrogen (TDN), and particulate organic carbon (POC), were submitted to Environment and Climate Change Canada’s National Laboratory for Analytical Testing (NLET, Burlington, Ontario) and prepared following their Standard Operating Procedures for filtration and preservation. Chlorophyll a (chla) concentrations were determined spectrophotometrically after extraction in acetone following the NLET method (NLET, 1997; Table 1).

Estimates of mean daily discharge from 2016 to 2019 were calculated from data available from the Water Survey of Canada¹ for the Thames River and from the U.S. National Water Information System² for the Maumee, Sandusky, St. Clair, and Detroit Rivers.

From 250 to 300 ml of water for DNA extraction was filtered in triplicate through 0.2 μm Polyethersulfone filters (Fisher Scientific) and stored at −80°C until further processing. All genomic procedures were carried out following a standard approach for genomic research conducted by the government of Canada (Edge et al., 2020). DNA was extracted using the Qiagen DNAEasy PowerSoil DNA isolation kits following instructions from the manufacturer. The V4–V5 hyper-variable regions of the 16S rRNA gene were amplified with primers 515F (GTGCCAGCMGCCGCGGTAA) and 926R (CCGYCAATTYMTTTRAGTTT). Libraries were prepared in the Energy, Mining and Environment Biotechnology Research Center of the National Research Council of Canada (Montreal, Quebec). All triplicate samples were sequenced on an Illumina MiSeq platform at the National Research Council (Saskatoon, Saskatchewan). The 16S rRNA gene sequences were analyzed in R Core Team (2018) following the dada2 pipeline (Callahan et al., 2016), using the high-performance computing environment of Shared Services Canada in Dorval, Quebec (Edge et al., 2020). First, non-biological sequences of the primers were removed with cutadapt (Martin, 2012). Then, for each sequencing plate, raw read quality profiles were assessed and the low-quality bases at the end of the read were trimmed with a truncQ score of 11, as suggested for large datasets. Sequences with a maximum expected error (maxEE) greater than 2 were removed and high-quality sequences were merged in an amplicon sequence variant (ASV) table. Chimeras were removed with the ‘removeBimeraDenovo’ command and taxonomy was assigned with the Silva database version 128 released in 2016 (Yilmaz et al., 2014). Additionally, a Barcode of Life Data System (BOLD) reference database that was developed using sequences of cyanobacterial and algal cultures, was used to classify sequences belonging to the Cyanobacteria phylum (Ivanova et al., 2019). Sequences corresponding to Archaea, Chloroplasts, and Eukaryotes were removed from the ASV table and identical sequences that only varied in length were collapsed with the ‘collapseNoMismatch’ command in dada2. The ASV table was rarefied a hundred times at 10,000 reads per sample with the ‘rarefy_even_depth’ command available in the package phyloseq version 1.36.0 (McMurdie and Holmes, 2013), and the average read counts of the 100 tables was used for downstream analyses. This process discarded 383 ASVs out of 68,230 and 50 samples, which when taking into account the number of samples that were collected and sequenced in triplicate, represented a loss of only 9 biological samples.

Dispersion plots for visualizing variation across systems, seasons, and years were generated with the ‘betadisper’ function available in the R package vegan version 2.6.2 (Oksanen et al., 2019). Non metric multidimensional scaling (NMDS) was generated using the command ‘metaMDS’ available in the R package vegan (Oksanen et al., 2019) and based on Bray-Curtis distances on the squared-root ASV table. To examine how environmental variables correlated with the community ordination, we selected the environmental variables measured in all systems (i.e., Cond, pH, temperature, TP, TDP, TN, TDN, POC, and chla) based on the following: a matrix of Spearman correlations among all environmental variables was calculated and those with high correlation (ρ > 0.8 and p < 0.01) with more than one other variable were removed. Hence, TDP and TN were removed for downstream analyses. The remaining environmental variables were projected onto the NMDS ordination using envfit function available in package vegan (Oksanen et al., 2019). Significant differences in microbial community composition as a function of the system sampled or the season were tested with PERMANOVA (Adonis function in package vegan) with 999 permutations.

Mantel correlations were used to assess correlations between geographic distances and structure of the microbial and cyanobacterial communities. Community matrices were first Hellinger transformed with the decostand command in vegan and then detrended (Borcard et al., 2018). Distances between each sampling site coordinate were calculated with the ‘pointDistance’ function available in R package raster version 3.6.3 (Hijmans, 2020) while accounting for the flow direction (from the most upstream to the most downstream site). Additionally, a variation partitioning analysis was employed to disentangle the proportion of variation observed in the microbial community that was due solely to environmental variables or spatial variability. This analysis was performed using the ‘varpart’ function and by transforming geographic coordinates into principal coordinates of neighbor matrices (PCNM) with the ‘pcnm’ function (both are available in the R package vegan; Oksanen et al., 2019). PCNM allows for the detection of any type of spatial patterns and is not restricted to linear ones (Borcard and Legendre, 2002). For this analysis, only environmental variables that were measured in all systems were used (Cond, pH, temperature, TP, TDP, and TDN), and samples containing missing data were removed (Table 1). All environmental variables and PCNMs were selected with the forward.sel function in the adespatial version 0.3.19 package (Dray et al., 2018) prior to being included in the variation partitioning analysis. Consequently, the analysis was performed on 231 samples (85 from the Upper Thames, 40 from the Lower Thames, 38 from Lake St. Clair, 12 from Detroit River, and 56 from Lake Erie).

Stacked bar plots of the microbial and cyanobacterial community composition, merged as a function of the system to improve plot clarity, were drawn with ggplot2 version 3.3.6 (Wickham, 2016). Sparse Correlations for Compositional data (SparCC) at the phylum level between the relative abundance of cyanobacterial and the most abundant bacterial groups, were produced with the sparcc command in the package SpiecEasi version 1.1.2 (Kurtz et al., 2022). Shared taxa across seasons and systems were represented with and upset plot from the UpSetR package version 1.4.0 (Gehlenborg, 2019).

A phylogenetic analysis of the 15 most abundant cyanobacterial ASVs was performed with sequences from the BOLD database (Ivanova et al., 2019) and from GenBank release 239 (Benson et al., 2018). Two Escherichia coli 16S sequences were used as outgroup (accession number HG917881 and HF584705). Sequences were first aligned with MUSCLE algorithm available in the software mega11. The best molecular model was tested with mega and a consensus neighbor-joining tree was then built, based on 1,000 trees, with associated bootstrap values using the best fit model Kimura 2 with gamma parameter.

Generalized additive models (GAMs) were used to assess the relationship between the proportion (from 0 to 1) of cyanobacterial genera relative abundance (arcsine and squared-root transformed) and the following environmental variables: temperature, pH, log-transformed TP, TDP, TN, TDN, and the ratio between TN and TP. GAMs have the advantage of mixing linear model terms with a smoothed and flexible model term and are therefore best suited for non-linear relationships because they can fit non-linear models to data (Hastie and Tibshirani, 1990). GAMs where constructed in R using the gam function in the mgcv version 1.8.36 package (Wood, 2011) and drawn with ggplot2.

The pH level tended to be highest during the summer in Lake St. Clair and Lake Erie (Table 1), and reached a maximum value of 10 in the Western Basin during the summer of 2017. Temperature and DO were more variable across seasons than systems, and temperature was, on average, the highest in Lake Erie surface water during the summer. There was a gradient of decreasing Cond and nutrient concentrations along the water continuum. In particular, the Thames River had higher average TP compared to the rest of the sampled systems, especially in winter for the Upper Thames and in spring for the Lower Thames (Table 1). The highest TP concentration was measured in January 2018 at the most upstream site of the Thames River, reaching 519 μgL⁻¹. TP decreased downstream in Lake St. Clair and averaged below 25 μgL⁻¹ in Lake Erie, even in the Western Basin, which experiences annual cyanobacterial blooms and normally corresponds to a mesotrophic state. However, a brief peak of 79.3 μgL⁻¹ TP was measured during spring 2017 in the Western Erie Basin. In fact, the average TP concentration in the Western Basin was only marginally higher than observed in the two other basins during the sampling period. A similar trend was observed for TDP, where concentrations were on average higher in the Thames River than in the rest of the aquatic continuum. No TN data were collected in the Thames River from 2016 to 2018, so reported values are only based on 2019 sampling. This showed a clear gradient of decreasing TN and TDN concentrations from the Thames River to Lake St. Clair and Lake Erie, where concentrations were on average one order of magnitude lower than in the Thames River (Table 1).

Chla, which can be used as a proxy for algal biomass, was on average higher in the Thames River and Lake St. Clair during the spring, and decreased in summer and fall in all systems. Maximum levels of chla were found in the Upper Thames River during the fall. In Lake Erie, chla levels were on average always higher in the Western Basin compared to the Central and Eastern basin.

The bio-informatics analysis recovered 67,847 ASVs present in 276 samples (85 in the Upper Thames, 55 in the Lower Thames, 62 in Lake St. Clair, 14 in Detroit River, and 60 in Lake Erie). Since there was a strong overlap among the centroids for the 4 years in dispersion plots (Supplementary Figure S1A), all further analyses were performed using the collated data for the 4 years.

Microbial community composition was significantly different among the different systems (PERMANOVA R² = 0.21, p < 0.01) and presented a longitudinal dissimilarity gradient along the hydrological continuum, as pointed out by two different multivariate ordinations (Figure 2A; Supplementary Figure S1B). Samples from the same system were more similar to each other with dissimilarity increasing from the Upper Thames River to Lake Erie’s Western Basin. The Central and Eastern Basins marginally clustered with the rest of the Lake Erie samples (Figure 2A). Cond, pH, Temp, POC, TP, and TDN were significantly correlated with the community spatial structuring (p < 0.01 for all variables). Our analyses indicated that the Thames River microbial community was influenced by higher nutrient concentrations, while those of Lake St. Clair and Lake Erie appeared to be mediated by higher pH, Cond, and temperature. Across the whole aquatic continuum, the microbial community showed distinct patterns across seasons (PERMANOVA R² = 0.15, p < 0.01), although there seemed to be more overlapping across season than system, for example, between summer and fall (Figure 2B; Supplementary Figure S1C).

The microbial and cyanobacterial communities showed the same type of spatial patterns, as inferred by the Mantel correlations, which suggested that both were significantly correlated with spatial extent. Both communities showed the highest positive and significant spatial autocorrelation across samples that were less than 7 km apart (Table 2). Spatial autocorrelation stayed positive and mostly significant for samples within a geographic distance of less than 79 km, suggesting that within this distance, communities were more similar to each other than by chance. Beyond this ~79 km radius, the Mantel r became negative, indicating that beyond this distance, communities were more distinct from each other than by chance. Variance partitioning analysis was used to disentangle the influence of spatial vs. environmental variables. Environmental variables that were measured in every system (TDN, Cond, Temp, TDP, pH, TP) and the 12 PCNM axes representing the spatial variation were selected by the forward selection as significant explaining factors (p < 0.01 for each variable) for the variation partitioning analysis. This analysis explained 35% of the microbial community composition, and attributed 24% (p < 0.01) of the variance in microbial community to the environmental and 23% (p < 0.01) to the spatial variables, with a 12% intersection between the two. This indicates that 12% (p < 0.01) of the variation was solely due to the measured environmental variable and 11% (p < 0.01) solely to the geographic location. Sixty-five percent of the variation was left unexplained.

The following proportions of phyla, genera, and ASVs represent average relative abundances of the 4 years of sampling with the standard deviation in parentheses. The phylum Proteobacteria dominated the Upper Thames River, representing 53% (9%) of the winter microbial community, and then decreasing in relative abundance through the water continuum as observed in each season (Figure 3A). A similar trend was observed for the Bacteriodetes, whereas the Actinobacteria and Cyanobacteria seemed to follow the opposite trend and increased in relative abundance from the Upper Thames River to Lake Erie. Cyanobacteria were further observed to be more abundant during summer and fall than in winter and spring. Here, the different cyanobacterial genera were classified using the BOLD database constructed with species representing a diverse array of cyanobacteria (Ivanova et al., 2019). Taxonomic assignment remained problematic for sequences corresponding to the genera Anabaena, Aphanizomenon, and Dolichospermum, which were then classified as the Anabaena-Aphanizomenon-Dolichospermum complex (AAD complex). During winter, the relative abundance of cyanobacteria was very low and only reached 0.3% of the microbial community in Lake Erie (only one winter sample available) and was mostly composed of Synechococcus and the AAD complex (60 and 40% of the cyanobacterial community, respectively; Figure 3B), while Planktothrix was dominant in the Upper and Lower Thames River (50% and 65% of the cyanobacterial community, respectively). In spring, the relative abundance of cyanobacteria reached a maximum of 1.4% (2%) in the Lower Thames River (Figure 3A), and was mostly composed of a mix between the genera Planktothrix and Synechococcus along the entire water continuum (Figure 3B). In summer, cyanobacterial relative abundance increased in the whole water continuum, reaching up to 17% (11%) in the Western Basin of Lake Erie (Figure 3A), and was mostly composed of the genera Synechoccocus and Microcystis, which increased in relative abundance from the Upper Thames River to the Western Basin of Lake Erie. At this season, Synechoccocus represented 57% of the cyanobacterial community in Lake St. Clair, 40% in the Detroit River, 50% in the Western Basin of Lake Erie, and 68% in Lake Erie Central and Eastern Basin, while Microcystis was the second most abundant genus with 25% of the cyanobacterial community in Lake St. Clair, 39% in the Detroit River, and 23% in Lake Erie Western Basin (Figure 3B). Microcystis was poorly represented in the central and eastern basin of Lake Erie, which were mostly dominated by Synechococcus. In contrast, Planktothrix was the most abundant cyanobacterial genus in the Upper Thames River, reaching 55% of the cyanobacterial community. During the fall, the cyanobacterial community stayed relatively abundant reaching 18% (14%) of the Lake Erie’s Western Basin, with Microcystis as the dominant cyanobacterial genus (ca. 45% of the cyanobacterial community), while Pseudanabaena dominated in the Central and Eastern Basins. Planktothrix was dominant in Lake St. Clair composing 32% of the cyanobacterial community.

The genetic connectivity was assessed amongst part of the system that were closely connected from upstream and downstream, and because of their poor connection with the rest of the water continuum, notably highlighted by the NMDS analysis (Figure 2A), Lake Erie’s Central and Eastern Basins were removed from subsequent analyses.

The highest number of ASVs was found in the Thames River, followed by Lake St. Clair (including the Detroit River) and Lake Erie (Western Basin; Figure 4). In each system, the number of ASVs was the highest during the summer, but overall, few were shared among systems and seasons. The data showed a ‘core microbiome’ of 287 ASVs shared across all systems and seasons, which represented less than 1% of all the 67,847 ASVs identified. The relative abundance of this ‘core microbiome’ varied from 10% of the community during the winter in the Upper Thames, to 77% during the spring in the Detroit River, and was composed of the same dominant phyla identified in Figure 3. Amongst those, the proportion of the cyanobacterial phyla in this core was from 0.04% of the microbial community during the winter in both Upper and Lower Thames, to 13% of the community during the fall in the Lower Thames. During spring and summer, the greatest number of ASVs were unique to the Thames River system and not found elsewhere along the water continuum. A smaller number of ASVs seemed to be specific to downstream Lake Erie Western Basin, and this number was the highest during spring (1823). The summer community in Lake Erie’s Western Basin had more ASVs in common with the summer community of Lake St. Clair and the Thames River or the fall community of Lake Erie than with the Western Basin of Lake Erie spring community.

The composition of the microbial community in Lake Erie’s Western Basin was highly connected to the rest of the watershed in terms of a high proportions of the reads (Figure 5A). A smaller proportion of ASVs (Figure 5B) in Lake Erie’s Western Basin were also detected in the upstream systems. The same pattern was observed for the cyanobacterial community (Figures 5C,D). In all seasons, the Thames River contributed the highest proportion of reads in the Western Basin of Lake Erie for both the whole microbial and the cyanobacterial community followed by Lake St. Clair and Lake Erie, while the contribution of the Detroit River was minimal. On the other hand, the greatest proportion of the bacterial and cyanobacterial ASV, at least in spring and fall, were first and only identified in the Lake Erie’s Western Basin.

Genetic connectivity within the watershed was further explored for the 15 most abundant cyanobacterial ASVs using a phylogenetic analysis. Those ASVs clustered with the genera Microcystis, Gloeocapsa, Planktothrix, Aphanizomenon, Anathece and Synechococcus (Figure 6A). Three ASVs (# 2, 5, and 10), clustering with representative of the genus Microcystis, seemed to dominate in different parts of the aquatic continuum and at different times (Figure 6B), and could generally be detected along the water continuum. During summer, ASV 2 started to become dominant in Lake St. Clair and reached 3.4% (3.8%) of the microbial community, followed by 3.7% (3.2%) in the Detroit River and 2.2% (2.5%) in Lake Erie’s Western Basin, it only accounted for less than 0.01% of relative abundance in the Upper Thames River. On the other hand, the other two ASVs clustering with Microcystis (ASV 5 and 10) seemed to become more abundant in Lake Erie during the fall, reaching 5% (8.3%) and 1.8% (1.4%) of the microbial community, respectively, but stayed in lower relative abundance in Lake St. Clair. ASV 12 clustered with the culture representative³ identified as Gloeocapsa, which peaked in relative abundance in the Western Basin during summer but represented only 0.5% (0.7%) of the microbial community. The next cluster in the phylogenetic tree contained two ASVs of the genus Planktothrix (ASV 1 and 3) that both dominated in the Thames River during summer and fall, especially ASV 1 which accounted for 9.9% (10.3%) of the microbial community in the Lower Thames River during fall. ASV 4 clustered closely with the sequence of the culture representative Aphanizomenon flos-aquae that was more abundant in the Thames River, reaching 4.4% (7.6%) of relative abundance in the Lower Thames during the fall. This ASV was also present in Lake St. Clair and Erie, although in lower relative abundance. Finally, the last cluster of the phylogenetic tree gathered sequences of 6 different ASVs clustering with the genus Synechococcus and two with the genus Anathece, both genera belonging to the Order Synechococcales. Those ASVs started to become more abundant in Lake St. Clair, followed by the Detroit River, and Lake Erie during summer and fall, and some were not detectable at all in the Upper Thames River (ASVs 7, 15 and 6).

Because of the non-linear nature of the relationships, Generalized Additive Models (GAMs) were used to test those correlations, and, although all the relationships presented here were significant (p < 0.01), the explanatory power of the R² remained relatively low (from 0.02 to 0.47), and no clear pattern was observed between the set of selected environmental variables and the relative abundance of dominant Cyanobacteria (Microcystis, Synechococcus, and Planktothrix). Compared to Planktothrix, both Microcystis and Synechoccocus genera seemed to peak at lower nutrient (N and P) concentrations which corresponded to concentrations found in Lake Erie in comparison to the Thames River (Supplementary Figure S2). It was noteworthy that across all sites, the relative abundance of both Microcystis and Synechoccocus tended to decrease as nutrient concentrations increased. Although significant and positively correlated, temperature had an extremely low R² for Microcystis (R² = 0.02) compared to Synechococcus (R² = 0.3). No clear relationship was observed between cyanobacterial relative abundance and pH. Finally, the relative abundance of Microcystis seemed to peak at a TN/TP mass ratio ≥ 30, which generally corresponds to the multi-annual average TN/TP mass ratio (35) of Lake Erie’s Western Basin, while the other systems had higher average TN/TP mass ratio values (Central and Eastern Basin of Lake Erie = 47, Lake St. Clair = 40, Detroit River = 41, Lower Thames = 48, and Upper Thames = 80). However, no clear trend was observed between Synechoccocus or Planktothrix relative abundances and nutrient mass ratios.

We observed longitudinal and seasonal changes in the aquatic microbiome within the Lake Erie watershed from the Upper Thames River to Lake St. Clair and Lake Erie. We observed a low relative abundance of cyanobacteria in winter and spring, followed by an increase in their relative abundance during summer and fall which is similar to observations elsewhere when cyanobacterial blooms tend to occur due to higher temperature optima and preference for strong stratification (Paerl and Huisman, 2009). Nevertheless, PERMANOVA analyses suggested that the spatial pattern was more influential than the seasonal pattern on the microbial community structure. While part of the same water continuum, the aquatic systems sampled here represents a gradient of environmental conditions and trophic status, with the Thames River draining nutrients from its predominantly agricultural watershed, which are gradually diluted downstream in Lake St. Clair and Erie. Lake Erie itself represents a longitudinal gradient in trophic status where its Western Basin is mesoeutrophic while the Central and Eastern Basins trend to oligomesotrophic (Dove and Chapra, 2015). The spatial directional pattern observed for the microbial community composition has already been observed in previous meta-community studies on boreal inland waters (Ruiz-Gonzalez et al., 2015), Arctic systems (Crump et al., 2012; Comte et al., 2018), and across the Great Lakes (Rozmarynowycz et al., 2019; Paver et al., 2020). Those studies, together with our findings on the most important tributary on the Canadian side, highlighted the importance of hydrological processes in shaping the microbial community structure in aquatic continuums.

The mantel correlations showed that, for both microbial and cyanobacterial communities, the highest spatial autocorrelation occurred within a radius of 7 km (Table 2). At a wider spatial scale, within a 79 km radius, communities remained more similar than by chance, indicating a certain degree of selection of the microbial community within similar types of system. This radius matches the limit of each system sampled (Upper and Lower Thames, Lake St. Clair, Detroit River and Lake Erie’s Western Basin), where the greatest distance between two samples of a similar system was 54 km between the two most distant Lower Thames samples. Beyond 79 km, a distance-decay pattern started to appear, which indicate dispersal limitation and system filtering as observed in other aquatic ecosystems (Logares et al., 2018). Even within Lake Erie, where despite having the shortest residence time of all the Great Lakes (ca 2.7 years; Quinn, 1992), there was a distinction in the microbial and cyanobacterial community structure between the Western Basin and the rest of the lake, probably driven by the difference in trophic status and hydrology. Collectively, the results indicated that both microbial and cyanobacterial community structures were significantly correlated with the spatial distances and that, beyond a threshold of 79 km which also corresponded to samples within the same system, communities switched from being more similar to more dissimilar than by chance.

Environmental variables including conductivity, temperature, pH, particulate organic matter and nutrients (TP and TDN) have been identified as variables which significantly influenced microbial community structure. Water chemistry parameters often shape aquatic system due to their direct impact on microbial metabolism. Nutrients, and by extension, trophic status, have a prime effect on microbial community structure (Lindström, 2000) and diversity, with higher richness associated with more eutrophic lakes (Kiersztyn et al., 2019), yet richness tends to decrease with N fertilization (Schulhof et al., 2020). As an integrator of landscape properties, pH is often identified as one of the main factors influencing microbial community structure (Lindström et al., 2005; Fierer and Jackson, 2006; Logue and Lindström, 2008; Niño-García et al., 2016). However, variance partitioning analysis of our data indicated that environmental variables per se only explained 12% of the variation observed within the microbial community, and this proportion reached 35% when accounting for both environmental and spatial patterns. This relatively low explanatory power is not uncommon when investigating factors that influence aquatic microbiomes (Marmen et al., 2020), and highlights the importance of examining other factors that can play key roles in microbial community assembly, for example, by including a more complete characterization of the hydrological properties of the systems, including the lag time (Lindström et al., 2005; Marmen et al., 2020), and taking into account the biotic interactions with co-occurring microbes or with higher and lower trophic levels (Louca et al., 2016).

The rest of the aquatic microbiome community interacts extensively with cyanobacteria, notably through mutualistic interactions (Eiler and Bertilsson, 2004; Woodhouse et al., 2016; Berg et al., 2018), toxin degradation (Mou et al., 2013), or predation (Rashidan and Bird, 2001); and bloom events may create a disturbance to this microbiome (Berry et al., 2017). The typical dominant freshwater bacterial phyla, i.e., Proteobacteria, Bacteroidetes, Actinobacteria, and Verrucomicrobia, are often found associated with cyanobacterial blooms (Eiler and Bertilsson, 2004; Bagatini et al., 2014; Tromas et al., 2017; Li et al., 2020), while Proteobacteria have been associated with two strains of cultured cyanobacteria (Microcystis and Cylindrospermopsis; Bagatini et al., 2014). However, in our case, Proteobacteria’s relative abundance was weakly and negatively correlated with cyanobacteria’s (SparCC ρ = −0.05, p < 0.01), and Proteobacteria were found in higher relative abundance in the Thames River compared to the downstream systems. Zhu et al. (2019) also noticed a significant negative correlation between the two phyla relative abundances, suggesting a potential competitive interaction. Negative correlations have been observed between cyanobacteria and Actinobacteria’s relative abundances (Ghai et al., 2014; Matson et al., 2020), with Actinobacteria occurrence being associated with less eutrophic systems (Haukka et al., 2006), although actinobacterial clades were abundant in eutrophic Lake Taihu (Tang et al., 2017). Yet, some taxa affiliated with the Actinobacteria phylum co-occur with cyanobacterial blooms (Woodhouse et al., 2016) or having the relative abundance of certain clade strongly correlated, either positively or negatively, to bloom occurrence (Berry et al., 2017). In our case, Actinobacteria’s relative abundance was positively correlated with cyanobacteria’s (SparCC ρ = 0.4, p < 0.01) and were found in higher relative abundances in the downstream Lakes St. Clair and Erie. Some taxa belonging to the Verrucomicrobia have been shown to display a higher diversity in the presence of cyanobacteria (Parveen et al., 2013) and are known to degrade algal polysaccharides and organic matter (Bagatini et al., 2014; Woodhouse et al., 2016). However, in our dataset, there was no significant relationship between relative abundances of Verrucomicrobia and Cyanobacteria. Bacteroidetes have also been shown to increase in absolute and relative abundance during or after a cyanobacterial bloom, especially the Sphingobacteria and Flavobacteria classes (Newton et al., 2011; Bagatini et al., 2014). However, we observed a weak but negative correlation between Bacteroidetes and Cyanobacteria’s relative abundances (SparCC ρ = −0.07, p < 0.01), which suggests that other factors not investigated in this study are important. In contrast, the relative abundance of Planctomycetes and Cyanobacteria were positively correlated (SparCC ρ = 0.5, p < 0.01), consistent with other reported associations between Planctomycetes and phytoplankton blooms (Eiler and Bertilsson, 2004). Because of their involvement in nutrient cycling, the aquatic microbiome is likely to play a crucial role in bloom occurrence and duration, and it has already been highlighted that the rest of the aquatic microbiome could be a better predictor of cHAB events than environmental factors (Tromas et al., 2017). More studies are needed to explore those interactions in further detail, notably by including eukaryotic microbes, and unveiling the mechanisms underlying these interactions.

Microscopy and molecular tools have been found to identify the same predominant taxa of cyanobacteria, with the exception of pico-cyanobacteria which are underestimated by microscopy (MacKeigan et al., 2022). Here, molecular tools allowed detection of Synechococcus as one of the dominant genera of cyanobacteria in Lake Erie’s watershed, along with Microcystis, Planktothrix, and Aphanizomenon, which were already identified in Lake Erie using microscopy (Millie et al., 2009; Chaffin et al., 2011; Allinger and Reavie, 2013; Davis et al., 2015). Other studies based on meta-barcoding results also identified Microcystis, Planktothrix, Synechococcus, Aphanizomenon, and Dolichospermum as dominant genera of cyanobacteria in river mouths (Sandusky and Maumee) and in the Western and Central Basins of Lake Erie (Berry et al., 2017; Salk et al., 2018; Jankowiak et al., 2019; Rozmarynowycz et al., 2019; Matson et al., 2020). Synechoccocus and Microcystis were found to co-occur as dominant genera in Lake St. Clair, Detroit River, and Lake Erie (Figure 3B), and this co-occurrence was already observed in Lake Erie (Berry et al., 2017; Rozmarynowycz et al., 2019) and in Lake Taihu (Ye et al., 2011). Synechococcus also initiated the 2014 cyanobacterial bloom in the Western Basin of Lake Erie and stayed abundant until the end of the bloom (Berry et al., 2017). In the Thames River, sequences corresponding to the genera Planktothrix and Aphanizomenon were found in high relative abundances (Figure 6B). To date, very few blooms have been reported in the Thames River, although there is increasing evidence that the river experiences recurrent blooms (Crevecoeur, Molina and Watson, pers. communication, Supplementary Figure S3). A bloom reported in the lower portion of the river around Chatham and identified by microscopy was dominated by Aphanizomenon flos aquae and Planktothrix agardhii (McKay et al., 2020), which is consistent at the genus level with what we found in our genomic data from 2016 to 2019. However, our data also suggest the presence of ASVs clustering with Planktothrix rubescens and P. suspensa (Figure 6A). As a comparison, Plantkothrix has also been reported as abundant on the US side of Lake Erie, specifically in the Sandusky Bay (Davis et al., 2015; Salk et al., 2018; Jankowiak et al., 2019). These observations suggest that Microcystis and Synechococcus might be well adapted to lake conditions and Planktothrix to more turbid river-type conditions. Our data also highlight the existence of complex assemblages of different strains of cyanobacteria that successively dominate the water continuum through space and time. Further investigation is needed to determine if these different strains and ASVs respond in the same way to environment triggers for toxin production.

While the connectivity between Lake Erie and its tributaries has been clearly identified in terms of P loading (Maccoux et al., 2016), the connectivity of the aquatic microbiome and the contribution of riverine Bacteria to the lake microbiome remain uncertain. The Thames River holds a high proportion of microbial taxa that were only identified in the upstream river but not in any of the downstream systems (Figure 4). Lake St. Clair appeared to act as a filter and is also influenced by other sources not taken into account here, such as the upstream Huron corridor. In summer and fall, there were more taxa shared between Lake Erie and the rest of the watershed than between Lake Erie in the spring. This highlights the seasonal succession and change in aquatic microbial community.

Within a watershed, the most commonly observed patterns of microbial community assembly are mass effect, i.e., massive advection of numerically dominant microbes by passive transport, and species sorting (Lindström and Bergström, 2004; Zhao et al., 2021). In our case, the majority of microbial reads found in the Western Basin of Lake Erie were likewise identified from the Thames River (Figure 5). This observation allows to unveil the impact of dispersal within the watershed in general, but this does not rule out other sources, notably amongst other more proximal tributaries, or with higher discharge, and/or the sediment, as ASVs can be shared across multiple watersheds (Urycki et al., 2022). Other studies have classified reads and ASVs as a function of the system in which they were first detected (Crump et al., 2012; Ruiz-Gonzalez et al., 2015; Matson et al., 2020), and our data clearly indicated that mass effect played an important role in the downstream microbial community assemblage. This phenomenon is prevalent when dispersion is high due to elevated levels of connectivity between the upstream river community and the downstream lake (Adams et al., 2014; Ruiz-Gonzalez et al., 2015). However, in lake systems, where water flow generally decreases and residence time increases, species sorting occurs (Zhao et al., 2021). Indeed, our data showed that certain cyanobacterial ASVs, mainly ASVs 2, 5 and 10, clustering with Microcystis isolates, were found in low relative abundance in the Upper Thames River, but became dominant in the Western Basin of Lake Erie (Figure 6), suggesting the lake conditions were selecting these cyanobacterial ASVs to become dominant. This genetic connectivity between upstream and downstream has already been highlighted for microcystin producers such as Microcystis, Planktothrix, and Anabaena in Lake Erie, Lake St. Clair, and Lake Ontario based on the mcyA gene responsible for microcystin production (Davis et al., 2014). Our data suggest that this connectivity extends beyond the Great Lakes to the upstream tributaries, and also across different seasons. As to whether Lake Erie Microcystis-dominated harmful blooms are seeded from the river, based on our data it seems unlikely that the Thames River is the sole source for those Microcystis ASVs given the extremely low relative abundance in which they were recovered. Additional sources likely contribute to the rapid development of summer and fall blooms in the Western Basin of Lake Erie, for example, recruitment from the sediment (Kitchens et al., 2018), or other tributaries (Maumee, Sandusky, and St. Clair Rivers), which still need to be investigated as both evidences supporting (Bridgeman et al., 2012) or refuting (Kutovaya et al., 2012) riverine seeding from the Maumee River have been documented.

In our study, no clear pattern was observed between cyanobacteria relative abundance and nutrient concentrations, although our results suggested that Planktothrix relative abundance peaked at higher N and P concentrations than Microcystis and Synechococcus (Supplementary Figure S2). Planktothrix also dominated in an experiment with nutrients (N, P) addition while, under the same conditions, Microcystis decreased in relative abundance (Jankowiak et al., 2019). A similar observation was made by Harke et al. (2016) based on shipboard incubations of samples collected from a bloom site on Lake Erie’s Western Basin. These authors found that Microcystis abundance significantly decreased with P enrichment, while Planktothrix and Dolichospermum (formerly Anabaena) dominated samples collected from the high P environment near the mouth of the Maumee River. On the other hand, Planktothrix dominance has also been observed in the Sandusky river plume, where it’s capacity to accumulate N through cyanophycin compounds provided a competitive advantage under N limiting conditions (Hampel et al., 2019).

Microcystis seems to be more efficient at scavenging P, giving it an advantage to dominate at lower nutrient concentrations (Gobler et al., 2016; Harke et al., 2016), but it is also more adapted to stable water conditions because of its ability to regulate its buoyancy (Dokulil and Teubner, 2000). In contrast, Synechococcus is generally dominant in oligotrophic lakes (Ruber et al., 2016). It has been argued that the ability of some cyanobacteria to fix N provides a competitive edge in eutrophic systems with low N:P ratios (Smith, 1983), however this tenet has been debated (Downing et al., 2001). Our data showed generally negative correlations between N concentrations and the relative abundance of some non-diazotroph genera like Microcystis and Synechococcus, which showed higher relative abundances at lower N:P ratios (Supplementary Figure S2). These genera could be complemented by the presence of diazotrophic cyanobacteria and Proteobacteria (Davis et al., 2015; Li et al., 2020). Therefore, the presence, for instance, of N2-fixing Aphanizomenon and Cylindrospermopsis along the water continuum, could play a role in promoting the growth of Microcystis and Synechococcus. More investigations are needed to address those questions but molecular tools have the potential to provide more insight into the role of the entire microbiome, often overlooked in early studies.

Increasing temperature and longer warmer periods likely to happen due to climate change will also have an effect on cHABS, because a more stable thermocline and reduction in epilimnion mixing can lead to better light competition outcomes for cyanobacteria due to their ability to adjust their buoyancy. For instance, when nutrient concentrations are high enough to sustain bloom formation, climatic variables may become crucial for predicting bloom onset and duration, as was observed by Zhang et al. (2012), who found that increasing temperature, hours of sunshine, and reduction of wind speed in Lake Taihu lead to earlier onset and longer bloom duration. However, Taihu is a subtropical lake, and because Lake Erie has a temperate climate, the variables influencing its microbiome might be very different. Since our study was conducted across different seasons, we were likely to identify an effect of temperature. The relationship was positive and significant for Microcystis, Synechoccocus, and Plantkothrix, but the R² of the GAMs were very low for Microcystis and Planktothrix (0.02 and 0.11, respectively); hence it is difficult to draw conclusions regarding the link between cyanobacteria relative abundance and temperature. However, climate projections predict that Lake Erie will have longer and more stable periods of stratification, which are likely to lead to increasing frequency and worsening impact of cHABs (Paerl and Huisman, 2008).