Biogeography and phenology of phytoplankton and other planktonic microbial eukaryote species are intimately associated with local nutrient and light conditions determined by mixing and stratification dynamics of the sunlit water column. In recent decades, global warming has fundamentally altered oceanic temperature and salinity, resulting in increasing stratification of the world upper ocean (Li et al., 2020). By inhibiting vertical mixing of the water column, stratification slows the exchange of nutrients and gases, which ultimately decreases marine productivity (Matear and Hirst, 2003; Behrenfeld et al., 2006; Breitburg et al., 2018). While low latitude oceans are mainly stratified by temperature, for Arctic and sub-Arctic seas, by far the most significant determinant of stratification is salinity (Carmack, 2007). Upper water column stratification is particularly pronounced in high latitude regions where rivers discharge, freshening the surface waters and increasing stratification. The strong water column stratification prevents vertical transfer of nutrient-rich deeper water into the euphotic zone (Arrigo et al., 1999; Tremblay and Gagnon, 2009; Ferland et al., 2011). Because of anthropogenic warming, river discharge is increasing along with an accelerated hydrological cycle, impacting Arctic regions including the large Arctic-sub-Arctic Hudson Bay Complex (Stadnyk et al., 2020). An additional concern is that warming conditions could also increase the probability of harmful algal events with impacts on human health and local food chains (Anderson et al., 2021).

Hudson Bay (HB) is the core of the Hudson Bay Complex that connects Arctic modified Pacific origin waters to the Labrador Sea through Hudson Strait (HS) and Ungava Bay (UB). The Hudson Bay catchment covers one third of Canada’s land mass and is dominated by large riverine inputs (Déry and Wood, 2005; Macdonald and Kuzyk, 2011). In the greater Hudson Bay a low nutrient, warmer and fresher summer mixed layer sits atop colder saltier waters that are more nutrient rich (Granskog et al., 2011). This surface layer ranges from 30-60 m depth over most of the HB and contains most of the freshwater load from the annual accumulated river runoff of the catchment as well as seasonal sea-ice melt. These large volumes of freshwater feed an anti-clockwise north flowing boundary current along the eastern coastal region of HB (Granskog et al., 2009). The warmer, less-saline and strongly stratified waters limit phytoplankton production of HB in late summer (Harvey et al., 1997; Ferland et al., 2011).

HS+UB regions are influenced by different hydrodynamic processes compared to HB. The HS+UB region drains a much smaller catchment basin and is characterized by strong tidal and meteorological mixing that weakens stratification and increases biological productivity (Drinkwater and Jones, 1987; Straneo and Saucier, 2008). However, although net annual productivity differences between HB and HS+UB are reported, other than historic species lists of mostly phytoplankton from microscopy (Poulin et al., 2011), little is known about the composition of other microbial eukaryotic communities in the two hydrologically contrasting systems. To our knowledge, there have been no direct comparisons of the summer protist communities that occur when stratification is greatest in HB, but nearly absent in the HS+UB region. Here, we addressed this data gap using molecular techniques to identify the protist taxa and in combination with physical oceanographic data sought to gain understanding for the role of stratification in species selection in this high latitude region. We focused our analysis on protist communities in the 3 µm to 50 µm size range, which includes nanoplankton (2-20 µm) and the slightly larger diatoms, dinoflagellates and ciliates, which are prey for zooplankton and commonly reported using light microscopy. As this size fraction would be biased towards taxa that are reported from microscopy, our results should facilitate transition to community-based monitoring by local indigenous communities that live along the coastal Ungava Bay Peninsula. However, as filtration is not perfect, typically some smaller and larger taxa are also retained on the filters, due to filter blockage and cell breakage, respectively, material captured on the filters would provide catalog of taxa suitable for comparison between different systems.

To address which organisms were likely more active as opposed those that persist under the different conditions (Morency et al., 2022), we used amplicon sequencing targeting both 18S rRNA genes (rDNA) and its transcriptional products (rRNA). Although rRNA cannot be directly used to assess the growth and activity of microbes, the variability of rRNA within a species may provide information on the potential of organisms to synthesize proteins in a particular environment (Blazewicz et al., 2013). In addition, given that harmful algal groups are prone to bloom under summer conditions, and climate change scenarios are of concern to local communities (McKenzie et al., 2020), we paid particular attention to the potential presence of harmful algal species in eastern coastal HB and HS+UB under their respective summer stratification regimes.

The study was carried out in eastern Hudson Bay, the southern reaches of Hudson Strait and along the coast of Ungava Bay (Nunavik, Québec, Canada) ( Table 1 ). Sampling was realized onboard the Canadian research icebreaker CCGS Amundsen during two cruises in early July 2017 and 2018. Seawater was collected directly from 12-L Niskin-type bottles mounted on a rosette system equipped with a relative nitrate (In-Situ Ultraviolet Spectrophotometer, ISUS, Satlantic), dissolved oxygen (Seabird SBE-43), chlorophyll fluorescence (Seapoint) and a fluorescent colored dissolved organic matter sensor (fCDOM; Wetlabs ECO). We targeted the surface layer (1 or 2 m) and the subsurface chlorophyll maximum (SCM). The depth of the subsurface chlorophyll maxima (SCM) was determined on the downward cast. For each depth, 6 L of seawater was first prefiltered through a 50 µm mesh to minimize microzooplankton and then filtered onto a 3 µm pore size polycarbonate filter (Millipore, Canada) so the nominal community collected would be in the 3 to 50 µm range. The polycarbonate filters were placed in 2 ml microtubes containing 1.8 ml of RNAlater (Thermo Fisher Scientific, USA). After 30 min at room temperature, filters were transferred to –80°C until nucleic acid extraction. Water samples for nitrate (NO₃ ^-), nitrite (NO₂ ^-), phosphate (PO₄ ^3-) and silicate (Si(OH)₄) measurements were collected into acid rinsed, then sample rinsed, 15-ml Falcon™ tubes (Corning, USA) after filtration through a 0.22 µm syringe filter. The nutrient samples were analyzed on board the ship using a Bran-Luebbe 3 autoanalyzer (Grasshoff et al., 1999). The stratification parameter was defined with the Brunt-Väisalä (or buoyancy) maximum frequency N² calculated for each station in Ocean Data View (v 5.2.1; Schlitzer, 2002). The more stratified the water column, the higher the N₂ value.

DNA and RNA were co-extracted from the filters using AllPrep DNA/RNA Mini kit (Qiagen, Germany) following the manufacturer protocol and V4 region of 18S rRNA gene (rDNA) and 18S rRNA (rRNA) was amplified and sequenced using a combination of universal forward E572F (5’-CYG CGG TAA TTC CAG CTC-3’) and reverse primers E1009R (5’-CRA AGA YGA TYA GAT ACC RT-3’) using the cycling conditions as described in Comeau et al. (2011). 26 Samples from 11 sampling events were sequenced at the Institut de Biologie Intégrative et des Systèmes, Université Laval, Sequencing Plateform ( Table 1 ). Overlapping paired end reads from the fastq files were processed within the qiime2 environment (Caporaso et al., 2010) using dada2 v2019.10 (Callahan et al., 2016). The quality of the forward and reverse reads was estimated using fastqc v0.12.0 (Andrews, 2010). For each run, forward and reverse reads were demultiplexed. Removing of primers, trimming of forward and reverse reads, denoising of low-quality reads, merging and removing of chimeras was performed using the denoise-paired command in dada2 (trim-left-f = 5, trim-left-r = 5, trunc-len-f = 285, trunc-len-r = 260, maxEE = 1). The two denoised runs were merged and taxonomy was assigned to each ASV in mothur using the PR2 database v4.12 (Guillou et al., 2013). Sequences affiliated to Metazoa and chloroplasts were removed from the analysis as well as sequences at the unclassified Phylum level using the R package Phyloseq v1.32.0 (McMurdie and Holmes, 2013). For dinoflagellate taxa detected as potentially harmful, taxonomy was refined by BLASTn against the NCBI nr database ( Table S2 ).

To correct for differential sequencing depth of the samples, data were transformed to relative abundance and ASVs below the threshold of 1 × 10^-5 total relative abundance were removed from the resulting matrix table. The rDNA and rRNA ASVs and reads were analysed on the same pipeline in Phyloseq, resulting in a relative abundance table of 856 ASVs. All the subsequent clustering analysis were run in the R environment. A square-root transformation was applied to the relative abundance table and a Bray-Curtis (BC) matrix was calculated. Hierarchical clustering was performed using the “ward.D2” method on the BC matrix. Distance-based redundancy analysis (db-RDA) of microbial communities was performed with standardized environmental variables using the capscale() function in vegan v2.5.6. We adjusted the R² results using Ezekiel’s formula with the RsquareAdj() function. The adjusted R² measures the unbiased amount of explained variation and was used to reduce the number of significant variables using a forward selection and 999 ANOVA permutations. Co-varying variables were removed from the final db-RDA calculation. Pearson correlation between environmental variables and relative abundance of taxa at the “Genus” classification level was calculated with the stats R package and visualized using ggcorrplot v0.1.4 (Kassambara, 2019). Taxa with significant correlation (P < 0.001) with at least one environmental parameter in the rDNA and rRNA dataset were retained for specific correlation analysis. Linear regression models were used on these selected taxa to further test the influence of environmental parameters on the distribution of rDNA and rRNA reads. To assess the distribution of the most abundant taxa, the relative abundance of ASVs were summed at the species level and zscore-transformed ((x-mean) sd^-1). When no known species were identified for the group, taxa were annotated at the highest taxonomic level. DNA and RNA sourced datasets were treated separately. To assess significant differences in mean relative abundance between the HB and HS+UB clusters, Student’s t-tests were performed with a P-value threshold of 0.05.

Marine microbial eukaryotes respond rapidly to changing environmental conditions and can have a direct influence on the energy and carbon transfer to higher trophic levels (Worden et al., 2015). Here we showed that the structure of protist assemblages recovered from 18S rDNA and rRNA distinguished the surface and SCM as two distinct communities in eastern HB, with a single more uniform HS+UB community. The differences in species composition between the eastern HB compared to northern and western stations surrounding the HS+UB were consistent with hydrological conditions in the two marine regions ( Figures 1 , 2 ).

For the stations that were sampled during both years, rRNA and rDNA showed little or no differences in the microbial community structure, suggesting our observations were characteristic of summer conditions in this part of the HB complex. Surface temperatures at st688 were the highest in 2017 in HS+UB, which may have influenced the clustering of the rRNA sample in the HB group. Overall, the strength of the stratification was a key driver of microbial eukaryotic community structure and distribution, especially in HB ( Table 1 ). The communities in eastern HB were vertically separated in keeping with greater stratification associated with a fresher warmer summer surface mixed layer above a well-established SCM ( Figures 1 , 2 ). The development of an SCM dominated by diatoms in summer is a characteristic feature of the HB (Roff and Legendre, 1986; Ferland et al., 2011; Matthes et al., 2021). The SCM was particularly well-defined in 2017 and was mainly associated with Thalassiosira spp. and Chaetoceros spp., two genera that also occur at the SCM after sea ice breakup in central HB (Jacquemot et al., 2022).

The low salinity and warm temperature in the surface waters of southeastern HB compared to HS+UB suggested a strong influence of river runoff in the region at the time of sampling, especially in 2018 ( Figure 1 ; Table 1 ). During our sampling in 2017 and 2018, most of the nutrients at the surface were depleted, with low inorganic nutrient conditions in the surface mixed layer ( Figure 1 ), consistent with the post spring bloom timing of our sample collection. The conditions were favorable for dinoflagellates which accounted for a high proportion of reads at the surface, particularly in the rDNA dataset where they exceeded 70% of the total reads at st732, st731, st730 and st736 ( Figures 3 ; S3A ). The presence of multiple gene copies of the 18S rRNA gene in dinoflagellates might result in an overestimation of their relative abundance in metabarcoding studies (Zhu et al., 2005; Medinger et al., 2010). However, variations in their relative abundance across different regions and depths can offer valuable insights into the preferential ecological niche of dinoflagellates. Compared to HS+UB, the dominance of reads associated with specific ASVs of dinoflagellates in HB suggest that some genera or species were particularly successful despite the low nutrient waters, but consistent with being primarily heterotrophic. For example, the marine stramenopiles from the MAST-1A clade, the Peridiniales Heterocapsa rotundata and the Gymnodiniales Gyrodinium sp. and Gymnodinium sp. had high relative representation in the surface waters of HB compared to the HS+UB ( Figure 4 ). Moreover, both rDNA and rRNA were higher for Heterocapsa and Gymnodinium in the HB than in the HS+UB, suggesting a high potential for growth and acclimatation of these two genera in the warmer waters of the HB surface ( Figure S3A ). Arctic Heterocapsa spp. have chloroplasts (Rintala et al., 2010) but the genus is considered mixotrophic (Millette et al., 2017) and able to photosynthesize and graze on bacteria. While some species of Gymnodiniales are mixotrophic (Li et al., 2000), Gyrodinium spp. are mostly heterotrophic dinoflagellates that graze on small phytoplankton and bacteria (Hoppenrath, 2017). Both genera are nearly always reported in Arctic waters in summer (Lovejoy et al., 2006; Kubiszyn and Wiktor, 2016; Onda et al., 2017). The significant correlation of Gymnodiniales, Gymnodinium and Heterocapsa relative abundance with increasing temperature and decreasing salinity suggest that the warmer and fresher waters in the surface mixed layer were favorable to the growth of these dinoflagellates ( Figures S1 , S2 ). This is in keeping with reports of salinity and temperature having the greatest influence on the growth rate and biomass accumulation of diverse dinoflagellate species (Kim et al., 2004; Dhib et al., 2013). In a recent study, Amiraux et al. (2022) showed that the higher concentrations of dinoflagellates in eastern HB resulted in lower but more diverse content of bioactive molecules being transferred to higher trophic levels. Because the open water season is becoming longer and freshwater inputs are expected to increase in HB, our results suggest heterotrophic groups such as MAST-1A and heterotrophic/mixotrophic dinoflagellates from the Order Peridiniales and Gymnodinales may become more prevalent in summer in coastal HB, with potential influence on the local marine food web.

This study provides the first extensive characterization of microbial eukaryotes in coastal Eastern Hudson Bay, Hudson Strait and the western side of Ungava Bay using a metabarcoding approach. We found that the stability of the water column drove the composition of protist assemblages at the surface with a dinoflagellate-dominated community in the stratified Hudson Bay and a community enriched in diatoms in the more well mixed waters of Hudson Strait. These results support the view that increasing stratification, under the influence of climate warming would favor a more heterotrophic/mixotrophic food web in coastal Hudson Bay. We showed that metabarcoding on the 18S rRNA gene is a suitable tool to detect the occurrence of several potential harmful microalgae. This study provides a time sensitive reference to track the potential introduction of new non-indigenous microbial species into the Hudson Bay System. Future efforts should move toward making DNA metabarcoding in near real time a standard practice to track harmful species in coastal Arctic regions in summer.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA627250.

LJ: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Visualization, Writing – original draft, Writing – review & editing. J-ET: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Validation, Writing – review & editing, Data curation. CM: Data curation, Formal analysis, Methodology, Validation, Writing – review & editing. CL: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.