Protist communities along freshwater–marine transition zones in Hudson Bay (Canada)

The environmental gradient from rivers to the sea can be sharp or more gradual depending on hydrodynamic processes and the residence times within estuaries (Uncles et al., 2002; Burchard et al., 2018). Bridging the river and the sea, estuaries are important boundaries between two ecosystems and constitute major gateways delivering freshwater into the marine system. Freshwater entering the coastal ocean influences physical, biogeochemical, and biological properties and processes of the recipient marine system. Within estuaries, salinity is variable and driven by circulation patterns that depend on wind, tides, estuarine geomorphology, and the volume and flow rate of riverine water entering the estuary (Simpson et al., 1990). Due to the coupling between tidal forcing and river runoff, estuaries are complex transition areas where pronounced gradients not only of salinity but also of temperature, organic matter, light, and nutrients influence bacterial (Bouvier and del Giorgio, 2002; Crump et al., 2004) and microbial eukaryote communities (Muylaert et al., 2009; Vigil et al., 2009; Bazin et al., 2014).

Discharge timing and intensity influence the maximum turbidity zone in many estuaries (Uncles et al., 2002; Burchard et al., 2018). On regulated rivers, discharge events depend on dam-controlled runoff and affect light and nutrient availability through suspended particulate matter (Domingues et al., 2012). During high-discharge events, suspended material and nutrient loading in the estuary increase, and the maximum turbidity zone moves seaward. Conversely, low-discharge events are associated with reduced nutrient inputs, sedimentation of suspended materials in the estuarine transition zone, increasing light availability for phytoplankton, and displacement of the remaining maximum turbidity zone toward the river. On unregulated rivers in the Arctic and sub-Arctic, seasonal discharge is usually closely tied to spring freshets when snow pack melt contributes to increased flow. All of these factors will likely affect estuarine plankton communities. Because of the inaccessibility of these regions during spring, when ice conditions are unstable, Arctic and sub-Arctic systems have been studied less compared to temperate estuaries. The moderately lower salinity of the Arctic Ocean and adjacent seas and the less abrupt salinity changes between the estuary and open water could also result in less marked community transitions between fresh and marine waters. In this study, we investigated river systems in Hudson Bay (Canada), a sub-Arctic inland sea that receives up to 800 km³ year^–1 of freshwater from 23 major rivers (McClelland et al., 2012). This river runoff to the coastal Hudson Bay freshens surface waters, which can range from salinities of 25 to 29. The fresher surface waters remain on top of more saline and deeper waters with salinities up to 33 (Granskog et al., 2011).

Specifically, we compared three river systems to examine freshwater–marine coupling on biological communities shortly after ice breakup. Many larger rivers in Hudson Bay are now under the influence of hydroelectric projects that change the timing and volume of freshwater brought by rivers into the coastal sea (Déry et al., 2011). The Nelson River is the largest contributor of freshwater to Western Hudson Bay and the most impacted by hydroelectric development (Déry et al., 2016; Stadnyk et al., 2019). Since 1976, approximately 75% of the annual flow of the Churchill River has been diverted into the catchment basin of the Nelson River for enhanced hydropower production (Newbury et al., 1984; Déry et al., 2016). As a consequence, the annual freshwater discharge of the Nelson has increased and the flow regime has evolved toward a flattened hydrograph, with the high spring discharge held in reservoirs and increased winter flow to meet higher hydroelectric demand (Déry et al., 2011; Déry et al., 2016; McCullough et al., 2019). Although seasonality of the Churchill River flow was only altered slightly, the annual discharge decreased significantly since the diversion. The unregulated Great Whale River on the eastern side of Hudson Bay contributes to freshwater discharge in the same range as the Churchill River (Bhiry et al., 2011; Déry et al., 2016; Stadnyk et al., 2019) and was a counterpoint to the Churchill system.

In Hudson Bay, freshwater inflows are considered a minor source of inorganic nitrogen to the coastal sea but may contribute indirectly to the primary production by influencing large-scale estuarine circulation and input of nutrient-rich deep waters (Granskog et al., 2009; Kuzyk et al., 2010). These observations, combined with the work of Heikkilä et al. (2014) on dinoflagellate cyst distribution, support the view that coastal Hudson Bay is highly productive compared to its center. However, except for a few studies investigating the main photosynthetic protist contributors to primary production (Harvey et al., 1997; Ferland et al., 2011; Lapoussière et al., 2013), little is known about the composition and structure of microbial eukaryote communities (hence referred to as protists) in coastal Hudson Bay.

The present study aims to make a first inventory of protist communities in the coastal and estuarine environments of the Hudson Bay sequencing the V4 region of 18 S rRNA (rRNA) and the 18S rRNA gene (rDNA). Our objective was to identify protist community composition along a salinity gradient in the Nelson, Churchill, and Great Whale River estuarine systems. Samples were collected along transects at the surface and below the fresher surface layer to capture diversity associated with the estuarine processes. Within this objective, we aimed to better understand the underlying environmental mechanisms structuring the communities in Arctic estuaries and to assess the influence of river discharge on the composition and distribution of microbial eukaryote assemblages.

The study was part of a larger multidisciplinary investigation to understand the contributions of climate change and hydroelectric regulation to the freshwater–marine coupling in the Hudson Bay System (BaySys; Barber, 2014). The main aim of the ship-based study was to sample close to the seasonal peak of freshwater discharge. The original mission was to have occurred in 2017, but due to ice and logistical constraints, the ship was not able to sample Western Hudson Bay. The Great Whale River was therefore sampled opportunistically in 2017 once the ship was able to enter Hudson Bay. The original ship-based study focusing on Western Hudson Bay was subsequently carried out in 2018.

Samples were collected from three rivers that discharge into Hudson Bay and adjacent coastal sites (Figure 1). All marine and outer estuarine samples were collected from the Canadian Coast Guard Icebreaker CCGS Amundsen using a Rosette system equipped with 12-L Niskin-type bottles and a conductivity, temperature, and depth (CTD) profiler (Sea-Bird SBE-911 CTD). The rosette was also equipped with relative nitrate (In-Situ Ultraviolet Spectrometer, ISUS, Satlantic), dissolved oxygen (Seabird SBE-43), chlorophyll fluorescence (Seapoint), fluorescent-colored dissolved organic matter (fCDOM; Wetlabs ECO), photosynthetically active radiation (PAR; Biospherical Instruments QCP3200), and transmissometer (WETlabs C-Star, SeaBird Scientific) sensors. The dissolved oxygen sensor was calibrated onboard against Winkler titrations. Discrete water samples for nutrients, chlorophyll a (chl a), microscopy, and nucleic acids were collected directly from 12-L Niskin-type bottles that were closed on the upward cast.

The Great Whale River system was sampled on July 7, 2017, with river samples collected directly from the edge of the river using clean 10-L carboys that had been rinsed in 5% HCL, sterile water, and sample water. The Nelson and Churchill River systems were sampled from June 28, 2018, to July 4, 2018 (Table S1). The river water samples were collected directly into clean carboys from a zodiac or a barge deployed from the CCGS Amundsen. River conductivity and temperature were obtained using a SeaBird 19 + CTD probe. Turbidity in the Nelson and Churchill Rivers was measured with a transmissometer (WETlabs C-Star, SeaBird Scientific) deployed from a barge.

Samples for chl a, nucleic acids, flow cytometry (FCM), and particulate organic carbon (POC) along the Nelson and Churchill River transects were collected mostly from the surface using a bucket except for stations NE-D and NE-E (Figure 1), where near bottom water from 5- and 7-m depths, respectively, was collected with a submersible pump (Cyclone®) mounted on a telescope pole.

For all samples, 6 L of water was subsampled for nucleic acids. This water was sequentially filtered through a 50-µm nylon mesh, a 47-mm diameter and 3-µm polycarbonate filter, and finally through a 0.22-µm Sterivex™ unit (Millipore Canada Ltd, Mississauga, Ontario, Canada). Large (50–3 µm) and small (3–0.22 µm) fraction filters were retained for sequencing. The 3-µm filters were folded and placed in 1.5-mL tubes with RNAlater™ (ThermoFisher). The RNAlater was added to the Sterivex units, and the large and small fraction samples were stored at –80 °C until nucleic acid extraction. Nutrients, FCM, and microscopy samples were collected from the same depths and sample bottles.

Water samples for nitrate (NO₃) and nitrite (NO₂; reported as the sum NO₃ + NO₂), phosphate (PO₄), and silicate (Si(OH)₄) measurements were collected into acid-rinsed 15-mL Falcon™ tubes (Corning) after the filtration through a 0.22-µm syringe filter and analyzed on board the ship using a Bran-Luebbe 3 autoanalyzer (Grasshoff et al., 1999). All FCM samples were fixed by adding 90 µl of 25% glutaraldhehyde to 1.8 mL of seawater. Preserved samples were left at 4 °C for 30 min in the dark, deep frozen in liquid nitrogen, and stored at –80 °C until laboratory analysis. Water samples for chl a and POC were collected in dark Polycarbonate Nalgene® bottles. For chl a, samples were filtered onboard the ship onto 25-mm GF/F (Whatman) filters using a vacuum pump and placed in 10 mL of 90% acetone at 5 °C for 18–24 h to extract the pigments. Fluorescence was measured using a 10AU Field Fluorometer (Turner Designs) before and after acidification with 5% HCL (Parsons et al., 1984). Chl a concentrations were determined using the equations by Holm-Hansen et al. (1965). For POC analysis, subsamples from the Nalgene bottles and a filtration blank (from filtered seawater) for each sampling station were filtered onto precombusted (450 °C for 5 h) 25-mm GF/F filters. Afterward, filters were wrapped in tinfoil and stored at –80 °C for later analysis using a COSTECH ECS 4010 Elemental Analyser coupled with a continuous-flow DeltaPlus XP Isotope Ratio Mass Spectrometer (IRMS, Thermo Electron Co, ThermoScientific, Waltham, MA, USA) at the University of British Columbia as in Glaz et al. (2014). Microscopy samples were filtered onto black 0.8-µm polycarbonate filters, stained with 4′,6-Diamidino-2-phenylindole dihydrochloride (DAPI), mounted in immersion oil onto glass slides and stored at –20 °C until they were examined (Joli et al., 2018).

DNA and RNA samples were co-extracted from the filters using AllPrep DNA/RNA Mini kit (Qiagen, Germany) following the suggested protocol. For RNA purification, an additional DNAase step was carried out to remove possible DNA contamination. Absence of residual DNA was confirmed by Polymerase Chain Reaction (PCR) with eukaryotic primers and the cleaned RNA as template showing no amplification after 30 cycles. The RNA was then converted to cDNA using the High Capacity Reverse Transcription Kit (ThermoFisher, USA). The V4 region of 18S rRNA gene (rDNA) and 18S rRNA (rRNA) was amplified to construct libraries using a combination of universal forward E572F (CYG CGG TAA TTC CAG CTC) and reverse primers E1009 R (CRA AGA YGA TYA GAT ACC RT; Comeau et al., 2016). Amplicons were tagged for multiplexing with MiSeq®-specific linking primers, and equimolar concentrations of amplicons were sequenced on an Illumina MiSeq® by the “Plateforme d’Analyses Génomiques” (IBIS, Université Laval, Quebec, Canada). Raw paired-end reads have been merged with BBMmerge (Bushnell et al., 2017) and deposited in NCBI under the BioProject accession number PRJNA627250.

Microbial cell concentrations were measured on a BD Accuri™ C6 flow cytometer equipped with a Csampler (BD Biosciences); phytoplankton cell counts were based on chlorophyll red fluorescence and forward-scattered light (Marie et al., 2005) using 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) for 10 min with three wash and three agitation cycles between each sample. Flow rate was recalibrated daily with standard-size beads to normalize cell counts. Data were processed with BD CSampler Software. We confirmed the presence of the diatom Rhizosolenia and other cells from the DAPI-stained filters with the 40× objective and appropriate settings for high resolution using an SP8 Leica confocal laser scanning microscope. Potential cyanobacterial symbionts within the cells were searched under blue laser excitation (488 nm, Lambda scan 500–795 nm) targeting phycocyanin and phycoerythrin by using “lightning” settings on the microscope to enhance visualization.

Large (50–3 µm) and small (3–0.22 µm) fractions from 37 samples were sequenced (Table S1). Overlapping paired-end reads from the fastq files were merged using BBMerge algorithm (Bushnell et al., 2017). The data reads were quality filtered using vsearch (fastq_maxee = 0.5). Short sequences were identified in vsearch and removed from the analysis. USEARCH was used for chimera checking and Operational Taxonomic Unit (OTU) picking at a threshold of 98% similarity. The assignment of taxonomic identity was performed in mothur against the PR2 Database v4.9.0 (Guillou et al., 2013). Sequences affiliated to chloroplasts, Metazoa, and fungi were excluded from the diversity and protist proportional analysis. For information, fungi results are presented separately (Table S2). For the protist cross-comparisons between samples, small and large fraction communities were summed together, and singletons were removed from the analysis using the R package Phyloseq (McMurdie and Holmes, 2013).

To correct for the influence of rare species and differential sequencing depth, data standardization was performed with the Hellinger method using vegan (Oksanen et al., 2007) and OTUs below the threshold of 1 × 10^–5 total relative abundance were removed from the matrix table. This correction resulted in a relative abundance table of 6010 OTUs for rDNA (Table S3) and 4230 OTUs for rRNA (Table S4). A Bray–Curtis distance matrix was then calculated based on the relative abundance tables. Singletons were removed from the species table. Distance-based multivariate regression trees (db-MRT) were performed with the R package mvpart (De’Ath, 2002) to explore and predict the relationship between microbial community distance and environmental variables. This method allows partitioning of a group of samples where each successive partition from the Bray–Curtis distance predicts the influence of an environmental variable. Salinity, temperature, depth, total phytoplankton (from FCM), and nutrient concentrations (phosphate, nitrate + nitrite and silicate) were used as explanatory variables in the analysis. For each split, the method retains the variable and the associated value, which minimizes the resulting sum of within-group squared distance to the group mean for the response data. Six leaves were defined based on the minimum cross-validated relative error (CVRE) through 1,000 iterations. To test the robustness of our analysis, we compared this predictive model of communities with an unconstrained method of cluster analysis (nonmetric multidimensional scaling [NMDS]) to determine whether the clustering was similar between techniques. A Bray–Curtis distance matrix was performed on Hellinger-standardized density data based on the OTUs with the R package vegan. Environmental vectors were fitted on the ordination with the envfit function. The indicator value index IndVal was used to identify specific OTUs that characterized previously identified sample groups (Dufrêne and Legendre, 1997). Indicator value analysis was run with the R package indicspecies (v1.7.8; Cáceres and Legendre, 2009) using the multipatt function “indval.g,” and statistical significances of indicator taxa were tested by random permutations of stations (999 permutations). Only OTUs with a p value ≤ .001, indicator value (stat) ≥ 0.9, and high average relative proportions (total relative abundance ≥ 0.1%) were considered as significant, resulting in a final data set of 55 indicator OTUs. OTUs with low taxonomic assignment were further scrutinized with the “BLASTn” algorithm against the NBCI nr database. Indicator OTUs were represented on a ternary plot using the R package ggtern (Hamilton and Ferry, 2018).

To further increase the taxonomic resolution of estuarine diatom OTUs of Thalassiosira, Skeletonema, and Rhizosolenia, sequences were aligned with full-length 18S rRNA representative sequences using MUSCLE (Edgar, 2004). Representative sequences were selected based on blast results of OTU sequences on NCBI and previously published phylogenies. Reference trees were constructed with RAxML (Stamatakis, 2014) under the GTR + GAMMA model of substitution from 1,000 distinct randomized maximum parsimony starting trees (Yang, 1994). OTU queries were placed on the reference trees using the evolutionary placement algorithm in RAxML, and the expected distance between placement location value was calculated with PPlacer (Matsen et al., 2010). Final trees were visualized using the R package ggtree (Yu et al., 2017).

Within the Nelson River system, which included the Hayes River mouth (station HA-A), river waters were very warm (>20 °C) and fresh (salinity < 0.2; Table 1). The surface freshwater signal persisted and salinity increased to 20 at NE-45 (Figure 2). More typical marine salinities of >30 and colder temperatures (0 °C) were measured offshore. The deeper (32 m) marine waters at station NE-46 were even more saline (salinity: 31.7) and colder (–1.3 °C). The orientation of isopycnals supports salinity as the major parameter influencing stratification in the estuary (Figure S1). POC concentrations (Table 1) and beam transmission indicated a sharp turbidity front where the riverine water meets the seawater just before station NE-D and NE-E (Figure 3). The Churchill River water was also fresh with a low salinity of 0.1 but cooler than the Nelson River, with temperatures of approximately 10 °C at stations CH-A and CH-B. Surface salinity and temperature profiles along that transect showed a sharp freshwater gradient highlighting low riverine discharge (Figure 2). The water column was stratified at station CH-B defined by a warm (9.8 °C) and fresher (10.9) layer in the surface on top of a salty (29.4) and cold (3.2 °C) marine layer at 4 m, indicating a marine intrusion into the enclosed estuary (Table 1). Waters offshore of the Churchill River system were warmer than those offshore of the Nelson River system (Figure S1). For the Great Whale River, we recorded a salinity of 0 at the river station GW-A and 12.5 at station GW-B (temperatures not available). Surface water at station GW-C was relatively warm (9.6 °C) and less saline (25.3) compared to the more saline (31.7) and cold (–1.8 °C) deeper layer.

The highest concentrations of NO₃ + NO₂ were at the furthest upstream stations for the three transects with a maximum at the Nelson River station NE-A (0.44 µmol L^–1; Figure S2, Table 1). Concentrations were very low in the transition and coastal zones of the Great Whale and Churchill Rivers with a maximum at the surface station GW-C (0.07 µmol L^–1). Concentrations were higher in the Nelson River estuary, ranging from 0.03 to 0.43 µmol L^–1 of NO₃ + NO₂. Lowest concentrations of phosphate were recorded at all of the uppermost river stations (Figure S2) and then showed an increased seaward gradient with a maximum concentration at the surface of station GW-B (0.88 µmol L^–1). The nutrient molar ratio of N:P was <2 at all sites except for the Churchill and Great Whale Rivers, which were 7.16 and 9.54, respectively (Table 1). Highest silicate concentrations were observed in the Nelson River (38.5 µmol L^–1 at station NE-A) and decreased seaward. Silicate concentrations in the coastal waters on the Churchill and Great Whale River stations were lower with values <10 µmol L^–1 and a minimum of 0.42 µmol L^–1 at station CH-E (Figure S2).

A db-MRT and an NMDS analysis were conducted for both rDNA and rRNA OTU reads (Figure 4). Although no clear community patterns emerged based on rRNA, the db-MRT analysis based on rDNA separated six well-defined groups with low cumulative error (CVRE = 0.361, SE = 0.0893), showing that the model was a reasonable predictor of microbial community structure and associated environmental parameters (Figure 4). Salinity was the primary environmental predictor (R² = .66) of the explained variance, which first separated Groups 6 and 5 from the remaining samples. Group 6 consisted of river samples, and Group 5 tended to be more saline (S ≥ 1.3 and S < 11.35) and considered upper estuarine. Phosphate concentration was the best predictor separating Groups 5 and 6 with concentrations lowest further upstream. Salinity was also the main predictor between the lower estuarine Group 4 (S ≥ 11.5 and < 17.8) and three marine communities corresponding to the offshore sites of the three river systems (Groups 1, 2, and 3). Phosphate concentrations again were associated with the clustering of the marine stations with the Great Whale River system having the lowest values (Group 3; PO₄ < 0.33). NO₃ + NO₂ concentrations split the Western Hudson Bay marine samples into the Churchill Group 2 (NO₃ + NO₂ < 0.02) and predominantly Nelson Group 1 (NO₃ + NO₂ ≥ 0.02). The only geographic exception was the deeper sample at station CH-E, which was associated with Group 1. These results were consistent with clustering in the unconstrained NMDS analysis (Figure S3).

Relative proportions of 18S rDNA (Figure S4) and rRNA (Figures S5 and S6) reads along the transects revealed marked changes in microbial community composition from the rivers to the oceanic waters. Relative contributions and abundance of OTUs that contributed to the three salinity environments in the three systems were plotted on ternary diagrams and showed the sharp distinction of the river and sea communities in the three systems (Figure 5). The diagrams also showed a progressive replacement, in terms of percent community contribution, among shared OTUs between the estuary and sea taxa in all three systems and the river and estuary taxa in the Nelson and Churchill River systems. In addition, exclusive OTUs in the three salinity environments can be seen at the apex of the diagrams. Freshwater communities (Figure 4, Group 6) were characterized by higher proportions of rare taxa (indicated here as “others”) and unclassified eukaryotes. Although not classified as indicator taxa by our analysis, three offshore Great Whale River OTUs with a high proportion (0.8%–1.5%) of surface rDNA reads (Supplementary Data S1) were affiliated to the temperate Micromonas commoda clade (Figure S6B). The “riverine” community at all three sites was composed of typical freshwater and brackish organisms. In the Nelson River, these organisms were represented by Chrysophyceae clade C, freshwater Dinobryon spp., and freshwater representatives of the diatom genera Surirella and Navicula (Figure 6). The marine communities (Groups 1–3) were composed of higher proportions of Dinophyceae, Choanoflagellatea, Prymnesiophyceae, Telonemia, and Mamiellophyceae. Estuarine communities (Groups 4 and 5) both showed an increase in proportions of reads affiliated with heterotrophic taxa including Cercozoa and Katablepharidaceae (Figures S4–S6). Ciliates showed an increasing proportion of reads in the estuarine transition zone of all three river systems, but the composition of the estuarine ciliate communities differed among the estuaries (Figure S7A). In the Churchill River estuary, herbivorous ciliates Urotricha spp. and Didiniidae spp. increased in relative proportions while Mesodinium spp. and Strombidium spp. were recovered from the Great Whale River estuary. The most abundant OTU in the Great Whale River estuary was related to Mesodinium rubrum (Myrionecta rubra). The Nelson River estuarine groups included more diatoms with genera classified as Thalassiosira and Skeletonema. The Thalassiosira OTUs mapped to a node consisting of Thalassiosira pseudonana and Thalassiosira weissflogii (Figure S8), and the Skeletonema OTUs mapped to Skeletonema potamos (Figure S9).

To visualize the indicators of the three offshore groups from Figure 4, we plotted the relative proportion of indicator species (Figure 7; Table S5) based on the three groups (Groups 1–3) that segregated the three regions. This visualization highlighted species that were present in a single region and showed species present in two regions to the exclusion of a third. Churchill offshore indicator taxa (Group 2) included two representatives of the Marine Ochrophyta group MOCH-2 (Massana et al., 2014; OTUs 2730 and 745), several Cercozoan OTUs including Protaspa spp. (OTU 1771), the Bolidophyceae classified as Triparma strigata (OTU 1205), the haptophyte Chrysochromulina sp. (OTU 327), and the choanoflagellate Lagenoeca antarctica (OTU 1323). Offshore the Great Whale River (Group 3), which included both surface and Subsurface Chlorophyll Mamimum (SCM) offshore data, indicator OTUs were classified as Syndiniales (OTU 3280 and 913), and an OTU affiliated with the marine Stramenopiles MAST_4D (OTU 587; Figure 7).

The Nelson offshore indicators (Group 1) were Rhizosolenia spp. and an OTU matching Strombidiida B (PR2 classification). We examined the Nelson offshore indicators in more detail as more samples had been taken along that marine transect. The diatoms were placed on the Rhizosoleniales reference tree with most at a node of Rhizosolenia spp. An unclassified diatom (OTU 184) was placed on the lower branches of a reference phylogenetic tree (Figure S10). These Rhizosoleniales were highly specialized to the Nelson River system (at the apex of the ternary plot; Figure 7).

The Churchill (Group 2) and Nelson River (Group 1) coastal systems shared a number of taxa that were absent from the Great Whale River coastal region. The Churchill River coastal system also shared taxa with the offshore Great Whale River system (Group 3). No indicator taxa were shared between the Great Whale and the Nelson Rivers. The cryptophyte Teleaulax gracilis (OTU 2014) was the main photosynthetic OTU in the Western Hudson Bay (Nelson and Churchill River coastal region). More strikingly was the number of heterotrophic flagellates shared between the two sites, Ebriida sp. (OTU 383) and Minorisa sp. (OTU 199), Telonema (OTU 222, 555, 1098, 2554 and 818), Cryothecomonads (Protaspa; OTU 712 and 543), the choanoflagellate Stephanoecidae (OTU 896), and a thecofilosean in the Ventricleftida (OTU 1239). The shared taxa between the Churchill and Great Whale coastal region included different heterotrophic flagellates: MAST_1A (OTU 605), Picozoa—a centroheliozoan Pterocystida sp. (OTU 1554), and an environmental clade of Dinoflagellates (OTU 86, 21806 and 2032). In addition, there were several diatoms including Chaetoceros neogracile (OTU 1928) and Chaetoceros decipiens (OTU 7975) common to the Churchill (Group 2) and Great Whale (Group 3) Rivers.

The diatom genera in surface coastal waters were confirmed by confocal microscopy by examining the DAPI-stained microscopy filters under a range of excitation–emission settings. Major diatom taxa included Halochaete Chaetoceros spp. and Rhizosolenia. As with the indicator analysis, Rhizosolenia was particularly common in the Nelson River offshore samples (Figure S11a). Although we could clearly distinguish chloroplasts, no cyanobacteria with phycocyanin or phycoerythrin fluorescent signals were detected inside of the frustules, suggesting the absence of the potential N-fixing symbionts within the Rhizosolenia cells (Figure S11b).

Salinity was the major environmental parameter separating microbial communities along the three transects (Figure 4). The Nelson and Churchill River systems have very different freshwater and seawater mixing regimes. At the time of sampling, the Nelson River estuary had warm and fresh riverine water grading into cold and salty marine water. The gradual dissipation of freshwater and cooler surface waters suggested that the Nelson River plume extended up to 25 km offshore from the river mouth (Figure 2). The freshwater discharge rate, underlying geomorphology (a channel), and estuarine circulation led to the development of a strong maximum turbidity zone at the convergence of marine and riverine waters (Figures 2 and 3). In contrast, for the Churchill estuary, the salinity of the deeper layer was more uniform with open ocean salinity near the bottom of the sampled river region. This saline water below the fresh riverine water caused strong vertical stratification at the CH-B station with little evidence of freshwater offshore (Figures 2 and S1). The different salinity profiles reflected the strength of the river runoff, which differs greatly between these two Western Hudson Bay estuaries and would have been amplified by the ongoing Churchill River diversion into the Nelson catchment beginning in 1974 (Newbury et al., 1984; Déry et al., 2016).

The influence of salinity on the microbial communities was not surprising and suggests that the sub-Arctic communities in summer are driven by similar salinity constraints as reported in temperate estuaries for protist communities (Muylaert et al., 2009; Bazin et al., 2014; Lee et al., 2017; Filker et al., 2019). The microbial communities in the three river systems all contained characteristic freshwater “riverine” taxa. The proportion of the reads from freshwater organisms abruptly declined with increasing salinity, and freshwater eukaryotes were essentially absent in marine waters, indicating a clear dispersal barrier. Both the upper and lower estuarine communities contained a mix of freshwater, marine, and presumably euryhaline or lower-salinity (brackish) taxa. Taxonomic differences between marine and freshwater groups have been observed in bacteria, archaea, and unicellular eukaryotes consistent with evolutionary separation between marine and freshwater lineages (reviewed by Logares et al., 2009). Among parameters that could prevent local adaptation of riverine organisms to high-salinity conditions is a slower cellular division rate as a cost of osmoregulation (Demmig-Adams et al., 2017) exposing taxa to competition with locally adapted lineages and net losses by predation. On the other hand, marine or ice-associated species could be selected against in the warmer surface waters from the stratified estuaries as Arctic species tend to be temperature sensitive (Lovejoy et al., 2007; Daugbjerg et al., 2018). Thus, although the temperature was not considered as the most influential parameter separating estuarine and marine groups, the estuary could also represent a possible thermal barrier for the observed marine species as the temperature in the shallower river water tended to covary with salinity.

The three coastal communities were composed exclusively of marine taxa (Figure 5), but the species composition differed. Based on the Bray–Curtis distance, the offshore marine coastal communities were separated into distinct Nelson, Churchill, and Great Whale River clusters (Groups 1, 2, and 3, respectively; Figure 4). The only exception was the deeper water community at CH-E, which clustered predominantly with the Nelson River Group 1. Phosphate and NO₃ + NO₂ concentrations were the main factors separating the three coastal locations but accounted for <4% of the total explained variance from the db-MRT conducted using the whole data set. Similar to other regions of the Arctic Ocean in summer (Blais et al., 2012; Martin et al., 2013), low N:P ratios (NO₃ + NO₂:PO₄) were consistent with primary production limited by the availability of inorganic nitrogen (Table 1). Arctic nutrient concentrations often closely follow the development of the spring bloom, with nitrate concentrations drawn down in surface waters as the bloom progresses (Mei et al., 2003). In Western Hudson Bay in June 2018, an ice-edge bloom followed ice breakup in mid-May depleting nutrients in the surface mixed layer (Barbeido de Freitas et al., n.d.; Matthes et al., n.d,). The accepted phenology of Arctic phytoplankton is that assemblages move from a peak of diatom or Phaeocystis colonies, depending on the availability of silica, to small single-celled flagellates that rely on recycled nitrogen sources such as ammonium and urea (Mei et al., 2003; Simpson et al., 2008; Joli et al., 2017). At the time of sampling, NO₃ + NO₂ and silicate concentrations were 10-fold higher in the Nelson offshore samples compared to the Churchill samples (0.1 compared to 0.01). The lower nitrate concentrations in the offshore Churchill River were consistent with the opening of the northwestern Hudson Bay coastal polynya and a precocious bloom developing at the ice edge of the polynya. The bloom was evident in the Nelson River estuary a few weeks later when Rhizosolenia was found at the seaward station NE-46 (Figure 6).

Rhizosolenia spp. are found in both marine and brackish waters and can form huge blooms in temperate estuaries and coastal waters (Rousseau et al., 2002). They have been consistently reported in Arctic Seas over the past 100 years (Lebour, 1930), and some species have a pan-Arctic distribution (Lovejoy et al., 2002; Poulin et al., 2011). Although the molar N:P ratio was low, the relatively high proportions of Rhizosolenia in both rDNA and rRNA reads in the offshore stations of the Churchill and Nelson Rivers suggest that Rhizosolenia cells were active. Some Rhizosolenia spp. can bypass nitrogen limitation by forming an endosymbiotic association with the nitrogen-fixing cyanobacteria Richelia spp. (Villareal, 1990; Singler and Villareal, 2005; Zeev et al., 2008). Although we could clearly distinguish chloroplasts, we were not able to detect cyanobacterial symbionts within the Rhizosolenia cells using a confocal microscope (Figure S11). At the time of the Nelson River sampling, the limit of the mobile ice cover (defined as 50% ice concentration) was approximately 90 km offshore of the estuary and a large phytoplankton bloom was observed in the open water at the ice edge (LC Matthes, personal communication). Rhizosolenia blooms are also reported to be associated with ice-edge conditions and have been used to trace ice conditions and carbon flow through the polar marine trophic web (Goutte et al., 2013; Brown et al., 2018), and Rhizosolenia may have been part of an ice-edge bloom. The healthy chloroplasts seen under confocal microscopy and the presence of Rhizosolenia in the rRNA libraries suggest that light and nutrient levels continued to be favorable closer inshore where they may have arrived initially through estuarine circulation.

Most of the marine specialist taxa present in the very N-depleted offshore Churchill River water were heterotrophs (Figure 7). These included OTUs affiliated with Protaspa, a thecofilosean (Cerozoa; Howe et al., 2011). Related species have been reported as algal predators or parasites that occur after the spring diatom bloom in San Pedro Channel, CA (reported as Protapsis; Berdjeb et al., 2018), which would be consistent with the Churchill River system at the end of a bloom. In Arctic marine waters, a sister taxon in the same family, Cryothecomonas, was associated with waters affected by ice melt (Thaler and Lovejoy, 2012). Another identified indicator species was an environmental Marine Ochrophyta MOCH-2, which is one of the five MOCH groups described to date and most likely feeds on bacteria (Massana et al., 2014), as would the OTU with the closest match to the choanoflagellate L. antarctica, originally isolated from deep Antarctic waters (Nitsche et al., 2007). Two small photosynthetic protists were also recorded as indicators of Churchill River coastal waters: the Prymnesiophyceae Chyrsochromulina sp., a mixotrophic taxon reported in highly stratified or ice melt–freshened waters in the Arctic (Comeau et al., 2013; Lovejoy, 2014; Ardyna et al., 2017), and a Parmales (bolidophyte), with match closest to T. strigata (99.5%) and environmental sequences including those from the Arctic (Terrado et al., 2013). Non-silicious bolidophytes are reported to be mixotrophic, consuming bacteria (Piwosz et al., 2013), but whether the OTU represented a flagellate or silicified form cannot be stated without confirmation by electron microscopy (Yamada et al., 2020). In summary, following ice breakup about 1 month earlier (beginning of June 2018), most of the indicator species off the Churchill estuary were associated with postbloom conditions with increasing surface temperature, remaining mobile sea ice and low inorganic nutrients concentrations in the area at the time of sampling (July 4, 2018).

Two of the indicator taxa from the marine coastal Great Whale River (Figure 7) belonged to the marine parasitic order Syndiniales (Guillou et al., 2008). Because these parasites tend to infect dinoflagellates and other heterotrophs that occur later during the postbloom stage, their predominance is consistent with the offshore Great Whale River system being at a later bloom state on July 7 and supported by the earlier opening of Eastern Hudson Bay in 2017 relative to 2018 (Kirillov et al., 2020). The presence of M. rubrum along with fewer cryptophytes would also be consistent with being later in the bloom stage. A third Great Whale River marine OTU was placed within the uncultured marine stramenopile clade MAST_4. Members of MAST_4 are picoplanktonic protists widely reported in temperate waters but rare in Arctic samples (Massana et al., 2004, 2014). The clade was previously reported from western Arctic surface waters, however, during the then record warm summer in 2005 (Lovejoy and Potvin, 2011), which suggests that warmer water species may have been present in Hudson Bay at the time of sampling. In support of this notion was the high proportion of surface rDNA reads affiliated with the temperate M. commoda clade (Figure S6B). The M. commoda clade co-occurred with Micromonas polaris. M. polaris has a pan-Arctic distribution and has been reported throughout the year (Lovejoy et al., 2007; Joli et al., 2017). The presence of the more temperate taxa in the relatively warm surface water (9 °C) at GW-C suggests that warming of surface water favors non-Arctic populations of coastal protists in Hudson Bay, but as of yet, Arctic species are not excluded.

A distinct estuarine community was detected in the transition zones of all three river systems (Figure 5). The communities were a mix of specialist taxa, freshwater and marine species. The marine diatom species, mostly Chaetoceros spp. found in the Churchill and Great Whale River transition zones, tended to decrease toward the shore consistent with advection and loss of taxa shoreward. In the Nelson River transition zone, several diatom OTUs with closest matches to the euryhaline species S. potamos, T. pseudonana, and T. weissflogii (Hevia-Orube et al., 2016) were found at the intermediate salinities in both rRNA and rDNA libraries (Figures 6, S8, and S9). The same species contribute to phytoplankton blooms in temperate coastal environments (Cloern et al., 1985; Urrutxurtu et al., 2003; Quinlan and Phlips, 2007; O’Boyle and Silke, 2010). In particular, S. potamos is a common species in temperate rivers and lakes of Europe and North America but can grow at salinities up to 24 (Paasche, 1975; Torgan et al., 2009). The spatial distribution of these brackish diatoms in the Nelson River estuary coincided with silicate-rich waters and the maximum of POC concentrations measured at the surface (Figure 6).

Heterotrophic flagellates and ciliates have been reported from turbidity zones of major rivers flowing into the sea, for example, the Saint Lawrence (Lovejoy et al., 1993) and the Columbia River (Herfort et al., 2011). In Hudson Bay, we found distinct communities of heterotrophic protists in all three estuarine transition zones, consistent with a distinct euryhaline habitat. In the Nelson River estuary, high POC concentrations and low beam transmission signals (Figures 3 and 6) in the transition zone were indicative of a maximum turbidity zone. Both rDNA and rRNA Nelson River estuary data sets indicated higher proportion of phagotrophic thecofilosean Cercozoa including lineages related to Mataza sp. and Ventricleftida sp. (Figure S5). Many of these flagellated amoebae are described from soil and freshwater habitats but are not restricted to these environments. For example, the first described Mataza was isolated from a surface seawater sample in Tokyo Bay (Yabuki and Ishida, 2011), and related lineages have also been detected in both fresh and brackish waters in the estuary of the Vistula River in the Baltic Sea (Piwosz et al., 2018). Katablepharidales, another phagotrophic flagellate group, were also found in the brackish environment of the Nelson River, including 13 OTUs closely related to Katablepharis japonica, which is a broad-spectrum predator, feeding on heterotrophic bacteria, diverse phytoplankton, dinoflagellates, and ciliates (Kwon et al., 2017). The relatively high proportions of Katablepharis and other phagotrophs in the Nelson estuary suggest that local conditions in the estuarine environment favor the development of an autochthonous heterotrophic community in the estuarine transition zone (Figures S4 and S5). Moreover, the peak of rDNA Katablepharis reads at station NE-C coincided with the position of the turbidity front (Figures 3 and S4). A similar ecologically constrained distribution of Katablepharis was previously reported in the estuarine turbidity maximum zone of the Columbia River (Herfort et al., 2011; Kahn et al., 2014). Herfort et al. (2011) hypothesized that Katablepharis cells may be retained as cysts in the turbidity maximum; however, the co-occurrence with other bacterivores (Telonemia, Mataza, and Ventricleftida) in our rRNA data would be consistent with an active in situ grazer community (Figure S5).

Ciliates accounted for a greater proportion of reads in the estuarine transition zone of all three river systems, but with species differences found in the three systems. Ciliates are a significant component of the Arctic mixotrophic protist communities (Stoecker and Lavrentyev, 2018), and in the Great Whale River estuarine region, the most abundant OTU was related to M. rubrum (M. rubra). This ciliate, which is part of a species complex (Lasek-Nesselquist and Johnson, 2019), is routinely reported from temperate estuarine communities and has also been reported to form blooms at the ice–water interface during pack ice formation in the Arctic Ocean (Olsen et al., 2019). M. rubrum is dependent on sequestered chloroplasts from cryptophyte prey such as Teleaulax and Plagioselmis, which were also found in the same Great Whale River samples. Teleaulax itself is favored by low light conditions (Johnson et al., 2006), suggesting that there may have been a turbidity zone at this site; unfortunately, we were not able to measure CDOM or POC at the time. The Nelson River brackish conditions in the estuary tended to favor ciliates Tintinnidium spp., Didiniidae spp., and Litostomatea spp. These genera were previously reported from brackish environments, in the upper part of the transition zone with Strombidium spp. and in the lower part, as previously detected in other Arctic marine samples (Onda et al., 2017). The spatial repartition of the brackish community in the estuaries appears to be species-specific and suggests multiple ecological niches in the estuarine transition zones, especially in the Nelson River where two distinct estuarine communities were identified (Figures 6, S4, and S5).

The definition of the estuary as a boundary between freshwater and marine systems has been the subject of debate, with two types of transition zones proposed: ecocline and ecotone (Kent et al., 1997; Attrill and Rundle, 2002). For most of the protist groups, the patterns of community assemblage along the gradient resembled progressive change from fresh to brackish and from brackish to marine, with few species occurring across the whole transition (Figure 5). A gradient in protist assemblages from the river to the marine systems would support the view of the estuarine ecocline. However, we distinguished a distinct subset of estuarine species in all three systems, with the Nelson River, which was also the most intensively sampled, having the highest representation of estuarine taxa (Figure 5). These species distributions fit with the definition of an ecotone, which is a narrow ecological zone between two different and relatively homogeneous ecosystems. This community was composed of few species that may have an inherent ecological tolerance and resilience to stress induced by the high environmental variability (Elliott and Whitfield, 2011). From our observations, the definition of the estuarine boundary as an ecotone or an ecocline was variable, depending on the taxa studied, highlighting the estuarine transition zone as a combination of the two ecological concepts rather than a clear distinction between the two.

The formation of a transition zone in an estuary is the result of complex freshwater and marine interactions that are highly variable in time and space and differ substantially among estuaries (Uncles et al., 2002; Burchard et al., 2018). In temperate and sub-Arctic seas, estuarine fronts observed at the convergence of river runoff and tidal forcing can trap suspended sediments and planktonic organisms in a maximum turbidity zone, which is a common feature in estuaries (Frenette et al., 1995; Hetland and Hsu, 2013). Within this “hydrodynamic trapping” environment, conditions of high turbidity, low light, high organic matter accumulation, and high microbial activity should be advantageous for heterotrophic and mixotrophic plankton. In the Nelson River, the estuarine circulation led to the formation of a turbidity front between station NE-C and NE-E that coincided with the peaks observed for several estuarine species (Figures 7, S4, and S5). Heterotrophic flagellates such as Katablepharis or Cercozoa and ciliates were retained there and would sustain an active microbial food web by influencing bacterial dynamics, carbon remineralization, and nutrient regeneration through phagotrophy and grazing. Diatoms were also recovered in the Nelson River turbidity maximum zone, highlighting that some taxa were able to grow and adapt to low light conditions when nutrients, especially silicate, are available. In the Churchill and the Great Whale systems, the estuarine community was constrained in a narrower transition area and was mainly composed of mixotrophic ciliates. This difference in the composition and distribution of the estuarine community may be related to geomorphology and lower river discharge of the Churchill and Great Whale Rivers compared to the Nelson River. However, for the Great Whale River sampled in 2017, interannual variability cannot be ruled out since the river runoff and the start of the productive season might differ. Preliminary hydrodynamic modeling in the Nelson River estuary has shown that water could reside in the transition zone more than 30 days before flushing out (K Wong and K Sydor, Manitoba Hydro, personal communication). Such a period would be long enough to allow the development of an adapted autochthonous community that may move within the turbidity front depending on river runoff and tidal cycle.

The question remains as to how protist assemblages in coastal Hudson Bay might change in response to different river flow regimes, whether as a result of climate change or managed discharge. In temperate estuaries, ecological disturbance associated with freshwater regulation has been described in the Guadiana estuary, where abundance and biomass of diatoms and cyanobacteria decreased (Domingues et al., 2014), and in the Danube River, where a shift from a diatom to a flagellate-based community was attributed to a reduction of silica discharge (Humborg et al., 1997, Rocha et al., 2002). What we have shown here is that the three estuaries had different species assemblages, with the Nelson River system being more mixotroph and heterotroph-dominated in keeping with the presence of a maximum turbidity zone. The extent, location, and seasonal duration of this zone would be influenced by the timing of discharge and flow rate. Global climate models predict changes in precipitation that could increase annual discharge into Hudson Bay (Steiner et al., 2015); higher flow rate through the upstream reservoirs and increased permafrost melt will add uncertainty to current models (Hinzman et al., 2020). Increased flow and melting permafrost will also affect the input of labile organic material into the Nelson system (Wrona et al., 2016). Greater discharge could move the maximum turbidity zone offshore into deeper waters and limit the ability of the estuary to act as an organic carbon remineralization zone over the spring and summer. Organic nutrients could then be pushed offshore and promote harmful algal blooms that are often linked to high inputs of labile organic matter (Anderson et al., 2002). Already in other regions of the Arctic, there are reports of such species that may harm Arctic marine mammal populations (Joli et al., 2018). Historically, the high flow would have come in spring, but the managed Nelson system has dampened the seasonal signal with a flow rate that is conducive to maintaining a maximum turbidity zone, perhaps mitigating the potential for harmful blooms. Ongoing studies are needed to verify the relationship between species assemblages and changes in hydrography. Overall, the potential consequences for water quality and higher food webs indicate a need to monitor protist communities in Hudson Bay and evaluate potential harmful effects due to ongoing anthropogenic change.

This study provides the first inventory of estuarine and coastal microbial eukaryotes in Hudson Bay using molecular techniques. We found that estuarine circulation was a major driver of the dynamics and composition of microbial eukaryote communities in these sub-Arctic estuaries, leading to the formation of ecological niches. Due to the fast response of microbes to environmental change and the complexity of the estuarine circulation processes, direct consequences of river regulation on the microbial communities of the Nelson and Churchill Rivers are difficult to predict and would require temporal monitoring on both regulated and unregulated rivers. This study gives a first glimpse of the distribution and composition of estuarine and coastal protists of Hudson Bay, providing a baseline reference for understanding potential shifts in the microbial community composition in a changing environment.

Raw paired-end reads have been deposited on NCBI under the BioProject accession number PRJNA627250 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA627250). Conductivity, temperature, and depth/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:

Tables S1, S2, S5. Figures S1–S11. Supplemental materials for the main text. File: Supplementary_material.docx.

Table S3. Relative abundance table of 6010 rDNA OTUs, where each sequence has been classified according to the PR2 database (v4.9.0). File: Supplementary_Table_S3.DNA_OTUs_table_relative6010.xlsx.

Table S4. Relative abundance table of 4230 rRNA OTUs, where each sequence has been classified according to the PR2 database (v4.9.0). File: Supplementary_Table_S4.RNA_OTUs_table_relative4230.xlsx.

The authors thank P. Guillot for processing CTD data and G. Deslongchamps and J. Gagnon for collecting and processing nutrient data. Sequencing was from the Plateforme d’Analyses Genomiques (IBIS, Université Laval and support for data analysis from Compute Canada. The CTD data were collected on board 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. They are also grateful to the crew of the CCGS Amundsen for logistics and Marianne Potvin and Carlee Morency for laboratory work.

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 discovery and northern supplement grants to CL. The authors also acknowledge support from ArcticNet a Network of Centres of Excellence (project Arc3Bio) and from the Fonds de recherche du Québec Nature et Technologies (FRQNT) to Québec Océan. LJ received mobility grants from Sentinelle Nord (Canada First Research Excellence Fund) and Québec Océan. Conductivity, temperature, and depth data were made available through the Amundsen Science program, supported by the Canada Foundation for Innovative and the NSERC. The analyses were conducted using Compute Canada facilities.

The authors declare that they have no conflict of interest.

Contributed to the conception and design: LJ, DK, CL.

Contributed to acquisition of data: LJ, DK, LCM, JÉT, CJM.

Contributed to analysis and interpretation of data: LJ, DK, LCM, AV, JÉT, CJM, CL.

Commented and/or revised the article: All authors.

Wrote the manuscript: LJ and CL with comments by all authors.

Approved the submitted version for publication: All authors.

How to cite this article: Jacquemot, L, Kalenitchenko, D, Matthes, LC, Vigneron, A, Mundy, CJ, Tremblay, J-E, Lovejoy, C. 2021. Protist communities along freshwater–marine transition zones in Hudson Bay (Canada). Elementa: Science of Anthropocene 9(1). DOI: https://doi.org/10.1525/elementa.2021.00111.

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

Associate Editor: Jeff S. Bowman, Scripps Institution of Oceanography, University of California San Diego, La Jolla, CA, USA

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: BaySys