A microbial perspective on the local influence of Arctic rivers and estuaries on Hudson Bay (Canada)

Microbial assemblages are selected largely by environmental conditions, with a clear divide between typical freshwater and marine communities and a sensitivity of the transition from freshwater species to marine species to the sharpness of the salinity gradient (Herlemann et al., 2011). Within estuaries, salinity gradients are determined by freshwater flow rates and mixing regimes (Pritchard, 1989), and the diversity and spatial distribution of bacteria has been linked to salinity (Bouvier and del Giorgio, 2002; Crump et al., 2004, Kieft et al., 2020). Arctic and sub-Arctic systems are also affected by differences in source waters, including river catchments, which influence bacterial communities downstream. For example, dissolved organic carbon (DOC) from glacier versus nonglacier source waters modified the microbial communities in rivers flowing into Disko Bay, Greenland (Hauptmann et al., 2016). This finding is not surprising, as DOC characteristics would select for microbes with the metabolic capacity to degrade pools of source-specific molecules (Jørgensen et al., 2014; Sipler et al., 2017; Colatriano et al., 2018). In systems with similar DOC inputs, differences in communities would be determined by other biotic and abiotic factors. Ghosh and Bhadury (2019) examined estuarine microbial communities along a latitudinal gradient and suggested that temperature was a significant factor in determining microbial community structure. However, high latitude regions were not included in their study.

High latitude regions are of special interest because, in addition to seasonal temperature changes, sea-ice melt in spring contributes to the freshwater inventories of sub-Arctic and Arctic seas and their estuaries, influencing microbial communities (Campbell et al., 2018) with affects on the entire microbial food web (Kellogg et al., 2019). Added to seasonal sea-ice melt, which was relatively constant in the past, the Arctic and sub-Arctic seas are generally freshening with the loss of multiyear sea ice, melting permafrost, and increasing river flow as well as the increasing contributions of glacier melt water (IPCC, 2019). The increasing river flows and permafrost melt can also increase the amount of particulate material entering estuaries, which could alter microbial communities as particles are associated with characteristic bacterial taxa in both marine and freshwaters (Smith et al., 2013; Mohit et al., 2014). For Hudson Bay, inorganic nutrient and organic carbon concentrations are all greater in the rivers contributing to the bay compared to its seawater (Kuzyk et al., 2008), so that the marine waters of the bay are more oligotrophic than their respective estuaries for most of the year.

As a first step toward understanding freshwater-to-marine transitions by microbial communities in an Arctic context, we investigated 2 river-to-sea systems as part of a larger project (BaySys) on the impacts of climate change and river regulation in Hudson Bay (Lukovich et al., 2021). Hudson Bay is an inland sea, classified climatically as Arctic to sub-Arctic, that receives 522 to 760 km³ of freshwater annually from surrounding watersheds (Déry et al., 2011). For this study, we focused on 2 river systems along western Hudson Bay: The Nelson River, which is regulated to have a year-round high flow rate for stable electricity generation; and the much less regulated lower flow Churchill River. The Nelson and adjacent Hayes rivers both have headwaters near Lake Winnipeg and, along with the Churchill River, are part of a complex river network flowing through low relief prairie terrain (Stadnyk et al., 2020).

The Churchill River watershed has an area of approximately 284,000 km² (Clair et al., 2013). Its downstream region forms a semi-enclosed estuary 3 km wide and up to 13 km offshore (Kuzyk et al., 2008). In 1976, hydraulic control systems diverted 75% of the natural water flow rate from the Churchill River to supply a series of hydroelectric dams located on the Nelson River (Newbury et al., 1984). This water diversion from the Churchill River and the dam system along the Nelson River impacts sediment transport in these rivers (Duboc et al., 2017). Before this hydroelectric development, the Churchill River water flow rate was approximately 1,350 m³ s^–1; after the diversion, it decreased to approximately 600 m³ s^–1 (Duboc et al., 2017). Water from the Churchill diversion to the Nelson River increased the average discharge entering the Nelson River Delta from 3,570 m³ s^–1 to 4,380 m³ s^–1 (Prinsenberg, 1980). These anthropogenic disturbances alter the natural seasonal variations of the water flow rate in the Nelson River, flattening the average annual hydrographic curve (Déry et al., 2011). The seasonal flow pattern was not altered for the Churchill River, where in summer modest flow of warmer and fresher water from the Churchill River results in marked vertical stratification in the Churchill Estuary. This stratification modifies the distribution of nutrients and their availability to phytoplankton in the surface waters.

The Nelson River watershed has an area of 1,124,592 km². Its downstream region forms an open estuary (Clair et al., 2013). Compared to the Churchill River, the Nelson River has a higher concentration of total suspended solids (TSS; Godin et al., 2017) and can be turbid to the point of reducing light penetration in the water column (Bilotta and Brazier, 2008). Adjacent to the Nelson River, and flowing into the Nelson Estuary, the Hayes River just to the southeast has a watershed area of approximately 103,000 km² (Godin et al., 2017). The Hayes River has a lower TSS concentration than the Nelson River.

Our overarching goal was to describe the bacterial diversity in these 2 river-to-sea systems in the context of environmental conditions at a critical time of year following seasonal ice melt along western Hudson Bay. Given the major differences in the river flow rates and estuary type, we hypothesized that microbial community composition in the corresponding estuarine and coastal areas would diverge, compared to the riverine communities, which are derived from prairie terrain networks. In addition, the moderate salinities of Arctic and sub-Arctic seas could potentially favor marine microbes adapted to lower salinities in fresh-to-saltwater transition zones irrespective of river flows and tidal mixing. Although many studies to date have mostly relied on broad taxonomic categories to characterize microbial communities, we exploited results from high throughput amplicon sequencing to address the questions of whether microbes in the transition zones were a combination of marine and riverine taxa that enter the estuary or a result of autochthonous communities that developed within the estuary. Specifically, we used high throughput amplicon sequencing of the V4 region of the 16S rRNA gene (rDNA) to identify and follow taxa at the level of operational taxonomic units (OTUs), within specific habitats and size fractions from freshwater rivers, across 2 contrasting transition zones, and into adjacent coastal waters. The V4 region of 16S rRNA (rRNA) was also sequenced in the majority of samples and used to distinguish OTUs from purely detrital or aerial inputs and verify that the taxa found in a given habitat were potentially active there.

Samples were collected in July 2018 along 2 transects starting in the Nelson and Churchill rivers (Figure 1A). One additional sample was collected from the Hayes River, an unregulated river with its outflow near the mouth of the Nelson River. Sampling stations were located along a salinity and temperature gradient in each region, terminating in fully marine stations in Hudson Bay, as detailed in Jacquemot et al. (2021; Figure 1B). The transects included river surface samples collected from the river edge accessed by helicopter, with the Nelson and Hayes rivers (stations NE-A and HA-A, respectively) sampled on June 18, 2018, and the Churchill River (CH-A) sampled on June 28, 2018. All other Nelson area samples were collected from June 29 to July 1, 2018 (Figure 1C). Due to logistical constraints only one intermediate salinity site was sampled along the Churchill transect (Figure 1D), which was chosen based on salinity characteristics. The Churchill stations were sampled July 4–7, 2018 (Table S1).

Aside from the helicopter samples, shallow water upper estuarine samples were collected from either a zodiac or barge deployed from the research icebreaker CCGS Amundsen. All marine and outermost estuarine samples were collected using the ship rosette system, which was equipped with 12-L Niskin-type bottles and a conductivity, temperature, and depth (CTD) profiler (Sea-Bird SBE-911 CTD, Sea-Bird Scientific, Bellevue, WA). The rosette was also equipped with sensors for relative nitrate (In-Situ Ultraviolet Spectrometer, Satlantic, Halifax, NS, Canada), dissolved oxygen (Sea-Bird SBE-43), chlorophyll fluorescence (Seapoint, Exeter, NH), fluorescent colored dissolved organic matter (fCDOM; Wetlabs ECO, Philomath, OR), photosynthetically active radiation (PAR; Biospherical Instruments QCP3200, San Diego, CA), and transmissometry (WETlabs C-Star, Philomath, OR). For shallow waters, conductivity and temperature were obtained using a Sea-Bird 19+ CTD probe deployed from the zodiac or barge, and discrete water samples were collected by bucket or submersible pump (for the deeper samples) directly into clean carboys that were returned to the ship for processing. Ship samples were collected directly from 12-L Niskin-type bottles that were closed on the upward cast after identifying the pycnocline from the CTD profile (Jacquemot et al., 2021). For all but the upper river stations, water was collected from 2 different depths: surface (approximately 0 m) and deeper (between 4 and 11 m) just below the thermocline (bt), which was near the bottom for the shallow stations. For all stations and depths, water samples were collected for nutrients, chlorophyll a (Chl a), flow cytometry (FCM), and nucleic acids, as detailed in Jacquemot et al. (2021).

The bacterial concentration in each sample was estimated by FCM based on Marie et al. (1999) with modifications described below. All FCM samples were fixed by adding 90 µL of 25% glutaraldehyde to 1.8 mL of seawater. Preserved samples were refrigerated at 4°C for 30 min, deep-frozen in liquid nitrogen and stored at –80°C until laboratory analysis. For the bacterial and archaeal counts, referred to as bacterial counts for convenience, 200 µL of each sample were stained with 0.5 µL of Sybr™ Green I (Invitrogen™, ThemoFisher Scientific, Waltham, MA) from a 1000× stock solution. Samples were counted using a BD Accuri™ C6 flow cytometer (BD Biosciences, Franklin Lakes, NJ) equipped with a CSampler, using 96-well plates. Data acquisition was performed at slow flow rate (14 µL min^–1) for 5 min with 3 wash and 3 agitation cycles between each sample. To correct flow speed variation linked to salinity differences and standardize the counts, we added a known quantity of beads (BD Trucount^TM) to standards matching the salinity of each sample, with standards included for each separate run. The salinity at the NE-NS station was not available to correct the flow speed variation for its corresponding FCM samples, so we estimated the salinity based on the nearby NE-45 station (station locations in Figure 1C). The FCM samples from the NE-45 station were lost. For photosynthetic cyanobacteria and other phytoplankton, unstained samples were used. Fluorescence and scatter properties were measured with a 14.7 mW 640 nm Diode Red Laser and 20 mW 488 nm Solid State Blue Laser. For phytoplankton, data acquisition was performed at a fast flow rate (66 µL min^–1) for 10 min with 3 wash and 3 agitation cycles between each sample. Data were processed with BD CSampler software. The setup was only accurate for phycocyanin (PC)-rich cyanobacteria, which we take as a proxy for the cyanobacteria because the FCM system does not include a green laser and is not optimal for distinguishing PC-rich from phycoerythrin (PE)-rich cyanobacteria. We also note that PE-rich cyanobacteria could be confounded with cryptophytes containing PE pigments.

Nucleic acids were collected by filtering 6 L of water first through a 50-µm nylon mesh (prefilter) and then filtered sequentially by using a peristaltic pump. The particulate material was collected on a 47-mm diameter Millipore polycarbonate filter, with a pore size of 3 μm, then a 0.22-µm porosity Sterivex™ cartridge (ThermoFisher). The filtration strategy separated bacteria into 2 fractions: a large (>3 μm) fraction, considered as bacteria that were attached to particles or possibly conglomerated, and a small (0.22–3 μm) fraction, considered as free-living, nonattached bacteria (Mohit et al., 2014). All filters were preserved in RNAlater™ (ThermoFisher) for 30 min in the dark, flash-frozen in liquid nitrogen and then stored at –80°C until their nucleic acid extraction in the laboratory. DNA and RNA were extracted from all the samples using AllPrep DNA/RNA Mini Kit (Qiagen, Hilden, Germany). Retrotranscription of the extracted RNA was carried out using the High-Capacity Reverse Transcription kit (ThermoFisher). We amplified the V4 hypervariable region of the conserved 16S ribosomal RNA gene by PCR using universal primers (Forward (515F-Y) 5′ – GTG YCA GCM GCC GCG GTA A – 3′, Parada et al., 2016; Reverse (806RB) 5′ – GGA CTA CNV GGG TWT CTA AT – 3, Apprill et al.,2015).

Negative PCR controls were run along with samples to detect any contamination during the process. All DNA concentrations were measured spectrophotometrically using a Nanodrop or equivalent system at the Institut de Biologie Intégrative et des Systèmes (IBIS) genomic analysis platform, Université Laval (Québec, QC, Canada). The target minimal amplicon concentration of 20 ng µL^–1 after the first PCR amplification was purified using AMPure XP paramagnetic beads, and amplicon quality was verified by electrophoresis (1% agarose gel). A second PCR amplification was carried out on the purified PCR amplicons, which were tagged for multiplexing with MiSeq® specific linking primers. The metabarcoded amplicons were purified as above, and equimolar concentrations of amplicons were sequenced using an Illumina MiSeq system at the “Plateforme d’Analyses Génomiques” (IBIS, Université Laval, Quebec, Canada).

The bioinformatic pipeline to create the OTU table from all samples combined was carried on the Graham supercomputer (Compute Canada). Sequence quality was verified by using Fastqc program (Andrews, 2018). Forward and reverse reads were merged by using BBMerge algorithm (Bushnell et al., 2017). The merged sequences were filtered, by using VSEARCH program, to remove Illumina adaptors, rare sequences (<2 occurrences), and low-quality sequences over 300 bp (fastq_maxee = 0.5; Rognes et al., 2016). The OTU clustering, with a 97% similarity threshold, and chimera removal were performed by USEARCH (Edgar, 2010). Taxonomic assignations were carried out on the raw sequences in mothur (Schloss and Westcott, 2011) against the SILVA v132 (Quast et al., 2013) reference database. Chloroplast, mitochondrial, and metazoan sequences were removed manually. The main OTU table was divided into subtables to analyze the Nelson and Churchill transects separately. These 2 subtables were divided again to separate rRNA from rDNA samples. OTUs with a total occurrence <4 in all the samples from each subtable were deleted from the analysis. During the manipulation, there were technical errors with 2 samples (Nel_RNA_S_NE-A and Nel_RNA_L_NE-NS_bt), which were not included in the downstream analysis.

Statistical tests were then conducted in the R environment (R Core Team, 2019). To evaluate the dissimilarity level between all the samples, the Bray-Curtis dissimilarity index, based on OTU presence and relative abundance in each sample, was calculated using the vegan package, on the rDNA sequences (Oksanen et al., 2019). Specifically, to correct the influence of rare species, data were standardized by using the Hellinger transformation beforehand from the vegan package. A hierarchical cluster analysis was performed on the Bray-Curtis indexes by using the option “ward.D2” in the hclust function (R software). Relative abundance of taxonomic classes in each sample was calculated with vegan and dplyr packages (Wickham et al., 2019).

To identify environmental variables with the most influence on the composition and distribution of bacterial communities along the 2 transects, Spearman correlation tests were performed between salinity, temperature, depth, and Chl a, nitrite, nitrate, silicate and phosphate concentrations. Variables that covaried were defined as those with a Spearman correlation coefficient ≥0.7; in the case of covariables, only one was included in the adjusted R² calculations. The final variables tested were depth, salinity, phosphate, and nitrite concentration. The adjusted R² measures the unbiased amount of explained variation and was used to select significant variables using a forward selection and 499 ANOVA permutations. The adjusted R² was calculated with the ordiR2step function with the 1000 permutations in the R2permutations argument and a P value of 0.05 (vegan package; Blanchet et al., 2008). The Akaike information criterion (AIC) coefficient and the adjusted R² were used to find the most parsimonious model. A distance-based redundancy analysis (dbRDA) was conducted using the Hellinger standardized rDNA OTU subtable and the selected environmental variables (vegan package). NE-NS samples were not included in the dbRDA analysis because environmental data were not available.

To identify indicator species from different size fractions and habitats (sites) defined by the hierarchical cluster analysis, we used the indicator value index IndVal (Dufrêne and Legendre, 1997), calculated with the indicspecies package on relative abundances of rDNA OTU reads in each sample and the multipatt function “IndVal.g” with 999 permutations (Cáceres and Legendre, 2009). Hayes River (HA-A) samples were excluded from this analysis. Only the OTUs with a P value ≥0.05, a total mean relative abundance of ≥0.1% in the 3 habitats and presence in rRNA libraries from at least one sample in the corresponding habitat were considered. Indicator OTUs were considered specialists to a specific habitat, according to a ratio ≥0.8 of the mean relative abundance in a single habitat (riverine, estuarine, or marine) to the total mean relative abundance in all 3 habitats. Indicator species characteristic of more than one habitat were also identified following the same approach. The presence in at least one of the corresponding habitats in the rRNA data confirmed that the individual indicator OTUs were active in at least one habitat and not part of a larger detrital pool. A similar analysis was used to identify the size fraction indicator species. In the case of the size fraction indicator OTUs, we also verified if the same OTU was present in the corresponding rRNA data.

All the graphics were created within the R environment (R Core Team, 2019). The bar plots and bubble plots were generated in the ggplot2 package, ternary graphics, and Venn diagrams in ggtern and VennDiagram packages, respectively (Wickham, 2016; Hamilton and Ferry, 2018; Chen, 2018).

Nutrient and Chl a concentrations, as well as salinity and temperature including ship-based CTD vertical profiles, have been reported elsewhere (Jacquemot et al., 2021), with the data available at the Canadian Watershed Information Site (2018). Site-specific values for salinity, temperature, and Chl a, and nutrient concentrations for this study are given in Table S1.

Briefly, for the Churchill system, samples from furthest up-river station CH-A were freshwater (salinity of 0.1) and the temperature was 10.5°C. Temperature gradually decreased toward offshore (Figure S1), and the salinity gradient went from 10.9 to 30.7 over less than 20 km from the CH-B to CH-E stations (Figure 1B). Marine waters (salinity values >29) were detected below the thermocline, which also corresponded to a halocline with colder and saltier water beneath warmer fresher surface water. The Nelson River waters at the nearshore station NE-A were fresh with a salinity of 0.2 and had the highest water temperature recorded during the study (21.9°C). The adjacent Hayes River sample at station HA-A was also fresh (salinity of 0.1) and warm (20.7°C). The highest temperatures should be viewed with caution, as samples from surface water at the river edge might not have been representative of the main channel due to slack edge water and surface solar heating. The Nelson salinity gradient (0.1–31.4), as well as the temperature gradient from 14.3°C for river station NE-B sampled the same day from barge to 0.2°C, highlighted the more extensive freshwater plume along the Nelson transect, spreading more than 60 km away from the river mouth (Figure S1). The deepest station along this transect was NE-46 (Figure 1B and C).

Nutrient concentrations varied somewhat between the 2 study regions. Briefly, for the Churchill transect, highest concentrations for nitrite + nitrate were 0.32 µmol L^–1 at the river station CH-A. Nitrate was below the level of detection (0.02 μmol L^–1) at the surface for the remaining stations, with nitrite at <0.03 μmol L^–1. In contrast, phosphate ranged from 0.32 to 0.51 μmol L^–1 for all stations, except CH-A, with a concentration of 0.05 μmol L^–1. The silicate concentrations were generally high, especially at CH-A and in the surface sample from CH-B (18.11 and 17.72 μmol L^–1, respectively), compared to the coastal stations CH-C, CH-D, and CH-E which were <6.50 μmol L^–1 (Table S1). Compared to the Churchill transect, the Nelson transect had higher inorganic nitrogen concentrations at the upstream station NE-A, with a maximum of 0.44 μmol L^–1 nitrite + nitrate in the Nelson River. The maximum silicate concentration value was 38.5 μmol L^–1 at NE-A and decreased seaward but remained relatively high. In contrast, the phosphate concentrations were greater for the offshore stations, with a maximum value of 0.58 μmol L^–1 at NE-46. The Hayes River nutrient concentrations were lower than in the Nelson River at 0.13 μmol L^–1 for nitrite + nitrate, 34.6 μmol L^–1 for silicate, and 0.10 μmol L^–1 for phosphate (Table S1).

There were no trends across habitats in the total bacterial FCM concentrations, which varied from 0.9 to 5.3 × 10⁶ cells mL^–1 with the highest concentration from below the thermocline at the NE–WE3 station (Table S1). As the total bacterial FCM data did not satisfy all of the requirements to perform an ANOVA analysis to compare the FCM counts between the different habitats along the 2 transects, no further analysis was carried out.

After removing chimeras, nontarget sequences and rare OTUs (<4 occurrences) from the rDNA OTU table, we found 7,195 OTUs at a similarity level of 97%, of which 3,494 were in the Churchill transect (18 samples) and 2,768 were in the Nelson transect (42 samples) with 485 from the Hayes River (2 samples). In the rRNA subtables, we defined 2,575 OTUs in the Churchill transect subtable (18 samples) and 2,811 OTUs in the Nelson transect subtable (40 samples). Raw sequence data for both rRNA and rDNA are available at NCBI (2021) under PRJNA680849. Fasta files for OTUs mentioned in the text or tables are provided in Data File S1. The hierarchical cluster analysis, based on OTU presence and relative abundance in each sample, showed rDNA samples clearly grouped into 3 habitats: riverine, estuarine, and marine (Figure 2). The topology was evident in the distribution of diverse freshwater and marine representatives in the Gammaproteobacteria (Figure S2), which now encompasses groups previously classified as Betaproteobacteria (Hug et al., 2016). Among the 3 major habitats, some clustering was observed by size fraction within the marine samples from both the Churchill and Nelson transects and in the estuarine samples from the Nelson transect (Figure 2).

Significance of the environmental variables in determining community differences was tested using the ordination model by permutation test (ordiR2step function). Based on adjusted R² and P values (P values = 0.002), salinity, phosphate, and nitrite concentrations were significant. The AIC coefficient and the adjusted R² indicated that the most parsimonious model involved nitrite concentration (Table S2). A high F value in the salinity model indicated greater variation between the bacterial communities, accounting for 51.2% of the total variation (Figure 3A). Within each habitat (riverine, estuarine, marine), the samples were also separated corresponding to their geographic origin (Nelson versus Churchill; Figure 3A), with more unique OTUs along the Churchill transect (Figure 3B). Bacterial communities from the Hayes River were more similar to the ones from the Churchill River than the Nelson River (Figure 3A), with more shared OTUs (Figure 3C).

The bacterial communities included indicator species defined using IndVal for riverine, estuarine, and marine sites (Figure 4, Data File S2). For the Churchill transect, the most exclusive (specialist) indicators of a single habitat were seen near the apexes of the ternary plot (Figure 4A). Overall, 84 OTUs were classified as Churchill River indicators, 30 as estuarine indicators, and 61 as marine indicators. In addition, there were 22 OTUs that were indicative of both riverine and estuarine habitats and 29 indicators of marine plus estuarine sites (Figure 4A, Data File S2), but no indicators of riverine plus marine sites. These findings suggest that species in the estuary were nearly equally from the river, from marine sites and autochthonous in the Churchill region at the time of sampling.

The pattern from the Nelson system was similar for the single site indicators (Figure 4B, Data File S2), with 114 OTUs classified as Nelson River indicators, 26 as estuarine indicators, and 60 as marine indicators. There was a difference in the numbers of indicators for 2 sites, with only 6 OTUs indicative of the river plus estuary and 145 indicators of marine plus estuarine sites (Figure 4B, Data File S2). As with the Churchill system, there were no OTUs indicating river plus marine sites. The low number of river plus estuarine indicators and higher proportion of shared indicator for the estuary plus marine sites suggest that the Nelson estuarine community was mostly recruited from the marine system.

Cyanobacterial distributions were striking among the riverine and estuarine indicators (Figure 4), with a progressive OTU change along the river-to-estuary continuum in the Churchill system (green symbols in Figure 4A). The cyanobacteria from both rDNA and rRNA in the Nelson systems showed a logarithmic decrease with increasing salinity with a significant negative correlation with salinity (τ= –0.56, P value = 4.9 × 10^–12; Figures 5A and S3). FCM counts of the PC-rich cyanobacteria in the samples from the Nelson River transect mirrored the read counts (τ= 0.53, P value = 6.8 × 10^–11). The percentage of cyanobacterial reads per sample and the FCM counts from the Churchill transect tended to decrease with increasing salinity; however, the low number of sites sampled precluded statistical analysis (Figure 5B). At the level of OTUs, the cyanobacteria, in particular OTU 1, OTU 217, and OTU 611, were relatively common in the Nelson and Churchill riverine samples and decreased with increasing salinity in both systems (Figure S3). In the Churchill samples, this pattern was most evident from the rRNA read proportions, with reads for the rDNA more persistent offshore (Figure 5B).

Overall, the bacterial community from rDNA in all 3 rivers was characterized by a high relative abundance of Actinobacteria, mostly Sporichthyaceae (Figures S4 and S5), Gammaproteobacteria, especially from the family Burkholderiaceae in the order Betaproteobacteriales (Figure S2), and Cyanobacteria (Figures S4 and S5). Other major groups such as Chitinophagales, Verrucomicrobia, and SAR11 clade III (in the Alphaproteobacteria) were present, but locally more common in the Nelson River (Figures S4 and S5). Based on the rDNA, the Actinobacteria were the most diverse (having more OTUs) of these abundant major groups (Figures 6A and S5). At the putative species (OTU) level, the smaller rivers both had higher proportions of Sphingobacteriales, Cyanobacteria, and Betaproteobacteriales (Figures 6A, S2, S4, and S5). Specifically, OTUs retrieved at higher proportions (by read count) in the Churchill and Hayes river samples compared to the Nelson included OTU 2670 (Polynucleobacter) and OTU 8 (Rhodoferax). Although Actinobacteria were in all rivers, 3 of the Actinobacteria OTUs were almost exclusive to the Nelson-Hayes system (OTUs 159, 4092, and 548) and 2 (OTUs 13 and 944) were mostly from the Nelson River. Some SAR11 clade III taxa also appeared preferentially in the Nelson, including OTUs 713 and 624 (Figure 6A).

Many of the abundant OTUs, Flavobacteriales, Betaproteobacteriales, and some SAR11 taxa were seen in both estuaries (Figures 6B, S2, S6, and S7). A relatively greater proportion of classic marine Gammaproteobacteria, including Oceanospirillales, Cellvibrionales, and Thiomicrospirales were found in the estuaries (Figure S2). Differences between the 2 estuaries followed the clustering based on whole communities (Figure 2), separating the Nelson and Churchill transects. The Churchill estuarine samples had a higher relative abundance of Cyanobacteria and Verrucomicrobia than in estuarine samples from the Nelson (Figures S6 and S7).

Overall, the marine bacterial community was represented by high relative abundance of Gammaproteobacteria (Oceanospirillales, Cellvibrionales, Thiomicrospirales) and Alphaproteobacteria (SAR11 clade 1a and Rhodobacteraceae; Figures 6C, S2, S8, and S9). Archaea were more commonly found in the Nelson transect (Figure S9), but no archaeal OTUs were among the 5 most “abundant” estuarine OTUs identified per transect (Figure 6C). There were differences in relative abundance of several specific OTUs; for example, OTU 110 (Rhodobacteraceae, but with a best BLAST match of 98% to the strain Sulfitobacter sp. Haha2106) appeared preferentially in the Churchill marine samples (Figures 6C and S9). In the Nelson, Polaribacter_1 spp. (OTU 963), unclassified Flavobacteriaceae (OTUs 539 and 544), SUP05 cluster (OTU 1525), and SAR11 clade 1a (OTU 581) were prominent in marine samples (Figures 6C and S9) and in the estuarine samples (Figures 6B, S6, and S7).

Community clustering (Figure 2) suggested that fractionation by size was pronounced in the Nelson marine and estuarine samples and for the Churchill marine samples. We applied the IndVal analysis to identify the indicator species for the 2 size fractions (Table 1, Data File S3). No size fraction indicators were identified in the rivers or the Churchill Estuary.

For the marine samples and the Nelson Estuary, 64% of large fraction indicators were Bacteroidia, mostly Flavobacteriales (Table 1, Data File S3). The small size fraction community included Pseudomonadales and other Gammaproteobacteria (59%) and Alphaproteobacteria (37%; Table 1). The size fraction indicator OTUs were also recovered in rRNA from more than 50% of estuarine rRNA samples from the Nelson transect. Within the marine sample size fractions, 4 of the indicator taxa were shared between the Nelson and Churchill transects (Data File S3), and 4 were at both the Nelson marine and Nelson estuarine sites. There were no indicators shared among all 3 habitats.

The bacterial communities separated into 3 defined habitats, with the Bray-Curtis clustering following similar topology reported using the abiotic variables and microbial eukaryotic communities reported in Jacquemot et al. (2021). In that study, salinity and nutrient concentrations separated the same sites into riverine, estuarine, and marine clusters, with riverine sites most like each other. The overall separation of 3 habitats is consistent with the sub-Arctic estuaries being similar to lower latitude regions, with the development of characteristic estuarine communities (Crump et al., 1999). The biological difference between the 2 systems (Nelson versus Churchill) was more pronounced compared to Jacquemot et al. (2021), who found that surface microbial eukaryotes were similar in both estuaries. In addition, below the pycnocline, the Churchill microbial eukaryotes were characteristic of the more saline waters offshore. In contrast, our bacterial IndVal analysis showed that within the transition zone, communities separated mostly according to geography with little difference between depths. The differences between studies targeting 2 different domains of life could suggest greater sensitivity to modest salinity changes among eukaryotes compared to bacteria. Alternatively, the difference between eukaryotes and bacteria could be an artifact, with the 16S rRNA gene for bacteria less sensitive in detecting ecotypes compared to the 18S gene used for microbial eukaryotes. Even considering these potential biases, the overlap of marine and estuarine bacterial species suggests multiple euryhaline taxa able to tolerate wide salinity fluctuations (Figure 4). Having a wide salinity tolerance is in keeping with seasonal exposure to sea-ice melt and lower average salinity in the Arctic, and is consistent with the uniqueness of the mostly coastal Arctic interactome compared to the global microbial interactome (Chaffron et al., 2021).

One striking pattern was the differences at all taxonomic levels including phyla in the riverine samples compared to the marine samples (Figure 5). Most of the OTUs in the diverse river phyla were associated with freshwater species reported elsewhere, for example, Polynucleobacter spp., Rhodoferax, Actinobacteria (species in the Sporichthyaceae), Chitinophagales, and Candidatus Planktophila (Ghai et al., 2014; Coenye, 2014; Salcher et al., 2015; Cabello-Yeves et al., 2017; Tee et al., 2020). Other indicators of both the rivers and Churchill estuary were taxa from alpine and polar rivers, including Verrucomicrobiae and the actinobacterium Canditatus Planktophila limnetica (Karlov et al., 2011).

Polynucleobacter in the ecologically diverse family Burkholderiaceae (Coenye, 2014) includes free-living freshwater species and endosymbionts of the Euplotes, a genus of ciliates. However, Jacquemot et al. (2021) did not identify Euplotes in microbial eukaryote communities, and the higher relative abundance of Polynucleobacter in the small size fraction, especially in the Churchill and Hayes Rivers (Figure 5), suggests that it was a free-living species as reported for the Columbia River system (Crump et al., 1999). The photoheterotrophic, potentially nitrogen-fixing genus Rhodoferax was retrieved from the riverine habitat, but not in the estuaries suggesting that the species here was a freshwater ecotype.

Actinobacteria, mostly seen in the Nelson River, were primarily in the family Sporichthyaceace and the uncultivated CL500-29 clade. Actinobacteria were the least likely group to be detected in rRNA, which could be attributed to slow growth rate (Normand, 2006) or indicate passive transport or ground water infiltration and not in situ growth as seen for the cyanobacteria (Figure 5). Sporichthyaceae are reported mainly from terrestrial habitats, suggesting possible recruitment of these bacteria from the soil (Tamura, 2014), which could be consistent with higher flows and transport of subsurface water from the upstream catchment into the river within the Nelson system compared to the Churchill. CL500-29 has been reported from brackish regions of the Baltic Sea associated with relatively high quantities of terrestrial dissolved organic matter concentrations (Lindh et al., 2015). Other potential soil bacteria in the Nelson River included Terrimicrobium spp. (phylum Verrucomicrobia). To date there is only one described species in this genus, Terrimicrobium sacchariphilum, which is anaerobic and lives mainly in soils (Qiu et al., 2014). The similarity of the Hayes River and Churchill bacterial assemblages (Figures 2, 3, and S2) is in keeping with reduced prevalence of soil bacteria in the smaller, low flow rivers. Overall, the prevalence of potential soil microbes in the Nelson River with its larger catchment and higher flow rate suggests more connectivity between the river water and soils (Ruiz-González et al., 2015) compared to the Churchill.

Within the Nelson River samples were several Alphaproteobacteria SAR11 clade III, which were classified as indicators of the Nelson River but also were found in the Churchill estuary (Figures 6, S5, and S7). This distribution supports the notion that this clade could be euryhaline as suggested by Kraemer et al. (2020), who characterized SAR11 clades from the Beaufort Sea and found that the SAR11 P3.2 phylotype, belonging to SAR11 clade III, was only recovered in fresher Arctic surface water samples.

The Churchill estuarine bacterial community indicator taxa belonged mainly to Cyanobacteria, Verrucomicrobia, and SAR11 group III spp., which had riverine origins (Figure 4A). Based on the IndVal analysis, there were fundamental differences between the Nelson River estuary and the Churchill estuary, with the Nelson having few members of an estuarine bacterial community (Figure 4B) and instead a bacterial assemblage composed mainly of presumably euryhaline and marine species (Figures 4B and 6, Data File S2). In contrast, the Churchill Estuary included more riverine than marine taxa. Fortunato et al. (2013) noted that the Columbia River estuarine bacterial community was dominated by riverine taxa, which they attributed to a high riverine flow rate and short water residency time. In this case, although the Nelson River has a higher riverine input than the Churchill River, the presence of the marine-influenced bacterial community suggests that the primary selective forces were a longer residency time (Gueguen et al., 2016) and more gradual transition from fresh to marine waters favoring marine taxa brought into the estuary via estuarine circulation (Jacquemot et al., 2021).

As mentioned above, a number of river indicator species were found in the Churchill estuarine environment, but not in the offshore marine habitat in the Churchill system. We suggest that the low water flow rate, caused by the diversion of the Churchill River, combined with its semi-enclosed estuary type limits mixing between the riverine and marine waters, increasing stratification in the transition zone. The sharp salinity and temperature gradients between the estuarine and marine habitats could increase the probability of osmotic shock, limiting physiological adjustment via a stress response (Jayakumar et al., 2020) and ultimately resulting in the lack of brackish river taxa in the marine system. Tidal activity could reduce the exchange of taxa between the riverine and the marine environments during a slack tide in particular. In the Churchill region, the amplitude of the tide is approximately 145 cm (Ray, 2016); thus the effect of the tide on taxa exchange between the riverine and marine habitats was likely negligible as well. We sampled the systems over several days including from different tidal regimes along the Nelson, and the similarity of the estuarine samples suggests that tidal influence was minimal in that estuary.

Although at a community-wide clustering level the 2 estuaries differed, several indicator taxa were shared between estuaries and the marine samples, including Flavobacteria (unclassified and Polaribacter). Flavobacteria are found in many ecological niches and habitats, such as soils, sediments, and aquatic habitats. The marine members of the family Flavobacteriaceae are often found at the water surface, and some of these species could play a role in the degradation of macromolecules when they are associated with macroalgae, phytoplankton, and particulate material (McBride, 2014). Indeed, Krüger et al. (2019) reported that Polaribacter sp. and NS3a marine group sp. are involved in polysaccharide degradation during North Sea spring phytoplankton blooms.

The slightly fresher, colder offshore waters in the region of the Nelson transect (Table 1, Figure 1B) were due to more recent sea-ice melt (Jacquemot et al., 2021). “Pre-adapted” taxa retained by estuarine circulation would have been more competitive than river-sourced purely freshwater species, which were exposed to higher temperatures from the upstream summer conditions.

The marine taxa shared by both the Churchill and Nelson sites were mainly in the Flavobacteriales, Alphaproteobacteria (including SAR11 clade 1a spp.), and Gammaproteobacteria (Figures 4, S2, and S9) that are typical of marine environments (Seo et al., 2017; Kraemer et al., 2020). SAR11 clade 1a seems to be common in cold marine waters (Brown et al., 2012; Kraemer et al., 2020) and was relatively abundant in the surface layer of our marine samples in both transects. However, Kraemer et al. (2020) reported that SAR11 clade 1a was nearly absent from the surface of the Beaufort Gyre in the Western Arctic Ocean. As we did not sample below 11 m, deeper offshore profiles would be needed to characterize depth distributions in Hudson Bay. The timing and nutrient characteristics could also play a role, given that the Kraemer et al. study was carried out late in the season and in a much more stratified system. Although other relatively abundant marine OTUs included uncultured Nitrincolaceae, some of which are able to carry out heterotrophic nitrate reduction to ammonium, the majority of dominant taxa were aerobic heterotrophic taxa.

The low nitrite + nitrate concentrations in the offshore Churchill stations are consistent with an earlier algae bloom where phytoplankton would have drawn down inorganic nitrogen, including nitrate (Bristow et al., 2017; Matthes et al., 2021), and the microbial eukaryote community in the Churchill offshore samples was dominated by post-bloom heterotrophic taxa in a companion study (Jacquemot et al., 2021). In contrast, the Nelson marine samples were dominated by a diatom Rhizosolenia spp. blooming at the same time (Jacquemot et al., 2021). The contrast provided an opportunity to examine the bacterial communities in the post bloom (Churchill) and ongoing bloom (Nelson) marine sites. However, there was no consistent separation by the dominant taxa, with the possible exception of post-bloom Churchill marine samples having a higher proportion of the OTU 110 (Rhodobacteraceae with affiliations to Sulfitobacter sp.) compared to Nelson marine samples. Some Sulfitobacter spp. oxidize sulfite (Pujalte et al., 2014), and in the Arctic another strain is associated with inducing algal cell death and degradation of DMSP (Zeng et al., 2020), consistent with a post-bloom environment.

Among marine taxa also found in their corresponding estuarine bacterial community were taxa in the family Thioglobaceae that are reported from diverse habitats from the surface ocean to the lower oceanic crust (Li et al., 2020) and deep-sea vent animals (Georgieva et al., 2020). Thioglobaceae were indicators of both the estuarine and marine small size fraction, suggesting that at least in Hudson Bay, this family is dominated by free-living representatives (Nowinski et al., 2019; Li et al., 2020). Thioglobaceae were the most abundant taxa during a coastal autumn phytoplankton bloom in Monterey Bay, CA, USA (Nowinski et al., 2019), in keeping with the post-bloom environment of western Hudson Bay at the time of sampling (Matthes et al., 2021).

Because the majority of the reads of the Cyanobacteria were assigned to Cyanobium spp., which contains PC and not PE (Larsson et al., 2014), we are confident that the FCM counts were an accurate estimation of the total cyanobacteria. Cyanobacteria in our samples from the Churchill and Nelson transects followed a similar pattern to that found for the Mackenzie River–Beaufort Sea system (Waleron et al., 2007; Vallieres et al., 2008). In that study, 16S rDNA picocyanobacteria offshore were freshwater species likely transported by the Mackenzie River’s large fresher plume, with the same freshwater cyanobacteria found in some fresher offshore samples, but no cyanobacteria were found in the contiguous upwelling marine waters. Cyanobium is usually found in freshwater (Ernst et al., 2003; Shih et al., 2013) and in Hudson Bay the rRNA relative abundance of Cyanobium (Figure S3) decreased with increasing salinity in both transects. The sharper salinity gradient in the Churchill system may also have contributed to the rRNA to rDNA pattern of OTUs. In particular, the rDNA signal was more persistent moving away from the river and in contrast to the rapidly decreasing rRNA signal, which would be consistent with dead cells remaining in the water column (Blazewicz et al., 2013).

The distinct free-living bacterial community in the small fraction included SAR11 clade III spp. (Table 1), in keeping with the small size and presumed ecological niche of this clade as free-living taxa (Morris et al., 2002; Giovannoni, 2017). Thioglobaceae were indicators of small size fraction as reported elsewhere (Nowinski et al., 2019; Li et al., 2020). In the Churchill and Nelson marine samples, SAR11 clade Ia spp. were indicators of the small fraction in agreement with known distributions and ecology (Kraemer et al., 2020).

The Sphingobacteriales NS11-12 marine group spp. was an indicator for the large fraction (Table 1). The NS11-12 marine group (Sphingobacteriales) was previously reported to correlate positively with contaminants (Pb and Cu) in an urbanized marine coastal area (Coclet et al., 2019). While the full ecological implications remain unknown (Meziti et al., 2015), the association with the large fraction could be due to suspended inorganic particles, rich in minerals, and is consistent with seasonally acquired trace metals in polar bears in Western Hudson Bay (Bechshoft et al., 2020). Other marine large fraction bacteria included several Ulvibacter spp., which have been associated previously with diatoms and marine snow (Teeling et al., 2012). Smith et al. (2013) also noted Ulvibacter in the particulate fractions of the plume and hypoxic water samples from the Columbia River.

The size fraction indicator species analysis of the Nelson estuarine samples suggested several estuarine-specific particle-attached taxa (Data File S3). The most common were a Sphingobacteriales NS11-12 marine group species, the Actinobacteria (hgcl clade), and a Flavobacteriales related to uncultured Crocinitomicaceae, which is a chemoorganotrophic marine or freshwater family (Bowman, 2020). The Actinobacteria hgcl clade is normally associated with freshwater (Aguilar et al., 2019) and, considering the lower representation in the rRNA libraries (Supplementary Figure S5), could have been an indication of residual freshwater inflow like the cyanobacteria tracer.

Previously, in the Mackenzie River system in the western Arctic, Garneau et al. (2009) found dissimilar free-living and particle-attached bacterial communities in areas of high particulate organic material, which provides both a substrate for attachment and a source of complex organic matter that can be degraded by extracellular enzymes and used for growth (Vetter and Deming, 1999). In the estuary, these attached bacteria were likely a resource for the heterotrophic microbial eukaryotic community (based on 18S rDNA; Jacquemot et al., 2021) that rely on bacteria associated with the organic particulate material and suspended sediments. The post-bloom Churchill system would have resulted in the third type community of attached bacteria observed in the Churchill marine samples. In summary, bacteria in the Nelson Estuary and coastal area would have exploited particulates retained due to estuarine circulation and, on the seaward side, particulates associated with a phytoplankton bloom (Villacorte et al., 2015; Jacquemot et al., 2021).

The mostly marine bacteria found in the Nelson Estuary and high proportion of freshwater bacteria in the Churchill Estuary contrasts with the distinct estuarine protist communities reported by Jacquemot et al. (2021). The difference suggests that there may be additional environmental selective forces or biological constraints on bacterial communities. For example, given the high proportion of bacterial grazers in estuaries (Lovejoy et al., 1993; Jacquemot et al., 2021) there could be top-down control of bacteria with less opportunity for distinct communities to develop. In addition, phages could also act as agents influencing bacterial diversity in high latitude waters (Vaque et al., 2017).

The transition between freshwater and marine ecosystems can be abrupt or more gradual, depending on estuarine circulation, with the degree of estuarine circulation influenced by the magnitude of river flow, the underlying geomorphology and tidal cycles determining residence time (Frenette et al., 1995). Although hydrological data across most of the Arctic and sub-Arctic are limited, the Nelson and Churchill rivers are somewhat exceptional because of a history of flow regulation and monitoring by Manitoba Hydro beginning in 1975 (Newbury et al., 1984; Déry et al., 2011). The Churchill and Nelson rivers, which differ in terms of discharge volume and underlying geomorphology, are representative of 2 types of estuaries following the Pritchard estuarine classification (Pritchard, 1989). The Churchill Estuary, at the end of a more confined channel, is a stratified estuary with low riverine inputs and weak vertical mixing (Figures 1 and S1) and here led to a distinct river-influenced estuarine community. In contrast, the Nelson River fans out upon reaching Hudson Bay and, due to the hydroelectric diversions, has a comparatively greater discharge than the Churchill River. The Nelson Estuary is a partially mixed type estuary in the Prichard classification (Pritchard, 1989). During the time of sampling, a maximum turbidity zone was evident in the Nelson Estuary (Jacquemot et al., 2021) and, due to the higher river input, terrestrial dissolved organic matter enters and is retained in the estuary (Gueguen et al., 2016). In addition, estuarine circulation would ensure that the organic particulate material and suspended sediments, which are abundant in the Nelson system (Godin et al., 2017), would be suspended in the estuary. The higher particulate load and organic matter in the maximum turbidity zone in the Nelson Estuary and the Rhizosolenia-dominated phytoplankton bloom associated with marine samples would have added suitable substrates for specialized attached bacteria in both habitats.

Our study is based on a snapshot of a single time point for each site during a single season of the year. Given the high seasonal variability and potential for tidal and daily discharge interactions, research is needed across multiple temporal scales to fully describe the bacterial communities from the Nelson and Churchill rivers. Given that the flow rate of the Nelson is controlled with little seasonal change, we would predict less variability in the Nelson compared to the Churchill over the annual cycle, but this conjecture remains to be tested.

The bacterial communities separated into the 3 salinity-defined habitats: riverine, estuarine, and marine (Figure 2). Nutrient concentrations were also statistically associated with the community dissimilarities (Figure 3A, Table S2) contributing to the separation of the Churchill and Nelson systems. These results are in keeping with bacterial communities as sensitive indicators of environmental conditions that can provide insight into land, freshwater, and marine connectivity (Comte et al., 2018; Papale et al., 2020). In remote regions such as the Arctic and sub-Arctic, the study of processes that influence estuarine circulation, and thus whether the transition from freshwater to marine ecosystems is abrupt or gradual, has been limited due to poor accessibility and the difficulty of maintaining standard tidal and hydrological instrumentation. Moreover, the microbiological implications of the Churchill diversion into the Nelson including up to the present day are not known. This lack of knowledge and research can be ascribed partially to the engineering perspective that considers only macro effects of environmental perturbations and partially to a lack of suitable tools to resolve microbial diversity. The recent widespread adoption of DNA sequencing directly from the environment is set to overcome this second obstacle. Here we establish the utility of sequencing data to identify similarities and differences in bacterial communities in 2 very different river-to-sea transition zones in Hudson Bay.

This study provides the first detailed portrait of the composition and spatial distribution of bacterioplankton communities using high throughput amplicon sequencing in the Churchill and Nelson rivers, estuaries and nearshore Hudson Bay, and identifies the main environmental factors shaping these bacterial assemblages. Salinity and nutrient conditions in the Nelson and Churchill transects were associated with specific bacterioplankton communities. Churchill estuarine bacterial communities were mainly composed of likely halotolerant species from allochthonous river inputs and included habitat indicator species belonging to Cyanobacteria, Verrucomicrobia and SAR11 group III spp. The low flow rate, caused by the diversion of the Churchill River combined with its semi-enclosed estuary type, limits mixing between the riverine and marine waters, increasing stratification in the transition zone, and resulting in separate estuarine and marine communities. In contrast, the Nelson estuarine indicator species were marine species, consistent with retention and recirculating marine waters that were in balance with the higher flow rate of the regulated Nelson River. Furthermore, in the Nelson Estuary and in the marine habitat from both transects, the particle-attached and free-living bacterial communities differed, suggesting sufficient resources to maintain different niches.

From this first portrait of the bacterial communities from the regulated Nelson and Churchill rivers, the differences between estuaries suggest that modified flow regimes influenced the resulting bacterial communities. Although sampling intensity was lower in the much smaller Churchill system, compared to the larger Nelson, environmental factors in both estuaries were associated with distinct communities along the river-to-sea transition. In any estuary there is no fixed boundary between salt and freshwater, but the coupling between salinity and community provides a measure of how species are selected in the environment. This snapshot of 2 river-to-sea transitions highlights how estuaries can differ from each other and the need for more seasonal and interannual studies. Although the significance of the patterns observed in the context of ongoing climate change remains to be addressed, we speculate that climate change will contribute to an increasing frequency of rare events that will impact the timing and magnitude of river discharge and flow, and that bacterial community structure will reflect these impacts.

Raw paired-end reads have been deposited on NCBI under the BioProject accession number PRJNA680849 (https://www.ncbi.nlm.nih.gov/sra/PRJNA680849). CTD/Rosette data used in this study are publicly available on the Canadian Watershed Information Network (CanWIN): https://doi.org/10.5884/12713.

The supplemental files for this article can be found as follows:

Table S1. Environmental properties at the time of sampling.

Table S2A. Results of the ordination model permutation test performed on the variables used in the forward selection.

Table S2B. Biplot scores of dbRDA1 to dbRDA6 axes for constraining variables (dbRDA).

Table S2C. Spearman’s rank correlation coefficients calculated on environmental variables from the Nelson, Hayes, and Churchill available samples (bold values are >0.7).

Data File S1. Fasta files for OTUs mentioned in article.

Data File S2. Indicator species in each habitat from the Nelson and Churchill transects. The Index specifies the habitat of the significant OTU (River, 1; Estuary 2; Marine, 3; River + Estuary, 4; river + Marine, 5; Estuary + marine, 6). The statistic value (Stat) is the IndVal index. The sum of the mean relative abundances (%) of 1 OTU in the 3 habitats (Total_abund_%). The ratio of the mean relative abundance of the selected OTU in 1 habitat to the sum of mean relative abundances in the 3 habitats is represented in percent by “River_%”. “Estuary_%” or “Marine_%”.

Data File S3. Indicator species of the size fractions. Marine and estuarine samples from the Nelson transect and Marine Samples from the Churchill transect. The statistic value (Stat) is the IndVal index (func = “IndVal.g”), which represents the association strength between an OTU and its size fraction (0 = no association; 1 = strong association), and P value indicates significance level. The sum of the mean relative abundances (%) of one OTU in the 2 size fractions, from the estuarine samples (MRA%). The ratio of the mean relative abundance of the selected OTU in its corresponding size fraction to the sum of mean relative abundances in both size fractions is represented in percent by “Proportion.” Only OTUs with a Stat value ≥ 0.8, a P value ≥ 0.05, a MRA value ≥ 0.1%, and a Proportion value ≥ 65% have been selected.

Figure S1. Vertical temperature profiles from CTD data. The black dots with a white border indicate the position of the samples in the water column. Nelson transect is on the left (NE stations); Churchill transect, on the right (CH stations). Figure created with Ocean Data View Software.

Figure S2. Relative rDNA abundances of major categories of Gammaproteobacteria. The order was defined by the hierarchical cluster analysis based on the Bray-Curtis indices from the rDNA samples (Figure 1, main text). Large fraction samples are in solid symbols; small fraction, in open symbols.

Figure S3. Relative abundances of the 3 most abundant cyanobacteria OTUs along the 2 transects. Top panels are from the Nelson transect; lower panels are from the Churchill transect. All values are log10 of the reads from rDNA and rRNA.

Figure S4. Relative rDNA abundances of major bacterial groups in the riverine cluster samples. The riverine samples were defined by the hierarchical cluster analysis based on the Bray-Curtis indices from the rDNA samples (Figure 1, main text).

Figure S5. Most abundant OTUs (percentage of reads ≥1%) in Churchill, Hayes, and Nelson riverine samples. The riverine samples were defined by the hierarchical cluster analysis based on the Bray-Curtis indexes from the rDNA samples. OTUs not found in a given sample rRNA library (percentage of rDNA reads ≥1%) are indicated by an asterisk in the circle if not found in the specific sample. The “§” symbol above the Nel_DNA_S_NE-A_sf sample indicates that its corresponding rRNA sample was not available. A + sign in the habitat symbol on the right indicates that although the rDNA classified it as an indicator, the taxon was not recovered from the rRNA from that habitat.

Figure S6. Relative rDNA abundances of major bacterial groups in the estuarine cluster samples. The estuarine samples were defined by the hierarchical cluster analysis based on the Bray-Curtis indexes from the rDNA samples (Figure 1, main text).

Figure S7. Most abundant OTUs (percentage of reads ≥1%) in Churchill and Nelson estuarine samples. The estuarine samples were defined by the hierarchical cluster analysis based on the Bray-Curtis indexes from the rDNA with percentage of reads ≥1%. OTUs not found in a given sample rRNA library (percentage of rDNA reads ≥1%) are indicated by an asterisk in the circle.

Figure S8. Relative rDNA abundances of major bacterial groups in the marine cluster samples. The marine samples were defined by the hierarchical cluster analysis based on the Bray-Curtis indexes from the rDNA samples (Figure 1, main text).

Figure S9. Most abundant OTUs (percentage of reads ≥1%) in Churchill and Nelson marine samples. The marine samples were defined by the hierarchical cluster analysis based on the Bray-Curtis indexes from the rDNA samples. The “§” symbol above the Nel_DNA_L_NE-NS_bt sample indicates that its corresponding rRNA sample was not available.

We thank P. Guillot for processing CTD data, and G. Deslongchamps and J. Gagnon for collecting and processing nutrient data. Our gratitude goes to the Plateforme d’Analyses Genomiques (IBIS, Université Laval) for advice and sequencing support; data analysis was carried out using Compute Canada facilities. The CTD data were collected onboard the Canadian research icebreaker CCGS Amundsen as part of the Hudson Bay System Study (BaySys). Logistical support was provided by the Amundsen Science program, which is supported by the Canada Foundation for Innovation. We are grateful to the Canadian Coast Guard officers and crew of the CCGS Amundsen for logistics.

This research is a contribution to the BaySys project funded by Natural Science and Engineering Council of Canada (NSERC) and Manitoba Hydro with additional NSERC support in the form of an undergraduate student research award (USRA) to CM and discovery grants to CL. Additional funding came from Fonds de recherche du Québec Nature et Technologies (FRQNT) to Québec Océan, and ArcticNet, a Canadian Center of Research and Excellence. CTD data were made available through the Amundsen Science program, supported by the Canada Foundation for Innovative and NSERC.

The authors have no competing interests to declare.

• Contributed to the conception and design: CM, LJ, CL.

• Collected samples: LJ.

• Contributed to analysis and interpretation of data: CM, MP, LJ, CL.

• Wrote the manuscript: CM, CL, with input from all authors.

• Approved the submitted version for publication: All authors.

How to cite this article: Morency, C, Jacquemot, L, Potvin, M, Lovejoy, C. 2022. A microbial perspective on the local influence of Arctic rivers and estuaries on Hudson Bay (Canada). Elementa: Science of the Anthropocene 10(1). DOI: https://doi.org/10.1525/elementa.2021.00009

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Associate Editor: Jeff S. Bowman, Integrative Oceanography Division, Scripps Institution of Oceanography, La Jolla, CA, USA

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: The Hudson Bay System Study (BaySys)